applications/solvers/multiphase/compressibleInterDyMFoam/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())
  25     );
  26
  27     phi = phiU +
  28         (
  29             fvc::interpolate(interface.sigmaK())*
  30             fvc::snGrad(alpha1)*mesh.magSf()
  31           + fvc::interpolate(rho)*(g & mesh.Sf())
  32         )*rUAf;
  33
  34     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  35     {
  36         fvScalarMatrix pEqnIncomp
  37         (
  38             fvc::div(phi)
  39           - fvm::laplacian(rUAf, p)
  40         );
  41
  42         if
  43         (
  44             oCorr == nOuterCorr-1
  45             && corr == nCorr-1
  46             && nonOrth == nNonOrthCorr
  47         )
  48         {
  49             solve
  50             (
  51                 (
  52                     max(alpha1, scalar(0))*(psi1/rho1)
  53                   + max(alpha2, scalar(0))*(psi2/rho2)
  54                 )
  55                *pEqnComp()
  56               + pEqnIncomp,
  57                 mesh.solutionDict().solver(p.name() + "Final")
  58             );
  59         }
  60         else
  61         {
  62             solve
  63             (
  64                 (
  65                     max(alpha1, scalar(0))*(psi1/rho1)
  66                   + max(alpha2, scalar(0))*(psi2/rho2)
  67                 )
  68                *pEqnComp()
  69               + pEqnIncomp
  70             );
  71         }
  72
  73         if (nonOrth == nNonOrthCorr)
  74         {
  75             dgdt =
  76                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
  77                *(pEqnComp & p);
  78             phi += pEqnIncomp.flux();
  79         }
  80     }
  81
  82     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
  83     U.correctBoundaryConditions();
  84
  85     p.max(pMin);
  86
  87     rho1 = rho10 + psi1*p;
  88     rho2 = rho20 + psi2*p;
  89
  90     Info<< "max(U) " << max(mag(U)).value() << endl;
  91     Info<< "min(p) " << min(p).value() << endl;
  92
  93     // Make the fluxes relative to the mesh motion
  94     fvc::makeRelative(phi, U);
  95 }