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.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, "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
  82     (
  83         up
  84       / (
  85             sqrt(epsilon*tauEta)
  86           + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
  87         )
  88     );
  89
  90   //volScalarField l = 0.337*k*sqrt(k)/epsilon;
  91   //Reta *= max((l - dimensionedScalar("dl", dimLength, 1.5e-3))/l, 0.0);
  92
  93     // Calculate Xi flux
  94     // ~~~~~~~~~~~~~~~~~
  95     surfaceScalarField phiXi
  96     (
  97         phiSt
  98       - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
  99       + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
 100     );
 101
 102
 103     // Calculate mean and turbulent strain rates
 104     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 105
 106     volVectorField Ut(U + Su*Xi*n);
 107     volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
 108
 109     volScalarField sigmas
 110     (
 111         ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
 112       + (
 113             (n & n)*fvc::div(Su*n)
 114           - (n & fvc::grad(Su*n) & n)
 115         )*(Xi + scalar(1))/(2*Xi)
 116     );
 117
 118
 119     // Calculate the unstrained laminar flame speed
 120     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 121     volScalarField Su0(unstrainedLaminarFlameSpeed()());
 122
 123
 124     // Calculate the laminar flame speed in equilibrium with the applied strain
 125     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 126     volScalarField SuInf(Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01)));
 127
 128     if (SuModel == "unstrained")
 129     {
 130         Su == Su0;
 131     }
 132     else if (SuModel == "equilibrium")
 133     {
 134         Su == SuInf;
 135     }
 136     else if (SuModel == "transport")
 137     {
 138         // Solve for the strained laminar flame speed
 139         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 140
 141         volScalarField Rc
 142         (
 143             (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
 144             /(sqr(Su0 - SuInf) + sqr(SuMin))
 145         );
 146
 147         fvScalarMatrix SuEqn
 148         (
 149             fvm::ddt(rho, Su)
 150           + fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
 151           - fvm::Sp(fvc::div(phiXi), Su)
 152           ==
 153           - fvm::SuSp(-rho*Rc*Su0/Su, Su)
 154           - fvm::SuSp(rho*(sigmas + Rc), Su)
 155         );
 156
 157         SuEqn.relax();
 158         SuEqn.solve();
 159
 160         // Limit the maximum Su
 161         // ~~~~~~~~~~~~~~~~~~~~
 162         Su.min(SuMax);
 163         Su.max(SuMin);
 164     }
 165     else
 166     {
 167         FatalError
 168             << args.executable() << " : Unknown Su model " << SuModel
 169             << abort(FatalError);
 170     }
 171
 172
 173     // Calculate Xi according to the selected flame wrinkling model
 174     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 175
 176     if (XiModel == "fixed")
 177     {
 178         // Do nothing, Xi is fixed!
 179     }
 180     else if (XiModel == "algebraic")
 181     {
 182         // Simple algebraic model for Xi based on Gulders correlation
 183         // with a linear correction function to give a plausible profile for Xi
 184         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 185         Xi == scalar(1) +
 186             (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
 187            *XiCoef*sqrt(up/(Su + SuMin))*Reta;
 188     }
 189     else if (XiModel == "transport")
 190     {
 191         // Calculate Xi transport coefficients based on Gulders correlation
 192         // and DNS data for the rate of generation
 193         // with a linear correction function to give a plausible profile for Xi
 194         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 195
 196         volScalarField XiEqStar
 197         (
 198             scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
 199         );
 200
 201         volScalarField XiEq
 202         (
 203             scalar(1.001)
 204           + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
 205             *(XiEqStar - scalar(1.001))
 206         );
 207
 208         volScalarField Gstar(0.28/tauEta);
 209         volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
 210         volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
 211
 212         //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
 213
 214         // Solve for the flame wrinkling
 215         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 216         fvScalarMatrix XiEqn
 217         (
 218             fvm::ddt(rho, Xi)
 219           + fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
 220           - fvm::Sp(fvc::div(phiXi), Xi)
 221          ==
 222             rho*R
 223           - fvm::Sp(rho*(R - G), Xi)
 224           - fvm::Sp
 225             (
 226                 rho*max
 227                 (
 228                     sigmat - sigmas,
 229                     dimensionedScalar("0", sigmat.dimensions(), 0)
 230                 ),
 231                 Xi
 232             )
 233         );
 234
 235         XiEqn.relax();
 236         XiEqn.solve();
 237
 238         // Correct boundedness of Xi
 239         // ~~~~~~~~~~~~~~~~~~~~~~~~~
 240         Xi.max(1.0);
 241         Info<< "max(Xi) = " << max(Xi).value() << endl;
 242         Info<< "max(XiEq) = " << max(XiEq).value() << endl;
 243     }
 244     else
 245     {
 246         FatalError
 247             << args.executable() << " : Unknown Xi model " << XiModel
 248             << abort(FatalError);
 249     }
 250
 251     Info<< "Combustion progress = "
 252         << 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
 253         << endl;
 254
 255     St = Xi*Su;
 256 }