applications/solvers/combustion/PDRFoam/bEqn.H

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