applications/solvers/combustion/rhoReactingFoam/pEqn.H

   1 {
   2     rho = thermo.rho();
   3
   4     volScalarField rUA = 1.0/UEqn.A();
   5     U = rUA*UEqn.H();
   6
   7     if (transonic)
   8     {
   9         surfaceScalarField phiv =
  10             (fvc::interpolate(U) & mesh.Sf())
  11           + fvc::ddtPhiCorr(rUA, rho, U, phi);
  12
  13         phi = fvc::interpolate(rho)*phiv;
  14
  15         surfaceScalarField phid
  16         (
  17             "phid",
  18             fvc::interpolate(thermo.psi())*phiv
  19         );
  20
  21         fvScalarMatrix pDDtEqn
  22         (
  23             fvc::ddt(rho) + fvc::div(phi)
  24           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
  25         );
  26
  27         // Thermodynamic density needs to be updated by psi*d(p) after the
  28         // pressure solution - done in 2 parts. Part 1:
  29         thermo.rho() -= psi*p;
  30
  31         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  32         {
  33             fvScalarMatrix pEqn
  34             (
  35                 pDDtEqn - fvm::laplacian(rho*rUA, p)
  36             );
  37
  38             if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
  39             {
  40                 pEqn.solve(mesh.solver(p.name() + "Final"));
  41             }
  42             else
  43             {
  44                 pEqn.solve();
  45             }
  46
  47             if (nonOrth == nNonOrthCorr)
  48             {
  49                 phi += pEqn.flux();
  50             }
  51         }
  52
  53         // Second part of thermodynamic density update
  54         thermo.rho() += psi*p;
  55     }
  56     else
  57     {
  58         phi =
  59             fvc::interpolate(rho)
  60            *(
  61                 (fvc::interpolate(U) & mesh.Sf())
  62               + fvc::ddtPhiCorr(rUA, rho, U, phi)
  63             );
  64
  65         fvScalarMatrix pDDtEqn
  66         (
  67             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  68           + fvc::div(phi)
  69         );
  70
  71         // Thermodynamic density needs to be updated by psi*d(p) after the
  72         // pressure solution - done in 2 parts. Part 1:
  73         thermo.rho() -= psi*p;
  74
  75         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  76         {
  77             fvScalarMatrix pEqn
  78             (
  79                 pDDtEqn - fvm::laplacian(rho*rUA, p)
  80             );
  81
  82             if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
  83             {
  84                 pEqn.solve(mesh.solver(p.name() + "Final"));
  85             }
  86             else
  87             {
  88                 pEqn.solve();
  89             }
  90
  91             if (nonOrth == nNonOrthCorr)
  92             {
  93                 phi += pEqn.flux();
  94             }
  95         }
  96
  97         // Second part of thermodynamic density update
  98         thermo.rho() += psi*p;
  99     }
 100
 101     #include <finiteVolume/rhoEqn.H>
 102     #include <finiteVolume/compressibleContinuityErrs.H>
 103
 104     U -= rUA*fvc::grad(p);
 105     U.correctBoundaryConditions();
 106
 107     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 108 }