applications/solvers/combustion/rhoReactingFoam/pEqn.H

   1 {
   2     rho = thermo.rho();
   3
   4     // Thermodynamic density needs to be updated by psi*d(p) after the
   5     // pressure solution - done in 2 parts. Part 1:
   6     thermo.rho() -= psi*p;
   7
   8     volScalarField rUA = 1.0/UEqn.A();
   9     U = rUA*UEqn.H();
  10
  11     if (pimple.transonic())
  12     {
  13         surfaceScalarField phiv =
  14             (fvc::interpolate(U) & mesh.Sf())
  15           + fvc::ddtPhiCorr(rUA, rho, U, phi);
  16
  17         phi = fvc::interpolate(rho)*phiv;
  18
  19         surfaceScalarField phid
  20         (
  21             "phid",
  22             fvc::interpolate(thermo.psi())*phiv
  23         );
  24
  25         while (pimple.correctNonOrthogonal())
  26         {
  27             fvScalarMatrix pEqn
  28             (
  29                 fvc::ddt(rho) + fvc::div(phi)
  30               + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
  31               - fvm::laplacian(rho*rUA, p)
  32             );
  33
  34             pEqn.solve
  35             (
  36                 mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
  37             );
  38
  39             if (pimple.finalNonOrthogonalIter())
  40             {
  41                 phi += pEqn.flux();
  42             }
  43         }
  44     }
  45     else
  46     {
  47         phi =
  48             fvc::interpolate(rho)
  49            *(
  50                 (fvc::interpolate(U) & mesh.Sf())
  51               + fvc::ddtPhiCorr(rUA, rho, U, phi)
  52             );
  53
  54         while (pimple.correctNonOrthogonal())
  55         {
  56             fvScalarMatrix pEqn
  57             (
  58                 fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  59               + fvc::div(phi)
  60               - fvm::laplacian(rho*rUA, p)
  61             );
  62
  63             pEqn.solve
  64             (
  65                 mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
  66             );
  67
  68             if (pimple.finalNonOrthogonalIter())
  69             {
  70                 phi += pEqn.flux();
  71             }
  72         }
  73     }
  74
  75     // Second part of thermodynamic density update
  76     thermo.rho() += psi*p;
  77
  78     #include "rhoEqn.H"
  79     #include "compressibleContinuityErrs.H"
  80
  81     U -= rUA*fvc::grad(p);
  82     U.correctBoundaryConditions();
  83
  84     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
  85 }