applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H

   1 rho = thermo.rho();
   2
   3 volScalarField rUA = 1.0/UEqn().A();
   4 U = rUA*UEqn().H();
   5
   6 if (nCorr <= 1)
   7 {
   8     UEqn.clear();
   9 }
  10
  11 if (transonic)
  12 {
  13     surfaceScalarField phid
  14     (
  15         "phid",
  16         fvc::interpolate(psi)
  17        *(
  18             (fvc::interpolate(U) & mesh.Sf())
  19           + fvc::ddtPhiCorr(rUA, rho, U, phi)
  20         )
  21     );
  22     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
  23
  24     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  25     {
  26         fvScalarMatrix pEqn
  27         (
  28             fvm::ddt(psi, p)
  29           + fvm::div(phid, p)
  30           - fvm::laplacian(rho*rUA, p)
  31         );
  32
  33         if
  34         (
  35             oCorr == nOuterCorr-1
  36             && corr == nCorr-1
  37             && nonOrth == nNonOrthCorr
  38         )
  39         {
  40             pEqn.solve(mesh.solver("pFinal"));
  41         }
  42         else
  43         {
  44             pEqn.solve();
  45         }
  46
  47         if (nonOrth == nNonOrthCorr)
  48         {
  49             phi == pEqn.flux();
  50         }
  51     }
  52 }
  53 else
  54 {
  55     phi =
  56         fvc::interpolate(rho)*
  57         (
  58             (fvc::interpolate(U) & mesh.Sf())
  59         //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
  60         );
  61     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
  62
  63     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  64     {
  65         // Pressure corrector
  66         fvScalarMatrix pEqn
  67         (
  68             fvm::ddt(psi, p)
  69           + fvc::div(phi)
  70           - fvm::laplacian(rho*rUA, p)
  71         );
  72
  73         if
  74         (
  75             oCorr == nOuterCorr-1
  76          && corr == nCorr-1
  77          && nonOrth == nNonOrthCorr
  78         )
  79         {
  80             pEqn.solve(mesh.solver("pFinal"));
  81         }
  82         else
  83         {
  84             pEqn.solve();
  85         }
  86
  87         if (nonOrth == nNonOrthCorr)
  88         {
  89             phi += pEqn.flux();
  90         }
  91     }
  92 }
  93
  94 #include "rhoEqn.H"
  95 #include "compressibleContinuityErrs.H"
  96
  97 //if (oCorr != nOuterCorr-1)
  98 {
  99     // Explicitly relax pressure for momentum corrector
 100     p.relax();
 101
 102     rho = thermo.rho();
 103     rho.relax();
 104     Info<< "rho max/min : " << max(rho).value()
 105         << " " << min(rho).value() << endl;
 106 }
 107
 108 U -= rUA*fvc::grad(p);
 109 U.correctBoundaryConditions();
 110
 111 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 112
 113 bound(p, pMin);
 114
 115 // For closed-volume cases adjust the pressure and density levels
 116 // to obey overall mass continuity
 117 /*
 118 if (closedVolume)
 119 {
 120     p += (initialMass - fvc::domainIntegrate(psi*p))
 121         /fvc::domainIntegrate(psi);
 122 }
 123 */