FreeFOAM The Cross-Platform CFD Toolkit
bEqn.H
Go to the documentation of this file.
1 tmp<fv::convectionScheme<scalar> > mvConvection
2 (
3  fv::convectionScheme<scalar>::New
4  (
5  mesh,
6  fields,
7  phi,
8  mesh.divScheme("div(phi,ft_b_h_hu)")
9  )
10 );
11 
12 volScalarField Db("Db", turbulence->muEff());
13 
14 if (ign.ignited())
15 {
16  // Calculate the unstrained laminar flame speed
17  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
19 
20  const volScalarField& Xi = flameWrinkling->Xi();
21 
22  // progress variable
23  // ~~~~~~~~~~~~~~~~~
24  volScalarField c("c", 1.0 - b);
25 
26  // Unburnt gas density
27  // ~~~~~~~~~~~~~~~~~~~
28  volScalarField rhou = thermo.rhou();
29 
30  // Calculate flame normal etc.
31  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
32 
33  //volVectorField n = fvc::grad(b);
34  volVectorField n = fvc::reconstruct(fvc::snGrad(b,"snGrad(bprog)")*mesh.magSf());
35 
36  volScalarField mgb("mgb", mag(n));
37 
38  dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
39 
40  {
41  volScalarField bc = b*c;
42 
43  dMgb += 1.0e-3*
44  (bc*mgb)().weightedAverage(mesh.V())
45  /(bc.weightedAverage(mesh.V()) + SMALL);
46  }
47 
48  mgb += dMgb;
49 
50  surfaceVectorField Sfhat = mesh.Sf()/mesh.magSf();
52  nfVec += Sfhat*(fvc::snGrad(b,"snGrad(bprog)") - (Sfhat & nfVec));
53  nfVec /= (mag(nfVec) + dMgb);
54  surfaceScalarField nf("nf", mesh.Sf() & nfVec);
55  n /= mgb;
56 
57 # include <engine/StCorr.H>
58 
59  // Calculate turbulent flame speed flux
60  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
61  surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*St)*nf);
62 
63 # include "StCourantNo.H"
64 
65  Db = flameWrinkling->Db();
66 
67  // Create b equation
68  // ~~~~~~~~~~~~~~~~~
69  fvScalarMatrix bEqn
70  (
71  betav*fvm::ddt(rho, b)
72  + mvConvection->fvmDiv(phi, b)
73  + fvm::div(phiSt, b, "div(phiSt,bprog)")
74  - fvm::Sp(fvc::div(phiSt), b)
75  - fvm::laplacian(Db, b, "laplacian(DB,bprog)")
76  );
77 
78 
79  // Add ignition cell contribution to b-equation
80  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
81 # include <engine/ignite.H>
82 
83  // Solve for b
84  // ~~~~~~~~~~~
85  bEqn.solve();
86 
87  Info<< "min(bprog) = " << min(b).value() << endl;
88 
89  if (composition.contains("ft"))
90  {
91  volScalarField& ft = composition.Y("ft");
92 
93  Info<< "Combustion progress = "
94  << 100*(1.0 - b)().weightedAverage(mesh.V()*ft).value() << "%"
95  << endl;
96  }
97  else
98  {
99  Info<< "Combustion progress = "
100  << 100*(1.0 - b)().weightedAverage(mesh.V()).value() << "%"
101  << endl;
102  }
103 
104  // Correct the flame-wrinkling
105  flameWrinkling->correct(mvConvection);
106 
107  St = Xi*Su;
108 }
109 
110 // ************************ vim: set sw=4 sts=4 et: ************************ //