bEqn.H
Go to the documentation of this file.
1 if (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(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 "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)
56  ==
57  fvModels.source(rho, b)
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 
71 
72  bEqn.solve();
73 
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(up.dimensions(), 1e-8)
95  )
96  );
97 
98  // volScalarField l = 0.337*k*sqrt(k)/epsilon;
99  // Reta *= max((l - dimensionedScalar(dimLength, 1.5e-3))/l, 0);
100 
101  // Calculate Xi flux
102  // ~~~~~~~~~~~~~~~~~
103  surfaceScalarField phiXi
104  (
105  phiSt
107  (
108  fvc::laplacian(thermophysicalTransport->alphaEff(), b)/mgb
109  )*nf
110  + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
111  );
112 
113 
114  // Calculate mean and turbulent strain rates
115  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
116 
117  volVectorField Ut(U + Su*Xi*n);
118  volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
119 
120  volScalarField sigmas
121  (
122  ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
123  + (
124  (n & n)*fvc::div(Su*n)
125  - (n & fvc::grad(Su*n) & n)
126  )*(Xi + scalar(1))/(2*Xi)
127  );
128 
129 
130  // Calculate the unstrained laminar flame speed
131  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
133 
134 
135  // Calculate the laminar flame speed in equilibrium with the applied strain
136  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137  volScalarField SuInf(Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01)));
138 
139  if (SuModel == "unstrained")
140  {
141  Su == Su0;
142  }
143  else if (SuModel == "equilibrium")
144  {
145  Su == SuInf;
146  }
147  else if (SuModel == "transport")
148  {
149  // Solve for the strained laminar flame speed
150  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
151 
152  volScalarField Rc
153  (
154  (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
155  /(sqr(Su0 - SuInf) + sqr(SuMin))
156  );
157 
158  fvScalarMatrix SuEqn
159  (
160  fvm::ddt(rho, Su)
161  + fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
162  - fvm::Sp(fvc::div(phiXi), Su)
163  ==
164  - fvm::SuSp(-rho*Rc*Su0/Su, Su)
165  - fvm::SuSp(rho*(sigmas + Rc), Su)
166  + fvModels.source(rho, Su)
167  );
168 
169  SuEqn.relax();
170 
171  fvConstraints.constrain(SuEqn);
172 
173  SuEqn.solve();
174 
176 
177  // Limit the maximum Su
178  // ~~~~~~~~~~~~~~~~~~~~
179  Su.min(SuMax);
180  Su.max(SuMin);
181  }
182  else
183  {
184  FatalError
185  << args.executable() << " : Unknown Su model " << SuModel
186  << abort(FatalError);
187  }
188 
189 
190  // Calculate Xi according to the selected flame wrinkling model
191  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
192 
193  if (XiModel == "fixed")
194  {
195  // Do nothing, Xi is fixed!
196  }
197  else if (XiModel == "algebraic")
198  {
199  // Simple algebraic model for Xi based on Gulders correlation
200  // with a linear correction function to give a plausible profile for Xi
201  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202  Xi == scalar(1) +
203  (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
204  *XiCoef*sqrt(up/(Su + SuMin))*Reta;
205  }
206  else if (XiModel == "transport")
207  {
208  // Calculate Xi transport coefficients based on Gulders correlation
209  // and DNS data for the rate of generation
210  // with a linear correction function to give a plausible profile for Xi
211  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212 
213  volScalarField XiEqStar
214  (
215  scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
216  );
217 
218  volScalarField XiEq
219  (
220  scalar(1.001)
221  + (
222  scalar(1)
223  + (2*XiShapeCoef)
224  *(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
225  )*(XiEqStar - scalar(1.001))
226  );
227 
228  volScalarField Gstar(0.28/tauEta);
229  volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
230  volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
231 
232  // R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
233 
234  // Solve for the flame wrinkling
235  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
236  fvScalarMatrix XiEqn
237  (
238  fvm::ddt(rho, Xi)
239  + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
240  - fvm::Sp(fvc::div(phiXi), Xi)
241  ==
242  rho*R
243  - fvm::Sp(rho*(R - G), Xi)
244  - fvm::Sp
245  (
246  rho*max
247  (
248  sigmat - sigmas,
249  dimensionedScalar(sigmat.dimensions(), 0)
250  ),
251  Xi
252  )
253  + fvModels.source(rho, Xi)
254  );
255 
256  XiEqn.relax();
257 
258  fvConstraints.constrain(XiEqn);
259 
260  XiEqn.solve();
261 
263 
264  // Correct boundedness of Xi
265  // ~~~~~~~~~~~~~~~~~~~~~~~~~
266  Xi.max(1.0);
267  Info<< "max(Xi) = " << max(Xi).value() << endl;
268  Info<< "max(XiEq) = " << max(XiEq).value() << endl;
269  }
270  else
271  {
272  FatalError
273  << args.executable() << " : Unknown Xi model " << XiModel
274  << abort(FatalError);
275  }
276 
277  Info<< "Combustion progress = "
278  << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
279  << endl;
280 
281  St = Xi*Su;
282 }
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
const word & executable() const
Name of executable without the path.
Definition: argListI.H:36
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
tmp< GeometricField< Type, fvPatchField, volMesh > > SuSp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:102
fluidReactionThermo & thermo
Definition: createFields.H:28
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
engineTime & runTime
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const dimensionedScalar G
Newtonian constant of gravitation.
dimensionedScalar StCorr("StCorr", dimless, 1.0)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
const dimensionedScalar c
Speed of light in a vacuum.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
tmp< fv::convectionScheme< scalar > > mvConvection(fv::convectionScheme< scalar >::New(mesh, fields, phi, mesh.divScheme("div(phi,ft_b_ha_hau)")))
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Info<< "Creating field dpdt\"<< endl;volScalarField dpdt(IOobject("dpdt", runTime.timeName(), mesh), mesh, dimensionedScalar(p.dimensions()/dimTime, 0));Info<< "Creating field kinetic energy K\"<< endl;volScalarField K("K", 0.5 *magSqr(U));Info<< "Creating the unstrained laminar flame speed\"<< endl;autoPtr< laminarFlameSpeed > unstrainedLaminarFlameSpeed(laminarFlameSpeed::New(thermo))
Definition: createFields.H:90
dynamicFvMesh & mesh
Foam::fvConstraints & fvConstraints
phi
Definition: correctPhi.H:3
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\"<< endl;autoPtr< compressible::momentumTransportModel > turbulence(compressible::momentumTransportModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
fluidReactionThermophysicalTransportModel & thermophysicalTransport
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
const tmp< volScalarField::Internal > & Sp
Definition: alphaSuSp.H:7
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
U
Definition: pEqn.H:72
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
scalar epsilon
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
const tmp< volScalarField::Internal > & Su
Definition: alphaSuSp.H:6
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Foam::argList args(argc, argv)
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45