applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/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 (pZones.size() > 0)
  12     {
  13         // ddtPhiCorr not well defined for cases with porosity
  14         phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
  15     }
  16     else
  17     {
  18         phi =
  19             fvc::interpolate(rho)
  20            *(
  21                 (fvc::interpolate(U) & mesh.Sf())
  22               + fvc::ddtPhiCorr(rAU, rho, U, phi)
  23             );
  24     }
  25
  26     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  27     {
  28         fvScalarMatrix pEqn
  29         (
  30             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  31           + fvc::div(phi)
  32           - fvm::laplacian(rho*rAU, p)
  33         ==
  34             parcels.Srho()
  35           + pointMassSources.Su()
  36         );
  37
  38         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
  39         {
  40             pEqn.solve(mesh.solver("pFinal"));
  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     #include "rhoEqn.H"
  57     #include "compressibleContinuityErrs.H"
  58
  59     U -= rAU*fvc::grad(p);
  60     U.correctBoundaryConditions();
  61 }