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(T(fvc::grad(Ua)));
   7
   8         if (kineticTheory.on())
   9         {
  10             kineticTheory.solve(gradUaT);
  11             nuEffa = kineticTheory.mua()/rhoa;
  12         }
  13         else // If not using kinetic theory is using Ct model
  14         {
  15             nuEffa = sqr(Ct)*nutb + nua;
  16         }
  17
  18         volTensorField Rca
  19         (
  20             "Rca",
  21             ((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT
  22         );
  23
  24         if (kineticTheory.on())
  25         {
  26             Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
  27         }
  28
  29         surfaceScalarField phiRa
  30         (
  31             -fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
  32             /fvc::interpolate(alpha + scalar(0.001))
  33         );
  34
  35         UaEqn =
  36         (
  37             (scalar(1) + Cvm*rhob*beta/rhoa)*
  38             (
  39                 fvm::ddt(Ua)
  40               + fvm::div(phia, Ua, "div(phia,Ua)")
  41               - fvm::Sp(fvc::div(phia), Ua)
  42             )
  43
  44           - fvm::laplacian(nuEffa, Ua)
  45           + fvc::div(Rca)
  46
  47           + fvm::div(phiRa, Ua, "div(phia,Ua)")
  48           - fvm::Sp(fvc::div(phiRa), Ua)
  49           + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
  50          ==
  51         //  g                          // Buoyancy term transfered to p-equation
  52           - fvm::Sp(beta/rhoa*K, Ua)
  53         //+ beta/rhoa*K*Ub             // Explicit drag transfered to p-equation
  54           - beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
  55         );
  56
  57         UaEqn.relax();
  58     }
  59
  60     {
  61         volTensorField gradUbT(T(fvc::grad(Ub)));
  62         volTensorField Rcb
  63         (
  64             "Rcb",
  65             ((2.0/3.0)*I)*(k + nuEffb*tr(gradUbT)) - nuEffb*gradUbT
  66         );
  67
  68         surfaceScalarField phiRb
  69         (
  70             -fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
  71             /fvc::interpolate(beta + scalar(0.001))
  72         );
  73
  74         UbEqn =
  75         (
  76             (scalar(1) + Cvm*rhob*alpha/rhob)*
  77             (
  78                 fvm::ddt(Ub)
  79               + fvm::div(phib, Ub, "div(phib,Ub)")
  80               - fvm::Sp(fvc::div(phib), Ub)
  81             )
  82
  83           - fvm::laplacian(nuEffb, Ub)
  84           + fvc::div(Rcb)
  85
  86           + fvm::div(phiRb, Ub, "div(phib,Ub)")
  87           - fvm::Sp(fvc::div(phiRb), Ub)
  88
  89           + (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
  90          ==
  91         //  g                          // Buoyancy term transfered to p-equation
  92           - fvm::Sp(alpha/rhob*K, Ub)
  93         //+ alpha/rhob*K*Ua            // Explicit drag transfered to p-equation
  94           + alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
  95         );
  96
  97         UbEqn.relax();
  98     }
  99 }