applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H

   1 {
   2     volScalarField rUA = 1.0/UEqn.A();
   3     surfaceScalarField rUAf = fvc::interpolate(rUA);
   4
   5     tmp<fvScalarMatrix> pdEqnComp;
   6
   7     if (transonic)
   8     {
   9         pdEqnComp =
  10             (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
  11     }
  12     else
  13     {
  14         pdEqnComp =
  15             (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
  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)
  30           - ghf*fvc::snGrad(rho)
  31         )*rUAf*mesh.magSf();
  32
  33     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  34     {
  35         fvScalarMatrix pdEqnIncomp
  36         (
  37             fvc::div(phi)
  38           - fvm::laplacian(rUAf, pd)
  39         );
  40
  41         solve
  42         (
  43             (
  44                 max(alpha1, scalar(0))*(psi1/rho1)
  45               + max(alpha2, scalar(0))*(psi2/rho2)
  46             )
  47            *pdEqnComp()
  48           + pdEqnIncomp
  49         );
  50
  51         if (nonOrth == nNonOrthCorr)
  52         {
  53             dgdt =
  54                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
  55                *(pdEqnComp & pd);
  56             phi += pdEqnIncomp.flux();
  57         }
  58     }
  59
  60     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
  61     U.correctBoundaryConditions();
  62
  63     p = max
  64         (
  65             (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
  66             pMin
  67         );
  68
  69     rho1 = rho10 + psi1*p;
  70     rho2 = rho20 + psi2*p;
  71
  72     Info<< "max(U) " << max(mag(U)).value() << endl;
  73     Info<< "min(pd) " << min(pd).value() << endl;
  74 }