applications/solvers/combustion/PDRFoam/XiEqns

   1     // Calculate Xi according to the selected flame wrinkling model
   2     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   3
   4     // Calculate coefficients for Gulder's flame speed correlation
   5     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   6
   7     volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*k);
   8     volScalarField epsilonPlus = pow(uPrimeCoef, 3)*epsilon;
   9
  10     volScalarField tauEta = sqrt(mag(thermo->muu()/(rhou*epsilonPlus)));
  11     volScalarField Reta = up/
  12     (
  13         sqrt(epsilonPlus*tauEta)
  14       + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
  15     );
  16
  17     else if (XiModel == "algebraic")
  18     {
  19         // Simple algebraic model for Xi based on Gulders correlation
  20         // with a linear correction function to give a plausible profile for Xi
  21         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  22         volScalarField GEta = GEtaCoef/tauEta;
  23         volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
  24
  25         volScalarField R =
  26             GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
  27
  28         volScalarField XiEqStar = R/(R - GEta - GIn);
  29
  30         //- Calculate the unweighted b
  31         //volScalarField Rrho = thermo->rhou()/thermo->rhob();
  32         //volScalarField bbar = b/(b + (Rrho*(1.0 - b)));
  33
  34         Xi == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
  35     }
  36     else if (XiModel == "transport")
  37     {
  38         // Calculate Xi transport coefficients based on Gulders correlation
  39         // and DNS data for the rate of generation
  40         // with a linear correction function to give a plausible profile for Xi
  41         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  42         volScalarField GEta = GEtaCoef/tauEta;
  43         volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
  44
  45         volScalarField R =
  46             GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
  47
  48         volScalarField XiEqStar = R/(R - GEta - GIn);
  49
  50         volScalarField XiEq =
  51             1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
  52
  53         volScalarField G = R*(XiEq - 1.0)/XiEq;
  54
  55     // Calculate Xi flux
  56     // ~~~~~~~~~~~~~~~~~
  57     surfaceScalarField phiXi =
  58         phiSt
  59       + (
  60           - fvc::interpolate(fvc::laplacian(Dbf, b)/mgb)*nf
  61           + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
  62         );
  63
  64
  65         // Solve for the flame wrinkling
  66         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  67         solve
  68         (
  69             betav*fvm::ddt(rho, Xi)
  70           + mvConvection->fvmDiv(phi, Xi)
  71           + fvm::div(phiXi, Xi, "div(phiXi,Xi)")
  72           - fvm::Sp(fvc::div(phiXi), Xi)
  73          ==
  74             betav*rho*R
  75           - fvm::Sp(betav*rho*(R - G), Xi)
  76         );
  77
  78
  79         // Correct boundedness of Xi
  80         // ~~~~~~~~~~~~~~~~~~~~~~~~~
  81         Xi.max(1.0);
  82         Xi = min(Xi, 2.0*XiEq);
  83         Info<< "max(Xi) = " << max(Xi).value() << endl;
  84         Info<< "max(XiEq) = " << max(XiEq).value() << endl;
  85     }
  86     else
  87     {
  88         FatalError
  89             << args.executable() << " : Unknown Xi model " << XiModel
  90             << abort(FatalError);
  91     }