Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2013-08-12 12:03:08 +01:00
17 changed files with 995 additions and 15 deletions

View File

@ -5,7 +5,7 @@
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, 0.0)));
surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
tmp<surfaceScalarField> pPrimeByA;

View File

@ -349,7 +349,7 @@ Foam::RASModels::kineticTheoryModel::divDevRhoReff
void Foam::RASModels::kineticTheoryModel::correct()
{
// Local references
volScalarField alpha(max(this->alpha_, 0.0));
volScalarField alpha(max(this->alpha_, scalar(0)));
const volScalarField& rho = phase_.rho();
const surfaceScalarField& alphaPhi = this->alphaPhi_;
const volVectorField& U = this->U_;

View File

@ -88,7 +88,7 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
{
volScalarField aByaMax
(
cbrt(min(max(alpha, 1e-3), alphaMinFriction)/alphaMax)
cbrt(min(max(alpha, scalar(1e-3)), alphaMinFriction)/alphaMax)
);
return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax));

View File

@ -52,6 +52,9 @@ makeBaseTurbulenceModel
#include "kEpsilon.H"
makeRASModel(kEpsilon);
#include "mixtureKEpsilon.H"
makeRASModel(mixtureKEpsilon);
#include "LaheyKEpsilon.H"
makeRASModel(LaheyKEpsilon);

View File

@ -106,7 +106,7 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::readIfPresent()
(
"GeometricField<Type, PatchField, GeoMesh>::readIfPresent()"
) << "read option IOobject::MUST_READ or MUST_READ_IF_MODIFIED"
<< "suggests that a read constructor for field " << this->name()
<< " suggests that a read constructor for field " << this->name()
<< " would be more appropriate." << endl;
}
else if (this->readOpt() == IOobject::READ_IF_PRESENT && this->headerOk())

View File

