applications/solvers/multiphase/settlingFoam/kEpsilon.H

   1 if(turbulence)
   2 {
   3     if (mesh.changing())
   4     {
   5         y.correct();
   6     }
   7
   8     dimensionedScalar k0("k0", k.dimensions(), SMALL);
   9     dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), SMALL);
  10
  11     volScalarField divU = fvc::div(phi/fvc::interpolate(rho));
  12
  13     tmp<volTensorField> tgradU = fvc::grad(U);
  14     volScalarField G = 2*mut*(tgradU() && dev(symm(tgradU())));
  15     tgradU.clear();
  16
  17     volScalarField Gcoef =
  18         alphak*Cmu*k*(g & fvc::grad(rho))/(epsilon + epsilon0);
  19
  20 #   include "wallFunctions.H"
  21
  22     // Dissipation equation
  23     fvScalarMatrix epsEqn
  24     (
  25         fvm::ddt(rho, epsilon)
  26       + fvm::div(phi, epsilon)
  27       - fvm::laplacian
  28         (
  29             alphaEps*mut + mul, epsilon,
  30             "laplacian(DepsilonEff,epsilon)"
  31         )
  32      ==
  33         C1*G*epsilon/k
  34       - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
  35       - fvm::Sp(C2*rho*epsilon/k, epsilon)
  36     );
  37
  38 #   include "wallDissipation.H"
  39
  40     epsEqn.relax();
  41     epsEqn.solve();
  42
  43     bound(epsilon, epsilon0);
  44
  45
  46     // Turbulent kinetic energy equation
  47
  48     solve
  49     (
  50         fvm::ddt(rho, k)
  51       + fvm::div(phi, k)
  52       - fvm::laplacian(alphak*mut + mul, k, "laplacian(DkEff,k)")
  53      ==
  54         G
  55       - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
  56       - fvm::Sp(rho*epsilon/k, k)
  57     );
  58
  59     bound(k, k0);
  60
  61
  62     //- Re-calculate viscosity
  63     mut = rho*Cmu*sqr(k)/(epsilon + epsilon0);
  64
  65 #   include "wallViscosity.H"
  66 }
  67
  68 mu = mut + mul;