5 volScalarField c = scalar(1) - b;
9 volScalarField rhou = thermo.rhou();
11 // Calculate flame normal etc.
12 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
14 volVectorField n = fvc::grad(b);
16 volScalarField mgb = mag(n);
18 dimensionedScalar dMgb = 1.0e-3*
19 (b*c*mgb)().weightedAverage(mesh.V())
20 /((b*c)().weightedAverage(mesh.V()) + SMALL)
21 + dimensionedScalar("ddMgb", mgb.dimensions(), SMALL);
25 surfaceVectorField SfHat = mesh.Sf()/mesh.magSf();
26 surfaceVectorField nfVec = fvc::interpolate(n);
27 nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));
28 nfVec /= (mag(nfVec) + dMgb);
29 surfaceScalarField nf = (mesh.Sf() & nfVec);
33 # include <engine/StCorr.H>
35 // Calculate turbulent flame speed flux
36 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37 surfaceScalarField phiSt = fvc::interpolate(rhou*StCorr*Su*Xi)*nf;
41 mesh.surfaceInterpolation::deltaCoeffs()
42 *mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())
43 ).value()*runTime.deltaT().value();
45 Info<< "Max St-Courant Number = " << StCoNum << endl;
52 + mvConvection->fvmDiv(phi, b)
53 + fvm::div(phiSt, b, "div(phiSt,bprog)")
54 - fvm::Sp(fvc::div(phiSt), b)
55 - fvm::laplacian(turbulence->alphaEff(), b)
59 // Add ignition cell contribution to b-equation
60 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
61 # include <engine/ignite.H>
68 Info<< "min(b) = " << min(b).value() << endl;
71 // Calculate coefficients for Gulder's flame speed correlation
72 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
74 volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*turbulence->k());
75 //volScalarField up = sqrt(mag(diag(n * n) & diag(turbulence->r())));
77 volScalarField epsilon = pow(uPrimeCoef, 3)*turbulence->epsilon();
79 volScalarField tauEta = sqrt(thermo.muu()/(rhou*epsilon));
81 volScalarField Reta = up/
83 sqrt(epsilon*tauEta) + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
86 //volScalarField l = 0.337*k*sqrt(k)/epsilon;
87 //Reta *= max((l - dimensionedScalar("dl", dimLength, 1.5e-3))/l, 0.0);
91 surfaceScalarField phiXi =
93 - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
94 + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf;
97 // Calculate mean and turbulent strain rates
98 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
100 volVectorField Ut = U + Su*Xi*n;
101 volScalarField sigmat = (n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n);
103 volScalarField sigmas =
104 ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
106 (n & n)*fvc::div(Su*n)
107 - (n & fvc::grad(Su*n) & n)
108 )*(Xi + scalar(1))/(2*Xi);
111 // Calculate the unstrained laminar flame speed
112 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
113 volScalarField Su0 = unstrainedLaminarFlameSpeed()();
116 // Calculate the laminar flame speed in equilibrium with the applied strain
117 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118 volScalarField SuInf = Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01));
120 if (SuModel == "unstrained")
124 else if (SuModel == "equilibrium")
128 else if (SuModel == "transport")
130 // Solve for the strained laminar flame speed
131 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134 (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
135 /(sqr(Su0 - SuInf) + sqr(SuMin));
140 + fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
141 - fvm::Sp(fvc::div(phiXi), Su)
143 - fvm::SuSp(-rho*Rc*Su0/Su, Su)
144 - fvm::SuSp(rho*(sigmas + Rc), Su)
147 // Limit the maximum Su
148 // ~~~~~~~~~~~~~~~~~~~~
155 << args.executable() << " : Unknown Su model " << SuModel
156 << abort(FatalError);
160 // Calculate Xi according to the selected flame wrinkling model
161 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
163 if (XiModel == "fixed")
165 // Do nothing, Xi is fixed!
167 else if (XiModel == "algebraic")
169 // Simple algebraic model for Xi based on Gulders correlation
170 // with a linear correction function to give a plausible profile for Xi
171 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
173 (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
174 *XiCoef*sqrt(up/(Su + SuMin))*Reta;
176 else if (XiModel == "transport")
178 // Calculate Xi transport coefficients based on Gulders correlation
179 // and DNS data for the rate of generation
180 // with a linear correction function to give a plausible profile for Xi
181 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
183 volScalarField XiEqStar =
184 scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
186 volScalarField XiEq =
188 + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
189 *(XiEqStar - scalar(1.001));
191 volScalarField Gstar = 0.28/tauEta;
192 volScalarField R = Gstar*XiEqStar/(XiEqStar - scalar(1));
193 volScalarField G = R*(XiEq - scalar(1.001))/XiEq;
195 //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
197 // Solve for the flame wrinkling
198 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202 + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
203 - fvm::Sp(fvc::div(phiXi), Xi)
206 - fvm::Sp(rho*(R - G), Xi)
212 dimensionedScalar("0", sigmat.dimensions(), 0)
219 // Correct boundedness of Xi
220 // ~~~~~~~~~~~~~~~~~~~~~~~~~
222 Info<< "max(Xi) = " << max(Xi).value() << endl;
223 Info<< "max(XiEq) = " << max(XiEq).value() << endl;
228 << args.executable() << " : Unknown Xi model " << XiModel
229 << abort(FatalError);
232 Info<< "Combustion progress = "
233 << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
239 // ************************ vim: set sw=4 sts=4 et: ************************ //