applications/solvers/combustion/XiFoam/bEqn.H

   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
  14     volVectorField 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 = 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.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,b)")
  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 "ignite.H"
  62
  63
  64     // Solve for b
  65     // ~~~~~~~~~~~
  66     bEqn.solve();
  67
  68     Info<< "min(b) = " << 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     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 113     volScalarField Su0 = unstrainedLaminarFlameSpeed()();
 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 }