mixtureKEpsilon: Changed bounded non-conservative transport flux to be volumetric rather than mass-based

Avoids numerical problems caused by large density gradients, particularly around the interface.
Allows the use of higher-order differencing for km and epsilonm
This commit is contained in:
Henry
2015-03-20 17:16:05 +00:00
parent b6e803ab8b
commit bcb8b30cad

View File

@ -472,8 +472,8 @@ tmp<surfaceScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixFlux
surfaceScalarField rhogEfff(fvc::interpolate(rhogEff()));
return
fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_()))
*(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd);
(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd)
/(alphalf*rholEfff + alphagf*rhogEfff*fvc::interpolate(Ct2_()));
}
@ -516,14 +516,14 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::kSource() const
{
return fvm::Su(bubbleG(), km_());
return fvm::Su(bubbleG()/rhom_(), km_());
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::epsilonSource() const
{
return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_());
return fvm::Su(C3_*epsilonm_()*bubbleG()/(rhom_()*km_()), epsilonm_());
}
@ -647,14 +647,14 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(rhom, epsilonm)
fvm::ddt(epsilonm)
+ fvm::div(phim, epsilonm)
- fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), epsilonm)
- fvm::laplacian(DepsilonEff(rhom*nutm), epsilonm)
- fvm::Sp(fvc::div(phim), epsilonm)
- fvm::laplacian(DepsilonEff(nutm), epsilonm)
==
C1_*rhom*Gm*epsilonm/km
- fvm::SuSp(((2.0/3.0)*C1_)*rhom*divUm, epsilonm)
- fvm::Sp(C2_*rhom*epsilonm/km, epsilonm)
C1_*Gm*epsilonm/km
- fvm::SuSp(((2.0/3.0)*C1_)*divUm, epsilonm)
- fvm::Sp(C2_*epsilonm/km, epsilonm)
+ epsilonSource()
);
@ -669,14 +669,14 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kmEqn
(
fvm::ddt(rhom, km)
fvm::ddt(km)
+ fvm::div(phim, km)
- fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), km)
- fvm::laplacian(DkEff(rhom*nutm), km)
- fvm::Sp(fvc::div(phim), km)
- fvm::laplacian(DkEff(nutm), km)
==
rhom*Gm
- fvm::SuSp((2.0/3.0)*rhom*divUm, km)
- fvm::Sp(rhom*epsilonm/km, km)
Gm
- fvm::SuSp((2.0/3.0)*divUm, km)
- fvm::Sp(epsilonm/km, km)
+ kSource()
);