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 (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         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  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             if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
  35             {
  36                 pEqn.solve(mesh.solver(p.name() + "Final"));
  37             }
  38             else
  39             {
  40                 pEqn.solve();
  41             }
  42
  43             if (nonOrth == nNonOrthCorr)
  44             {
  45                 phi += pEqn.flux();
  46             }
  47         }
  48     }
  49     else
  50     {
  51         phi =
  52             fvc::interpolate(rho)
  53            *(
  54                 (fvc::interpolate(U) & mesh.Sf())
  55               + fvc::ddtPhiCorr(rUA, rho, U, phi)
  56             );
  57
  58         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  59         {
  60             fvScalarMatrix pEqn
  61             (
  62                 fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  63               + fvc::div(phi)
  64               - fvm::laplacian(rho*rUA, p)
  65             );
  66
  67             if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
  68             {
  69                 pEqn.solve(mesh.solver(p.name() + "Final"));
  70             }
  71             else
  72             {
  73                 pEqn.solve();
  74             }
  75
  76             if (nonOrth == nNonOrthCorr)
  77             {
  78                 phi += pEqn.flux();
  79             }
  80         }
  81     }
  82
  83     // Second part of thermodynamic density update
  84     thermo.rho() += psi*p;
  85
  86     #include "rhoEqn.H"
  87     #include "compressibleContinuityErrs.H"
  88
  89     U -= rUA*fvc::grad(p);
  90     U.correctBoundaryConditions();
  91
  92     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
  93 }