applications/solvers/multiphase/twoPhaseEulerFoam/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, epsilon,
  22             "laplacian(DepsilonEff,epsilon)"
  23         )
  24       ==
  25          C1*beta*G*epsilon/k
  26        - fvm::Sp(C2*beta*epsilon/k, epsilon)
  27     );
  28
  29     #include "wallDissipation.H"
  30
  31     epsEqn.relax();
  32     epsEqn.solve();
  33
  34     epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
  35
  36
  37     // Turbulent kinetic energy equation
  38     fvScalarMatrix kEqn
  39     (
  40         fvm::ddt(beta, k)
  41       + fvm::div(phib, k)
  42       - fvm::laplacian
  43         (
  44             alphak*nuEffb, k,
  45             "laplacian(DkEff,k)"
  46         )
  47       ==
  48         beta*G
  49       - fvm::Sp(beta*epsilon/k, k)
  50     );
  51     kEqn.relax();
  52     kEqn.solve();
  53
  54     k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
  55
  56     //- Re-calculate turbulence viscosity
  57     nutb = Cmu*sqr(k)/epsilon;
  58
  59     #include "wallViscosity.H"
  60 }
  61
  62 nuEffb = nutb + nub;