92 lines
3.1 KiB
Text
92 lines
3.1 KiB
Text
|
// Calculate Xi according to the selected flame wrinkling model
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
|
||
|
|
// Calculate coefficients for Gulder's flame speed correlation
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
|
||
|
|
volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*k);
|
||
|
|
volScalarField epsilonPlus = pow(uPrimeCoef, 3)*epsilon;
|
||
|
|
|
||
|
|
volScalarField tauEta = sqrt(mag(thermo->muu()/(rhou*epsilonPlus)));
|
||
|
|
volScalarField Reta = up/
|
||
|
|
(
|
||
|
|
sqrt(epsilonPlus*tauEta)
|
||
|
|
+ dimensionedScalar("1e-8", up.dimensions(), 1e-8)
|
||
|
|
);
|
||
|
|
|
||
|
|
else if (XiModel == "algebraic")
|
||
|
|
{
|
||
|
|
// Simple algebraic model for Xi based on Gulders correlation
|
||
|
|
// with a linear correction function to give a plausible profile for Xi
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
volScalarField GEta = GEtaCoef/tauEta;
|
||
|
|
volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
|
||
|
|
|
||
|
|
volScalarField R =
|
||
|
|
GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
|
||
|
|
|
||
|
|
volScalarField XiEqStar = R/(R - GEta - GIn);
|
||
|
|
|
||
|
|
//- Calculate the unweighted b
|
||
|
|
//volScalarField Rrho = thermo->rhou()/thermo->rhob();
|
||
|
|
//volScalarField bbar = b/(b + (Rrho*(1.0 - b)));
|
||
|
|
|
||
|
|
Xi == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
|
||
|
|
}
|
||
|
|
else if (XiModel == "transport")
|
||
|
|
{
|
||
|
|
// Calculate Xi transport coefficients based on Gulders correlation
|
||
|
|
// and DNS data for the rate of generation
|
||
|
|
// with a linear correction function to give a plausible profile for Xi
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
volScalarField GEta = GEtaCoef/tauEta;
|
||
|
|
volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
|
||
|
|
|
||
|
|
volScalarField R =
|
||
|
|
GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
|
||
|
|
|
||
|
|
volScalarField XiEqStar = R/(R - GEta - GIn);
|
||
|
|
|
||
|
|
volScalarField XiEq =
|
||
|
|
1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
|
||
|
|
|
||
|
|
volScalarField G = R*(XiEq - 1.0)/XiEq;
|
||
|
|
|
||
|
|
// Calculate Xi flux
|
||
|
|
// ~~~~~~~~~~~~~~~~~
|
||
|
|
surfaceScalarField phiXi =
|
||
|
|
phiSt
|
||
|
|
+ (
|
||
|
|
- fvc::interpolate(fvc::laplacian(Dbf, b)/mgb)*nf
|
||
|
|
+ fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
|
||
|
|
);
|
||
|
|
|
||
|
|
|
||
|
|
// Solve for the flame wrinkling
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
solve
|
||
|
|
(
|
||
|
|
betav*fvm::ddt(rho, Xi)
|
||
|
|
+ mvConvection->fvmDiv(phi, Xi)
|
||
|
|
+ fvm::div(phiXi, Xi, "div(phiXi,Xi)")
|
||
|
|
- fvm::Sp(fvc::div(phiXi), Xi)
|
||
|
|
==
|
||
|
|
betav*rho*R
|
||
|
|
- fvm::Sp(betav*rho*(R - G), Xi)
|
||
|
|
);
|
||
|
|
|
||
|
|
|
||
|
|
// Correct boundedness of Xi
|
||
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
|
Xi.max(1.0);
|
||
|
|
Xi = min(Xi, 2.0*XiEq);
|
||
|
|
Info<< "max(Xi) = " << max(Xi).value() << endl;
|
||
|
|
Info<< "max(XiEq) = " << max(XiEq).value() << endl;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
FatalError
|
||
|
|
<< args.executable() << " : Unknown Xi model " << XiModel
|
||
|
|
<< abort(FatalError);
|
||
|
|
}
|