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