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