This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/solvers/combustion/PDRFoam/XiEqns

91 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);
}