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())*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         if
  42         (
  43             oCorr == nOuterCorr-1
  44             && corr == nCorr-1
  45             && nonOrth == nNonOrthCorr
  46         )
  47         {
  48             solve
  49             (
  50                 (
  51                     max(alpha1, scalar(0))*(psi1/rho1)
  52                   + max(alpha2, scalar(0))*(psi2/rho2)
  53                 )
  54                *pEqnComp()
  55               + pEqnIncomp,
  56                 mesh.solver(p.name() + "Final")
  57             );
  58         }
  59         else
  60         {
  61             solve
  62             (
  63                 (
  64                     max(alpha1, scalar(0))*(psi1/rho1)
  65                   + max(alpha2, scalar(0))*(psi2/rho2)
  66                 )
  67                *pEqnComp()
  68               + pEqnIncomp
  69             );
  70         }
  71
  72         if (nonOrth == nNonOrthCorr)
  73         {
  74             dgdt =
  75                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
  76                *(pEqnComp & p);
  77             phi += pEqnIncomp.flux();
  78         }
  79     }
  80
  81     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
  82     U.correctBoundaryConditions();
  83
  84     p.max(pMin);
  85
  86     rho1 = rho10 + psi1*p;
  87     rho2 = rho20 + psi2*p;
  88
  89     Info<< "max(U) " << max(mag(U)).value() << endl;
  90     Info<< "min(p) " << min(p).value() << endl;
  91
  92     // Make the fluxes relative to the mesh motion
  93     fvc::makeRelative(phi, U);
  94 }