applications/solvers/combustion/PDRFoam/kEpsilon.H

   1 if(turbulence)
   2 {
   3     volScalarField divU = fvc::div(Uf & mesh.Sf());
   4
   5     tmp<volTensorField> tgradU = fvc::grad(Uf);
   6     volScalarField G = 2*mut*(tgradU() && dev(symm(tgradU())));
   7     tgradU.clear();
   8
   9     // Add the blockage generation term so that it is included consistently
  10     // in both the k and epsilon equations
  11     volScalarField GR = rho*mag(U)*(U & CT & U);
  12
  13 #   include "wallFunctions.H"
  14
  15     // Dissipation equation
  16     fvScalarMatrix epsEqn
  17     (
  18         betav*fvm::ddt(rho, epsilon)
  19       + fvm::div(phi, epsilon)
  20       - fvm::laplacian(fvc::interpolate(alphaEps*muEff), epsilon)
  21       ==
  22         C1*(betav*G + GR)*epsilon/k
  23       - fvm::SuSp((2.0/3.0*C1)*betav*rho*divU, epsilon)
  24       - fvm::Sp(C2*betav*rho*epsilon/k, epsilon)
  25     );
  26
  27 #   include "wallDissipation.H"
  28
  29     epsEqn.solve();
  30     bound(epsilon, dimensionedScalar("0", epsilon.dimensions(), 1.0e-15));
  31
  32
  33     // Turbulent kinetic energy equation
  34     solve
  35     (
  36         betav*fvm::ddt(rho, k)
  37       + fvm::div(phi, k)
  38       - fvm::laplacian(fvc::interpolate(alphak*muEff), k)
  39       ==
  40         betav*G + GR
  41       - fvm::SuSp(2.0/3.0*betav*rho*divU, k)
  42       - fvm::Sp(betav*rho*epsilon/k, k)
  43     );
  44
  45     bound(k, dimensionedScalar("0", k.dimensions(), 0.0));
  46
  47     //- Re-calculate turbulence viscosity
  48     mut = Cmu*rho*sqr(k)/epsilon;
  49
  50 #   include "wallViscosity.H"
  51
  52 }
  53
  54 muEff = mut + thermo->mu();