applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

   1 fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
   2 fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
   3
   4 {
   5     {
   6         volTensorField gradUaT = fvc::grad(Ua)().T();
   7         volTensorField Rca
   8         (
   9             "Rca",
  10             ((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT
  11         );
  12
  13         if (kineticTheory.on())
  14         {
  15             Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
  16         }
  17
  18         surfaceScalarField phiRa =
  19             -fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
  20             /fvc::interpolate(alpha + scalar(0.001));
  21
  22         UaEqn =
  23         (
  24             (scalar(1) + Cvm*rhob*beta/rhoa)*
  25             (
  26                 fvm::ddt(Ua)
  27               + fvm::div(phia, Ua, "div(phia,Ua)")
  28               - fvm::Sp(fvc::div(phia), Ua)
  29             )
  30
  31           - fvm::laplacian(nuEffa, Ua)
  32           + fvc::div(Rca)
  33
  34           + fvm::div(phiRa, Ua, "div(phia,Ua)")
  35           - fvm::Sp(fvc::div(phiRa), Ua)
  36           + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
  37          ==
  38         //  g                          // Buoyancy term transfered to p-equation
  39           - fvm::Sp(beta/rhoa*K, Ua)
  40         //+ beta/rhoa*K*Ub             // Explicit drag transfered to p-equation
  41           - beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
  42         );
  43
  44         UaEqn.relax();
  45     }
  46
  47     {
  48         volTensorField gradUbT = fvc::grad(Ub)().T();
  49         volTensorField Rcb
  50         (
  51             "Rcb",
  52             ((2.0/3.0)*I)*(k + nuEffb*tr(gradUbT)) - nuEffb*gradUbT
  53         );
  54
  55         surfaceScalarField phiRb =
  56             -fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
  57             /fvc::interpolate(beta + scalar(0.001));
  58
  59         UbEqn =
  60         (
  61             (scalar(1) + Cvm*rhob*alpha/rhob)*
  62             (
  63                 fvm::ddt(Ub)
  64               + fvm::div(phib, Ub, "div(phib,Ub)")
  65               - fvm::Sp(fvc::div(phib), Ub)
  66             )
  67
  68           - fvm::laplacian(nuEffb, Ub)
  69           + fvc::div(Rcb)
  70
  71           + fvm::div(phiRb, Ub, "div(phib,Ub)")
  72           - fvm::Sp(fvc::div(phiRb), Ub)
  73
  74           + (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
  75          ==
  76         //  g                          // Buoyancy term transfered to p-equation
  77           - fvm::Sp(alpha/rhob*K, Ub)
  78         //+ alpha/rhob*K*Ua            // Explicit drag transfered to p-equation
  79           + alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
  80         );
  81
  82         UbEqn.relax();
  83     }
  84 }