@ -200,7 +200,7 @@ NicenoKEqn<BasicTurbulenceModel>::phaseTransferCoeff() const
return
(
max(alphaInversion_ - alpha, 0.0)
max(alphaInversion_ - alpha, scalar(0))
*rho
*min
(

View File

@ -137,7 +137,7 @@ continuousGasKEqn<BasicTurbulenceModel>::phaseTransferCoeff() const
return
(
max(alphaInversion_ - alpha, 0.0)
max(alphaInversion_ - alpha, scalar(0))
*rho
*min
(

View File

@ -83,6 +83,16 @@ LaheyKEpsilon<BasicTurbulenceModel>::LaheyKEpsilon
)
),
C3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cp",
this->coeffDict_,
this->C2_.value()
)
),
Cmub_
(
dimensioned<scalar>::lookupOrAddToDict
@ -110,6 +120,7 @@ bool LaheyKEpsilon<BasicTurbulenceModel>::read()
{
alphaInversion_.readIfPresent(this->coeffDict());
Cp_.readIfPresent(this->coeffDict());
C3_.readIfPresent(this->coeffDict());
Cmub_.readIfPresent(this->coeffDict());
return true;
@ -207,7 +218,7 @@ LaheyKEpsilon<BasicTurbulenceModel>::phaseTransferCoeff() const
return
(
max(alphaInversion_ - alpha, 0.0)
max(alphaInversion_ - alpha, scalar(0))
*rho
*min(gasTurbulence.epsilon()/gasTurbulence.k(), 1.0/U.time().deltaT())
);
@ -244,7 +255,7 @@ tmp<fvScalarMatrix> LaheyKEpsilon<BasicTurbulenceModel>::epsilonSource() const
const volScalarField phaseTransferCoeff(this->phaseTransferCoeff());
return
alpha*rho*this->C2_*this->epsilon_*bubbleG()/this->k_
alpha*rho*this->C3_*this->epsilon_*bubbleG()/this->k_
+ phaseTransferCoeff*gasTurbulence.epsilon()
- fvm::Sp(phaseTransferCoeff, this->epsilon_);
}

View File

@ -110,6 +110,7 @@ protected:
dimensionedScalar alphaInversion_;
dimensionedScalar Cp_;
dimensionedScalar C3_;
dimensionedScalar Cmub_;

View File

@ -123,7 +123,7 @@ void continuousGasKEpsilon<BasicTurbulenceModel>::correctNut()
volScalarField thetal(liquidTurbulence.k()/liquidTurbulence.epsilon());
volScalarField thetag((1.0/(18*liquid.nu()))*sqr(gas.d()));
volScalarField expThetar(exp(min(thetal/thetag, 50.0)));
volScalarField expThetar(exp(min(thetal/thetag, scalar(50))));
volScalarField omega(sqr(expThetar - 1)/(sqr(expThetar) - 1));
nutEff_ = omega*liquidTurbulence.nut();
@ -163,7 +163,15 @@ continuousGasKEpsilon<BasicTurbulenceModel>::nuEff() const
{
volScalarField blend
(
max(min((this->alpha_ - 0.5)/(alphaInversion_ - 0.5), 1.0), 0.0)
max
(
min
(
(this->alpha_ - scalar(0.5))/(alphaInversion_ - 0.5),
scalar(1)
),
scalar(0)
)
);
return tmp<volScalarField>
@ -210,7 +218,7 @@ continuousGasKEpsilon<BasicTurbulenceModel>::phaseTransferCoeff() const
return
(
max(alphaInversion_ - alpha, 0.0)
max(alphaInversion_ - alpha, scalar(0))
*rho
*min
(

View File

@ -0,0 +1,698 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "mixtureKEpsilon.H"
#include "bound.H"
#include "twoPhaseSystem.H"
#include "fixedValueFvPatchFields.H"
#include "inletOutletFvPatchFields.H"
#include "fvmSup.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
mixtureKEpsilon<BasicTurbulenceModel>::mixtureKEpsilon
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
eddyViscosity<RASModel<BasicTurbulenceModel> >
(
type,
alpha,
rho,
U,
alphaPhi,
phi,
transport,
propertiesName
),
liquidTurbulencePtr_(NULL),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
this->coeffDict_,
0.09
)
),
C1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C1",
this->coeffDict_,
1.44
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
this->coeffDict_,
1.92
)
),
C3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C3",
this->coeffDict_,
C2_.value()
)
),
Cp_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cp",
this->coeffDict_,
0.25
)
),
sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmak",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.3
)
),
k_
(
IOobject
(
IOobject::groupName("k", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
epsilon_
(
IOobject
(
IOobject::groupName("epsilon", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
)
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
}
template<class BasicTurbulenceModel>
wordList mixtureKEpsilon<BasicTurbulenceModel>::epsilonBoundaryTypes
(
const volScalarField& epsilon
) const
{
const volScalarField::GeometricBoundaryField& ebf = epsilon.boundaryField();
wordList ebt = ebf.types();
forAll(ebf, patchi)
{
if (isA<fixedValueFvPatchScalarField>(ebf[patchi]))
{
ebt[patchi] = fixedValueFvPatchScalarField::typeName;
}
}
return ebt;
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correctInletOutlet
(
volScalarField& vsf,
const volScalarField& refVsf
) const
{
volScalarField::GeometricBoundaryField& bf = vsf.boundaryField();
const volScalarField::GeometricBoundaryField& refBf =
refVsf.boundaryField();
forAll(bf, patchi)
{
if
(
isA<inletOutletFvPatchScalarField>(bf[patchi])
&& isA<inletOutletFvPatchScalarField>(refBf[patchi])
)
{
refCast<inletOutletFvPatchScalarField>
(bf[patchi]).refValue() =
refCast<const inletOutletFvPatchScalarField>
(refBf[patchi]).refValue();
}
}
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::initMixtureFields()
{
if (rhom_.valid()) return;
// Local references to gas-phase properties
const volScalarField& kg = this->k_;
const volScalarField& epsilong = this->epsilon_;
// Local references to liquid-phase properties
mixtureKEpsilon<BasicTurbulenceModel>& turbc = this->liquidTurbulence();
const volScalarField& kl = turbc.k_;
const volScalarField& epsilonl = turbc.epsilon_;
Ct2_.set
(
new volScalarField
(
IOobject
(
"Ct2",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
Ct2()
)
);
rhom_.set
(
new volScalarField
(
IOobject
(
"rhom",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
rhom()
)
);
km_.set
(
new volScalarField
(
IOobject
(
"km",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mix(kl, kg),
kl.boundaryField().types()
)
);
correctInletOutlet(km_(), kl);
epsilonm_.set
(
new volScalarField
(
IOobject
(
"epsilonm",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mix(epsilonl, epsilong),
epsilonBoundaryTypes(epsilonl)
)
);
correctInletOutlet(epsilonm_(), epsilonl);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool mixtureKEpsilon<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel> >::read())
{
Cmu_.readIfPresent(this->coeffDict());
C1_.readIfPresent(this->coeffDict());
C2_.readIfPresent(this->coeffDict());
C3_.readIfPresent(this->coeffDict());
Cp_.readIfPresent(this->coeffDict());
sigmak_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correctNut()
{
this->nut_ = Cmu_*sqr(k_)/epsilon_;
this->nut_.correctBoundaryConditions();
}
template<class BasicTurbulenceModel>
mixtureKEpsilon<BasicTurbulenceModel>&
mixtureKEpsilon<BasicTurbulenceModel>::liquidTurbulence() const
{
if (!liquidTurbulencePtr_)
{
const volVectorField& U = this->U_;
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
liquidTurbulencePtr_ =
&const_cast<mixtureKEpsilon<BasicTurbulenceModel>&>
(
U.db().lookupObject<mixtureKEpsilon<BasicTurbulenceModel> >
(
IOobject::groupName
(
turbulenceModel::propertiesName,
liquid.name()
)
)
);
}
return *liquidTurbulencePtr_;
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::Ct2() const
{
const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
const volScalarField& alphag = this->alpha_;
volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
volScalarField beta
(
(6*this->Cmu_/(4*sqrt(3.0/2.0)))
*alphag*fluid.drag(gas).K(magUr)/liquid.rho()
*(liquidTurbulence.k_/liquidTurbulence.epsilon_)
);
volScalarField Ct0((3 + beta)/(1 + beta + 2*gas.rho()/liquid.rho()));
volScalarField fAlphad((180 + (-4.71e3 + 4.26e4*alphag)*alphag)*alphag);
return sqr(1 + (Ct0 - 1)*exp(-fAlphad));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rholEff() const
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
return fluid.otherPhase(gas).rho();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhogEff() const
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
return
gas.rho()
+ fluid.Cvm()*fluid.otherPhase(gas).rho();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhom() const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return alphal*rholEff() + alphag*rhogEff();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mix
(
const volScalarField& fc,
const volScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return (alphal*rholEff()*fc + alphag*rhogEff()*fd)/rhom_();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixU
(
const volScalarField& fc,
const volScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return
(alphal*rholEff()*fc + alphag*rhogEff()*Ct2_()*fd)
/(alphal*rholEff() + alphag*rhogEff()*Ct2_());
}
template<class BasicTurbulenceModel>
tmp<surfaceScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixFlux
(
const surfaceScalarField& fc,
const surfaceScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
surfaceScalarField alphalf(fvc::interpolate(alphal));
surfaceScalarField alphagf(fvc::interpolate(alphag));
surfaceScalarField rholEfff(fvc::interpolate(rholEff()));
surfaceScalarField rhogEfff(fvc::interpolate(rhogEff()));
return
fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_()))
*(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd);
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
{
const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
// Lahey model
tmp<volScalarField> bubbleG
(
Cp_
*sqr(liquid)*liquid.rho()
*(
pow3(magUr)
+ pow(fluid.drag(gas).K(magUr)*gas.d()/liquid.rho(), 3.0/4.0)
*pow(magUr, 9.0/4.0)
)
*gas
/gas.d()
);
// Simple model
// tmp<volScalarField> bubbleG
// (
// Cp_*sqr(liquid)*gas*fluid.drag(gas).K(magUr)*sqr(magUr)
// );
return bubbleG;
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::kSource() const
{
return fvm::Su(bubbleG(), km_());
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::epsilonSource() const
{
return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_());
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correct()
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
// Only solve the mixture turbulence for the gas-phase
if (&gas != &fluid.phase1())
{
// This is the liquid phase but check the model for the gas-phase
// is consistent
this->liquidTurbulence();
return;
}
if (!this->turbulence_)
{
return;
}
// Initialise the mixture fields if they have not yet been constructed
initMixtureFields();
// Local references to gas-phase properties
const surfaceScalarField& phig = this->phi_;
const volVectorField& Ug = this->U_;
const volScalarField& alphag = this->alpha_;
volScalarField& kg = this->k_;
volScalarField& epsilong = this->epsilon_;
volScalarField& nutg = this->nut_;
// Local references to liquid-phase properties
mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const surfaceScalarField& phil = liquidTurbulence.phi_;
const volVectorField& Ul = liquidTurbulence.U_;
const volScalarField& alphal = liquidTurbulence.alpha_;
volScalarField& kl = liquidTurbulence.k_;
volScalarField& epsilonl = liquidTurbulence.epsilon_;
volScalarField& nutl = liquidTurbulence.nut_;
// Local references to mixture properties
volScalarField& rhom = rhom_();
volScalarField& km = km_();
volScalarField& epsilonm = epsilonm_();
eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();
// Update the effective mixture density
rhom = this->rhom();
// Mixture flux
surfaceScalarField phim("phim", mixFlux(phil, phig));
// Mixture velocity divergence
volScalarField divUm
(
mixU
(
fvc::div(fvc::absolute(phil, Ul)),
fvc::div(fvc::absolute(phig, Ug))
)
);
tmp<volScalarField> Gc;
{
tmp<volTensorField> tgradUl = fvc::grad(Ul);
Gc = tmp<volScalarField>
(
new volScalarField
(
this->GName(),
nutl*(tgradUl() && dev(twoSymm(tgradUl())))
)
);
tgradUl.clear();
// Update k, epsilon and G at the wall
kl.boundaryField().updateCoeffs();
epsilonl.boundaryField().updateCoeffs();
}
tmp<volScalarField> Gd;
{
tmp<volTensorField> tgradUg = fvc::grad(Ug);
Gd = tmp<volScalarField>
(
new volScalarField
(
this->GName(),
nutg*(tgradUg() && dev(twoSymm(tgradUg())))
)
);
tgradUg.clear();
// Update k, epsilon and G at the wall
kg.boundaryField().updateCoeffs();
epsilong.boundaryField().updateCoeffs();
}
// Mixture turbulence generation
volScalarField Gm(mix(Gc, Gd));
// Mixture turbulence viscosity
volScalarField nutm(mixU(nutl, nutg));
// Update the mixture k and epsilon boundary conditions
km == mix(kl, kg);
bound(km, this->kMin_);
epsilonm == mix(epsilonl, epsilong);
bound(epsilonm, this->epsilonMin_);
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(rhom, epsilonm)
+ fvm::div(phim, epsilonm)
- fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), epsilonm)
- fvm::laplacian(DepsilonEff(rhom*nutm), epsilonm)
==
C1_*rhom*Gm*epsilonm/km
- fvm::SuSp(((2.0/3.0)*C1_)*rhom*divUm, epsilonm)
- fvm::Sp(C2_*rhom*epsilonm/km, epsilonm)
+ epsilonSource()
);
epsEqn().relax();
epsEqn().boundaryManipulate(epsilonm.boundaryField());
solve(epsEqn);
bound(epsilonm, this->epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kmEqn
(
fvm::ddt(rhom, km)
+ fvm::div(phim, km)
- fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), km)
- fvm::laplacian(DkEff(rhom*nutm), km)
==
rhom*Gm
- fvm::SuSp((2.0/3.0)*rhom*divUm, km)
- fvm::Sp(rhom*epsilonm/km, km)
+ kSource()
);
kmEqn().relax();
solve(kmEqn);
bound(km, this->kMin_);
km.correctBoundaryConditions();
volScalarField Cc2(rhom/(alphal*rholEff() + alphag*rhogEff()*Ct2_()));
kl = Cc2*km;
kl.correctBoundaryConditions();
epsilonl = Cc2*epsilonm;
epsilonl.correctBoundaryConditions();
liquidTurbulence.correctNut();
Ct2_() = Ct2();
kg = Ct2_()*kl;
kg.correctBoundaryConditions();
epsilong = Ct2_()*epsilonl;
epsilong.correctBoundaryConditions();
nutg = Ct2_()*(liquidTurbulence.nu()/this->nu())*nutl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,252 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::RASModels::mixtureKEpsilon
Group
grpRASTurbulence
Description
Standard k-epsilon turbulence model
The default model coefficients correspond to the following:
\verbatim
mixtureKEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 -0.33;
sigmak 1.0;
sigmaEps 1.3;
}
\endverbatim
SourceFiles
mixtureKEpsilon.C
\*---------------------------------------------------------------------------*/
#ifndef mixtureKEpsilon_H
#define mixtureKEpsilon_H
#include "RASModel.H"
#include "eddyViscosity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class mixtureKEpsilon Declaration
\*---------------------------------------------------------------------------*/
template<class BasicTurbulenceModel>
class mixtureKEpsilon
:
public eddyViscosity<RASModel<BasicTurbulenceModel> >
{
// Private data
mutable mixtureKEpsilon<BasicTurbulenceModel> *liquidTurbulencePtr_;
// Private Member Functions
// Disallow default bitwise copy construct and assignment
mixtureKEpsilon(const mixtureKEpsilon&);
mixtureKEpsilon& operator=(const mixtureKEpsilon&);
//- Return the turbulence model for the other phase
mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence() const;
protected:
// Protected data
// Model coefficients
dimensionedScalar Cmu_;
dimensionedScalar C1_;
dimensionedScalar C2_;
dimensionedScalar C3_;
dimensionedScalar Cp_;
dimensionedScalar sigmak_;
dimensionedScalar sigmaEps_;
// Fields
volScalarField k_;
volScalarField epsilon_;
// Mixture fields
autoPtr<volScalarField> Ct2_;
autoPtr<volScalarField> rhom_;
autoPtr<volScalarField> km_;
autoPtr<volScalarField> epsilonm_;
// Protected Member Functions
wordList epsilonBoundaryTypes(const volScalarField& epsilon) const;
void correctInletOutlet
(
volScalarField& vsf,
const volScalarField& refVsf
) const;
void initMixtureFields();
virtual void correctNut();
tmp<volScalarField> Ct2() const;
tmp<volScalarField> rholEff() const;
tmp<volScalarField> rhogEff() const;
tmp<volScalarField> rhom() const;
tmp<volScalarField> mix
(
const volScalarField& fc,
const volScalarField& fd
) const;
tmp<volScalarField> mixU
(
const volScalarField& fc,
const volScalarField& fd
) const;
tmp<surfaceScalarField> mixFlux
(
const surfaceScalarField& fc,
const surfaceScalarField& fd
) const;
tmp<volScalarField> bubbleG() const;
virtual tmp<fvScalarMatrix> kSource() const;
virtual tmp<fvScalarMatrix> epsilonSource() const;
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& nutm) const
{
return tmp<volScalarField>
(
new volScalarField
(
"DkEff",
nutm/sigmak_
)
);
}
//- Return the effective diffusivity for epsilon
tmp<volScalarField> DepsilonEff(const volScalarField& nutm) const
{
return tmp<volScalarField>
(
new volScalarField
(
"DepsilonEff",
nutm/sigmaEps_
)
);
}
public:
typedef typename BasicTurbulenceModel::alphaField alphaField;
typedef typename BasicTurbulenceModel::rhoField rhoField;
typedef typename BasicTurbulenceModel::transportModel transportModel;
//- Runtime type information
TypeName("mixtureKEpsilon");
// Constructors
//- Construct from components
mixtureKEpsilon
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName = turbulenceModel::propertiesName,
const word& type = typeName
);
//- Destructor
virtual ~mixtureKEpsilon()
{}
// Member Functions
//- Re-read model coefficients if they have changed
virtual bool read();
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return epsilon_;
}
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "mixtureKEpsilon.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -205,7 +205,7 @@ void Foam::patchInjectionBase::setPositionAndCell
const point pf(tri.randomPoint(rnd));
// position perturbed away from face (into domain)
const scalar a = rnd.position(0.1, 0.5);
const scalar a = rnd.position(scalar(0.1), scalar(0.5));
const vector& pc = mesh.cellCentres()[cellOwner];
const vector d = mag(pf - pc)*patchNormal_[faceI];

View File

@ -197,7 +197,7 @@ updateCoeffs()
const vector& myRayId = dom.IRay(rayId).d();
const scalarField& Ir = dom.Qin();
const scalarField& Ir = dom.Qin().boundaryField()[patchI];
forAll(Iw, faceI)
{

View File

@ -19,10 +19,16 @@ simulationType RAS;
RAS
{
RASModel continuousGasKEpsilon;
RASModel mixtureKEpsilon; //continuousGasKEpsilon;
turbulence on;
printCoeffs on;
// mixtureKEpsilonCoeffs
// {
// Cp 1;
// C3 1;
// }
}
// ************************************************************************* //

View File

@ -19,7 +19,7 @@ simulationType RAS;
RAS
{
RASModel LaheyKEpsilon;
RASModel mixtureKEpsilon; //LaheyKEpsilon;
turbulence on;
printCoeffs on;

View File

@ -40,6 +40,7 @@ divSchemes
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(k|epsilon).*\)" Gauss limitedLinear 1;
"div\(phim,(k|epsilon)m\)" Gauss limitedLinear 1;
"div\(\(\(alpha.*nuEff.*\)*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
}