applications/solvers/multiphase/compressibleInterFoam/pEqn.H

   1 {
   2     volScalarField rUA = 1.0/UEqn.A();
   3     surfaceScalarField rUAf = fvc::interpolate(rUA);
   4
   5     tmp<fvScalarMatrix> pEqnComp;
   6
   7     if (transonic)
   8     {
   9         pEqnComp =
  10             (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
  11     }
  12     else
  13     {
  14         pEqnComp =
  15             (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
  16     }
  17
  18
  19     U = rUA*UEqn.H();
  20
  21     surfaceScalarField phiU
  22     (
  23         "phiU",
  24         (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
  25     );
  26
  27     phi = phiU +
  28         (
  29             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
  30           + fvc::interpolate(rho)*(g & mesh.Sf())
  31         )*rUAf;
  32
  33     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  34     {
  35         fvScalarMatrix pEqnIncomp
  36         (
  37             fvc::div(phi)
  38           - fvm::laplacian(rUAf, p)
  39         );
  40
  41         solve
  42         (
  43             (
  44                 max(alpha1, scalar(0))*(psi1/rho1)
  45               + max(alpha2, scalar(0))*(psi2/rho2)
  46             )
  47            *pEqnComp()
  48           + pEqnIncomp
  49         );
  50
  51         if (nonOrth == nNonOrthCorr)
  52         {
  53             dgdt =
  54                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
  55                *(pEqnComp & p);
  56             phi += pEqnIncomp.flux();
  57         }
  58     }
  59
  60     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
  61     U.correctBoundaryConditions();
  62
  63     p.max(pMin);
  64
  65     rho1 = rho10 + psi1*p;
  66     rho2 = rho20 + psi2*p;
  67
  68     Info<< "max(U) " << max(mag(U)).value() << endl;
  69     Info<< "min(p) " << min(p).value() << endl;
  70 }