applications/solvers/multiphase/compressibleInterFoam/pEqn.H

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