applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/kEpsilon.H

   1 if(turbulence)
   2 {
   3     if (mesh.changing())
   4     {
   5         y.correct();
   6     }
   7
   8     tmp<volTensorField> tgradUb = fvc::grad(Ub);
   9     volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
  10     tgradUb.clear();
  11
  12 #   include "wallFunctions.H"
  13
  14     // Dissipation equation
  15     fvScalarMatrix epsEqn
  16     (
  17         fvm::ddt(beta, epsilon)
  18       + fvm::div(phib, epsilon)
  19       - fvm::laplacian
  20         (
  21             alphaEps*nuEffb,
  22             epsilon,
  23             "laplacian((alphaEps*nuEffb),epsilon)"
  24         )
  25       ==
  26          C1*beta*G*epsilon/k
  27        - fvm::Sp(C2*beta*epsilon/k, epsilon)
  28     );
  29
  30 #   include "wallDissipation.H"
  31
  32     epsEqn.relax();
  33     epsEqn.solve();
  34
  35     epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
  36
  37
  38     // Turbulent kinetic energy equation
  39     fvScalarMatrix kEqn
  40     (
  41         fvm::ddt(beta, k)
  42       + fvm::div(phib, k)
  43       - fvm::laplacian(alphak*nuEffb, k, "laplacian((alphak*nuEffb),k)")
  44       ==
  45         beta*G
  46       - fvm::Sp(beta*epsilon/k, k)
  47     );
  48     kEqn.relax();
  49     kEqn.solve();
  50
  51     k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
  52
  53     //- Re-calculate turbulence viscosity
  54     nutb = Cmu*sqr(k)/epsilon;
  55
  56 #   include "wallViscosity.H"
  57
  58 }
  59
  60 nuEffa = sqr(Ct)*nutb + nua;
  61 nuEffb = nutb + nub;