applications/solvers/multiphase/bubbleFoam/UEqns.H

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