applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H

   1 rho = thermo.rho();
   2 rho = max(rho, rhoMin);
   3 rho = min(rho, rhoMax);
   4 rho.relax();
   5
   6 volScalarField rAU(1.0/UEqn().A());
   7 U = rAU*UEqn().H();
   8
   9 if (pimple.nCorr() <= 1)
  10 {
  11     UEqn.clear();
  12 }
  13
  14 if (pimple.transonic())
  15 {
  16     surfaceScalarField phid
  17     (
  18         "phid",
  19         fvc::interpolate(psi)
  20        *(
  21             (fvc::interpolate(U) & mesh.Sf())
  22           + fvc::ddtPhiCorr(rAU, rho, U, phi)
  23         )
  24     );
  25     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
  26
  27     for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
  28     {
  29         fvScalarMatrix pEqn
  30         (
  31             fvm::ddt(psi, p)
  32           + fvm::div(phid, p)
  33           - fvm::laplacian(rho*rAU, p)
  34         );
  35
  36         pEqn.solve
  37         (
  38             mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  39         );
  40
  41         if (nonOrth == pimple.nNonOrthCorr())
  42         {
  43             phi == pEqn.flux();
  44         }
  45     }
  46 }
  47 else
  48 {
  49     phi =
  50         fvc::interpolate(rho)*
  51         (
  52             (fvc::interpolate(U) & mesh.Sf())
  53           + fvc::ddtPhiCorr(rAU, rho, U, phi)
  54         );
  55     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
  56
  57     for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
  58     {
  59         // Pressure corrector
  60         fvScalarMatrix pEqn
  61         (
  62             fvm::ddt(psi, p)
  63           + fvc::div(phi)
  64           - fvm::laplacian(rho*rAU, p)
  65         );
  66
  67         pEqn.solve
  68         (
  69             mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  70         );
  71
  72         if (nonOrth == pimple.nNonOrthCorr())
  73         {
  74             phi += pEqn.flux();
  75         }
  76     }
  77 }
  78
  79 #include "rhoEqn.H"
  80 #include "compressibleContinuityErrs.H"
  81
  82 // Explicitly relax pressure for momentum corrector
  83 p.relax();
  84
  85 rho = thermo.rho();
  86 rho = max(rho, rhoMin);
  87 rho = min(rho, rhoMax);
  88 rho.relax();
  89 Info<< "rho max/min : " << max(rho).value()
  90     << " " << min(rho).value() << endl;
  91
  92 U -= rAU*fvc::grad(p);
  93 U.correctBoundaryConditions();
  94
  95 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);