bEqn.H
Go to the documentation of this file.
1if (ign.ignited())
2{
3 // progress variable
4 // ~~~~~~~~~~~~~~~~~
5 volScalarField c("c", scalar(1) - b);
6
7 // Unburnt gas density
8 // ~~~~~~~~~~~~~~~~~~~
9 volScalarField rhou(thermo.rhou());
10
11 // Calculate flame normal etc.
12 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
13
14 volVectorField n("n", fvc::grad(b));
15
16 volScalarField mgb(mag(n));
17
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);
22
23 mgb += dMgb;
24
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));
30 n /= mgb;
31
32
33 #include "StCorr.H"
34
35 // Calculate turbulent flame speed flux
36 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37 surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
38
39 scalar StCoNum = max
40 (
41 mesh.surfaceInterpolation::deltaCoeffs()
42 *mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())
43 ).value()*runTime.deltaTValue();
44
45 Info<< "Max St-Courant Number = " << StCoNum << endl;
46
47 // Create b equation
48 // ~~~~~~~~~~~~~~~~~
49 fvScalarMatrix bEqn
50 (
51 fvm::ddt(rho, b)
52 + mvConvection->fvmDiv(phi, b)
53 + fvm::div(phiSt, b)
54 - fvm::Sp(fvc::div(phiSt), b)
55 - fvm::laplacian(turbulence->alphaEff(), b)
56 ==
58 );
59
60
61 // Add ignition cell contribution to b-equation
62 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
63 #include "ignite.H"
64
65
66 // Solve for b
67 // ~~~~~~~~~~~
68 bEqn.relax();
69
70 fvOptions.constrain(bEqn);
71
72 bEqn.solve();
73
74 fvOptions.correct(b);
75
76 Info<< "min(b) = " << min(b).value() << endl;
77
78
79 // Calculate coefficients for Gulder's flame speed correlation
80 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
81
82 volScalarField up(uPrimeCoef*sqrt((2.0/3.0)*turbulence->k()));
83 //volScalarField up(sqrt(mag(diag(n * n) & diag(turbulence->r()))));
84
85 volScalarField epsilon(pow(uPrimeCoef, 3)*turbulence->epsilon());
86
87 volScalarField tauEta(sqrt(thermo.muu()/(rhou*epsilon)));
88
89 volScalarField Reta
90 (
91 up
92 / (
93 sqrt(epsilon*tauEta)
94 + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
95 )
96 );
97
98 //volScalarField l = 0.337*k*sqrt(k)/epsilon;
99 //Reta *= max((l - dimensionedScalar("dl", dimLength, 1.5e-3))/l, 0.0);
100
101 // Calculate Xi flux
102 // ~~~~~~~~~~~~~~~~~
103 surfaceScalarField phiXi
104 (
105 phiSt
106 - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
107 + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
108 );
109
110
111 // Calculate mean and turbulent strain rates
112 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
113
114 volVectorField Ut(U + Su*Xi*n);
115 volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
116
117 volScalarField sigmas
118 (
119 ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
120 + (
121 (n & n)*fvc::div(Su*n)
122 - (n & fvc::grad(Su*n) & n)
123 )*(Xi + scalar(1))/(2*Xi)
124 );
125
126
127 // Calculate the unstrained laminar flame speed
128 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129 volScalarField Su0(unstrainedLaminarFlameSpeed()());
130
131
132 // Calculate the laminar flame speed in equilibrium with the applied strain
133 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134 volScalarField SuInf(Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01)));
135
136 if (SuModel == "unstrained")
137 {
138 Su == Su0;
139 }
140 else if (SuModel == "equilibrium")
141 {
142 Su == SuInf;
143 }
144 else if (SuModel == "transport")
145 {
146 // Solve for the strained laminar flame speed
147 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
148
149 volScalarField Rc
150 (
151 (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
152 /(sqr(Su0 - SuInf) + sqr(SuMin))
153 );
154
155 fvScalarMatrix SuEqn
156 (
157 fvm::ddt(rho, Su)
158 + fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
159 - fvm::Sp(fvc::div(phiXi), Su)
160 ==
161 - fvm::SuSp(-rho*Rc*Su0/Su, Su)
162 - fvm::SuSp(rho*(sigmas + Rc), Su)
163 + fvOptions(rho, Su)
164 );
165
166 SuEqn.relax();
167
168 fvOptions.constrain(SuEqn);
169
170 SuEqn.solve();
171
172 fvOptions.correct(Su);
173
174 // Limit the maximum Su
175 // ~~~~~~~~~~~~~~~~~~~~
176 Su.min(SuMax);
177 Su.max(SuMin);
178 }
179 else
180 {
181 FatalError
182 << args.executable() << " : Unknown Su model " << SuModel
183 << abort(FatalError);
184 }
185
186
187 // Calculate Xi according to the selected flame wrinkling model
188 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
189
190 if (XiModel == "fixed")
191 {
192 // Do nothing, Xi is fixed!
193 }
194 else if (XiModel == "algebraic")
195 {
196 // Simple algebraic model for Xi based on Gulders correlation
197 // with a linear correction function to give a plausible profile for Xi
198 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
199 Xi == scalar(1) +
200 (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
201 *XiCoef*sqrt(up/(Su + SuMin))*Reta;
202 }
203 else if (XiModel == "transport")
204 {
205 // Calculate Xi transport coefficients based on Gulders correlation
206 // and DNS data for the rate of generation
207 // with a linear correction function to give a plausible profile for Xi
208 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
209
210 volScalarField XiEqStar
211 (
212 scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
213 );
214
215 volScalarField XiEq
216 (
217 scalar(1.001)
218 + (
219 scalar(1)
220 + (2*XiShapeCoef)
221 *(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
222 )*(XiEqStar - scalar(1.001))
223 );
224
225 volScalarField Gstar(0.28/tauEta);
226 volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
227 volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
228
229 //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
230
231 // Solve for the flame wrinkling
232 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233 fvScalarMatrix XiEqn
234 (
235 fvm::ddt(rho, Xi)
236 + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
237 - fvm::Sp(fvc::div(phiXi), Xi)
238 ==
239 rho*R
240 - fvm::Sp(rho*(R - G), Xi)
241 - fvm::Sp
242 (
243 rho*max
244 (
245 sigmat - sigmas,
246 dimensionedScalar(sigmat.dimensions(), Zero)
247 ),
248 Xi
249 )
250 + fvOptions(rho, Xi)
251 );
252
253 XiEqn.relax();
254
255 fvOptions.constrain(XiEqn);
256
257 XiEqn.solve();
258
259 fvOptions.correct(Xi);
260
261 // Correct boundedness of Xi
262 // ~~~~~~~~~~~~~~~~~~~~~~~~~
263 Xi.max(1.0);
264 Info<< "max(Xi) = " << max(Xi).value() << endl;
265 Info<< "max(XiEq) = " << max(XiEq).value() << endl;
266 }
267 else
268 {
269 FatalError
270 << args.executable() << " : Unknown Xi model " << XiModel
271 << abort(FatalError);
272 }
273
274 Info<< "Combustion progress = "
275 << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
276 << endl;
277
278 St = Xi*Su;
279}
tmp< fv::convectionScheme< scalar > > mvConvection(fv::convectionScheme< scalar >::New(mesh, fields, phi, mesh.divScheme("div(phi,ft_b_ha_hau)")))
#define R(A, B, C, D, E, F, K, M)
dimensionedScalar StCorr("StCorr", dimless, 1.0)
label n
Y[inertIndex] max(0.0)
fv::options & fvOptions
surfaceScalarField & phi
const word & executable() const noexcept
Name of executable without the path.
Definition: argListI.H:51
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
U
Definition: pEqn.H:72
dynamicFvMesh & mesh
engineTime & runTime
scalar epsilon
compressible::turbulenceModel & turbulence
zeroField Su
Definition: alphaSuSp.H:1
Foam::argList args(argc, argv)
Info<< "Creating the unstrained laminar flame speed\n"<< endl;autoPtr< laminarFlameSpeed > unstrainedLaminarFlameSpeed(laminarFlameSpeed::New(thermo))
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11