FreeFOAM The Cross-Platform CFD Toolkit
bEqn.H
Go to the documentation of this file.
1 if (ign.ignited())
2 {
3  // progress variable
4  // ~~~~~~~~~~~~~~~~~
5  volScalarField c = scalar(1) - b;
6 
7  // Unburnt gas density
8  // ~~~~~~~~~~~~~~~~~~~
9  volScalarField rhou = thermo.rhou();
10 
11  // Calculate flame normal etc.
12  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
13 
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();
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 <engine/StCorr.H>
34 
35  // Calculate turbulent flame speed flux
36  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38 
39  scalar StCoNum = max
40  (
41  mesh.surfaceInterpolation::deltaCoeffs()
42  *mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())
43  ).value()*runTime.deltaT().value();
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, "div(phiSt,bprog)")
54  - fvm::Sp(fvc::div(phiSt), b)
55  - fvm::laplacian(turbulence->alphaEff(), b)
56  );
57 
58 
59  // Add ignition cell contribution to b-equation
60  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
61  #include <engine/ignite.H>
62 
63 
64  // Solve for b
65  // ~~~~~~~~~~~
66  bEqn.solve();
67 
68  Info<< "min(bprog) = " << min(b).value() << endl;
69 
70 
71  // Calculate coefficients for Gulder's flame speed correlation
72  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
73 
74  volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*turbulence->k());
75  //volScalarField up = sqrt(mag(diag(n * n) & diag(turbulence->r())));
76 
77  volScalarField epsilon = pow(uPrimeCoef, 3)*turbulence->epsilon();
78 
79  volScalarField tauEta = sqrt(thermo.muu()/(rhou*epsilon));
80 
81  volScalarField Reta = up/
82  (
83  sqrt(epsilon*tauEta) + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
84  );
85 
86  //volScalarField l = 0.337*k*sqrt(k)/epsilon;
87  //Reta *= max((l - dimensionedScalar("dl", dimLength, 1.5e-3))/l, 0.0);
88 
89  // Calculate Xi flux
90  // ~~~~~~~~~~~~~~~~~
91  surfaceScalarField phiXi =
92  phiSt
93  - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
94  + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf;
95 
96 
97  // Calculate mean and turbulent strain rates
98  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99 
100  volVectorField Ut = U + Su*Xi*n;
101  volScalarField sigmat = (n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n);
102 
103  volScalarField sigmas =
104  ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
105  + (
106  (n & n)*fvc::div(Su*n)
107  - (n & fvc::grad(Su*n) & n)
108  )*(Xi + scalar(1))/(2*Xi);
109 
110 
111  // Calculate the unstrained laminar flame speed
112  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114 
115 
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));
119 
120  if (SuModel == "unstrained")
121  {
122  Su == Su0;
123  }
124  else if (SuModel == "equilibrium")
125  {
126  Su == SuInf;
127  }
128  else if (SuModel == "transport")
129  {
130  // Solve for the strained laminar flame speed
131  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132 
133  volScalarField Rc =
134  (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
135  /(sqr(Su0 - SuInf) + sqr(SuMin));
136 
137  fvScalarMatrix SuEqn
138  (
139  fvm::ddt(rho, Su)
140  + fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
141  - fvm::Sp(fvc::div(phiXi), Su)
142  ==
143  - fvm::SuSp(-rho*Rc*Su0/Su, Su)
144  - fvm::SuSp(rho*(sigmas + Rc), Su)
145  );
146 
147  SuEqn.relax();
148  SuEqn.solve();
149 
150  // Limit the maximum Su
151  // ~~~~~~~~~~~~~~~~~~~~
152  Su.min(SuMax);
153  Su.max(SuMin);
154  }
155  else
156  {
157  FatalError
158  << args.executable() << " : Unknown Su model " << SuModel
159  << abort(FatalError);
160  }
161 
162 
163  // Calculate Xi according to the selected flame wrinkling model
164  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
165 
166  if (XiModel == "fixed")
167  {
168  // Do nothing, Xi is fixed!
169  }
170  else if (XiModel == "algebraic")
171  {
172  // Simple algebraic model for Xi based on Gulders correlation
173  // with a linear correction function to give a plausible profile for Xi
174  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175  Xi == scalar(1) +
176  (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
177  *XiCoef*sqrt(up/(Su + SuMin))*Reta;
178  }
179  else if (XiModel == "transport")
180  {
181  // Calculate Xi transport coefficients based on Gulders correlation
182  // and DNS data for the rate of generation
183  // with a linear correction function to give a plausible profile for Xi
184  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185 
186  volScalarField XiEqStar =
187  scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
188 
189  volScalarField XiEq =
190  scalar(1.001)
191  + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
192  *(XiEqStar - scalar(1.001));
193 
194  volScalarField Gstar = 0.28/tauEta;
195  volScalarField R = Gstar*XiEqStar/(XiEqStar - scalar(1));
196  volScalarField G = R*(XiEq - scalar(1.001))/XiEq;
197 
198  //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
199 
200  // Solve for the flame wrinkling
201  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202  fvScalarMatrix XiEqn
203  (
204  fvm::ddt(rho, Xi)
205  + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
206  - fvm::Sp(fvc::div(phiXi), Xi)
207  ==
208  rho*R
209  - fvm::Sp(rho*(R - G), Xi)
210  - fvm::Sp
211  (
212  rho*max
213  (
214  sigmat - sigmas,
215  dimensionedScalar("0", sigmat.dimensions(), 0)
216  ),
217  Xi
218  )
219  );
220 
221  XiEqn.relax();
222  XiEqn.solve();
223 
224  // Correct boundedness of Xi
225  // ~~~~~~~~~~~~~~~~~~~~~~~~~
226  Xi.max(1.0);
227  Info<< "max(Xi) = " << max(Xi).value() << endl;
228  Info<< "max(XiEq) = " << max(XiEq).value() << endl;
229  }
230  else
231  {
232  FatalError
233  << args.executable() << " : Unknown Xi model " << XiModel
234  << abort(FatalError);
235  }
236 
237  Info<< "Combustion progress = "
238  << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
239  << endl;
240 
241  St = Xi*Su;
242 }
243 
244 // ************************ vim: set sw=4 sts=4 et: ************************ //