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 rAU(1.0/UEqn.A());
   9     U = rAU*UEqn.H();
  10
  11     if (pimple.transonic())
  12     {
  13         surfaceScalarField phiv
  14         (
  15             (fvc::interpolate(U) & mesh.Sf())
  16           + fvc::ddtPhiCorr(rAU, rho, U, phi)
  17         );
  18
  19         phi = fvc::interpolate(rho)*phiv;
  20
  21         surfaceScalarField phid
  22         (
  23             "phid",
  24             fvc::interpolate(thermo.psi())*phiv
  25         );
  26
  27         fvScalarMatrix pDDtEqn
  28         (
  29             fvc::ddt(rho) + fvc::div(phi)
  30           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
  31         );
  32
  33         for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
  34         {
  35             fvScalarMatrix pEqn
  36             (
  37                 pDDtEqn
  38               - fvm::laplacian(rho*rAU, p)
  39             );
  40
  41             pEqn.solve
  42             (
  43                 mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  44             );
  45
  46             if (nonOrth == pimple.nNonOrthCorr())
  47             {
  48                 phi += pEqn.flux();
  49             }
  50         }
  51     }
  52     else
  53     {
  54         phi =
  55             fvc::interpolate(rho)
  56            *(
  57                 (fvc::interpolate(U) & mesh.Sf())
  58               + fvc::ddtPhiCorr(rAU, rho, U, phi)
  59             );
  60
  61         fvScalarMatrix pDDtEqn
  62         (
  63             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  64           + fvc::div(phi)
  65         );
  66
  67         for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
  68         {
  69             fvScalarMatrix pEqn
  70             (
  71                 pDDtEqn
  72               - fvm::laplacian(rho*rAU, p)
  73             );
  74
  75             pEqn.solve
  76             (
  77                 mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  78             );
  79
  80             if (nonOrth == pimple.nNonOrthCorr())
  81             {
  82                 phi += pEqn.flux();
  83             }
  84         }
  85     }
  86
  87     // Second part of thermodynamic density update
  88     thermo.rho() += psi*p;
  89
  90     #include "rhoEqn.H"
  91     #include "compressibleContinuityErrs.H"
  92
  93     U -= rAU*fvc::grad(p);
  94     U.correctBoundaryConditions();
  95
  96     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
  97 }