applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

   1 {
   2     volScalarField rUA = 1.0/UEqn.A();
   3     surfaceScalarField rUAf = fvc::interpolate(rUA);
   4
   5     U = rUA*UEqn.H();
   6
   7     surfaceScalarField phiU
   8     (
   9         "phiU",
  10         (fvc::interpolate(U) & mesh.Sf())
  11       + fvc::ddtPhiCorr(rUA, rho, U, phi)
  12     );
  13
  14     adjustPhi(phiU, U, pd);
  15
  16     phi = phiU +
  17         (
  18             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
  19           - ghf*fvc::snGrad(rho)
  20         )*rUAf*mesh.magSf();
  21
  22
  23     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
  24     const volScalarField& vDotcP = vDotP[0]();
  25     const volScalarField& vDotvP = vDotP[1]();
  26
  27     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  28     {
  29         fvScalarMatrix pdEqn
  30         (
  31             fvc::div(phi) - fvm::laplacian(rUAf, pd)
  32           + (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
  33         );
  34
  35         pdEqn.setReference(pdRefCell, pdRefValue);
  36
  37         if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
  38         {
  39             pdEqn.solve(mesh.solutionDict().solver(pd.name() + "Final"));
  40         }
  41         else
  42         {
  43             pdEqn.solve(mesh.solutionDict().solver(pd.name()));
  44         }
  45
  46         if (nonOrth == nNonOrthCorr)
  47         {
  48             phi += pdEqn.flux();
  49         }
  50     }
  51
  52     p = pd + rho*gh;
  53
  54     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
  55     U.correctBoundaryConditions();
  56 }