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         solve
 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         // Limit the maximum Su
 148         // ~~~~~~~~~~~~~~~~~~~~
 149         Su.min(SuMax);
 150         Su.max(SuMin);
 151     }
 152     else
 153     {
 154         FatalError
 155             << args.executable() << " : Unknown Su model " << SuModel
 156             << abort(FatalError);
 157     }
 158
 159
 160     // Calculate Xi according to the selected flame wrinkling model
 161     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 162
 163     if (XiModel == "fixed")
 164     {
 165         // Do nothing, Xi is fixed!
 166     }
 167     else if (XiModel == "algebraic")
 168     {
 169         // Simple algebraic model for Xi based on Gulders correlation
 170         // with a linear correction function to give a plausible profile for Xi
 171         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 172         Xi == scalar(1) +
 173             (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
 174            *XiCoef*sqrt(up/(Su + SuMin))*Reta;
 175     }
 176     else if (XiModel == "transport")
 177     {
 178         // Calculate Xi transport coefficients based on Gulders correlation
 179         // and DNS data for the rate of generation
 180         // with a linear correction function to give a plausible profile for Xi
 181         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 182
 183         volScalarField XiEqStar =
 184             scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
 185
 186         volScalarField XiEq =
 187             scalar(1.001)
 188           + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
 189            *(XiEqStar - scalar(1.001));
 190
 191         volScalarField Gstar = 0.28/tauEta;
 192         volScalarField R = Gstar*XiEqStar/(XiEqStar - scalar(1));
 193         volScalarField G = R*(XiEq - scalar(1.001))/XiEq;
 194
 195         //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
 196
 197         // Solve for the flame wrinkling
 198         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 199         solve
 200         (
 201             fvm::ddt(rho, Xi)
 202           + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
 203           - fvm::Sp(fvc::div(phiXi), Xi)
 204          ==
 205             rho*R
 206           - fvm::Sp(rho*(R - G), Xi)
 207           - fvm::Sp
 208             (
 209                 rho*max
 210                 (
 211                     sigmat - sigmas,
 212                     dimensionedScalar("0", sigmat.dimensions(), 0)
 213                 ),
 214                 Xi
 215             )
 216         );
 217
 218
 219         // Correct boundedness of Xi
 220         // ~~~~~~~~~~~~~~~~~~~~~~~~~
 221         Xi.max(1.0);
 222         Info<< "max(Xi) = " << max(Xi).value() << endl;
 223         Info<< "max(XiEq) = " << max(XiEq).value() << endl;
 224     }
 225     else
 226     {
 227         FatalError
 228             << args.executable() << " : Unknown Xi model " << XiModel
 229             << abort(FatalError);
 230     }
 231
 232     Info<< "Combustion progress = "
 233         << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
 234         << endl;
 235
 236     St = Xi*Su;
 237 }