applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

   1 {
   2     rho = thermo.rho();
   3
   4     volScalarField rAU = 1.0/UEqn.A();
   5     U = rAU*UEqn.H();
   6
   7     if (pZones.size() > 0)
   8     {
   9         // ddtPhiCorr not well defined for cases with porosity
  10         phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
  11     }
  12     else
  13     {
  14         phi =
  15             fvc::interpolate(rho)
  16            *(
  17                 (fvc::interpolate(U) & mesh.Sf())
  18               + fvc::ddtPhiCorr(rAU, rho, U, phi)
  19             );
  20     }
  21
  22     {
  23         fvScalarMatrix pDDtEqn
  24         (
  25             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
  26           + fvc::div(phi)
  27         );
  28
  29         // Thermodynamic density needs to be updated by psi*d(p) after the
  30         // pressure solution - done in 2 parts. Part 1:
  31         thermo.rho() -= psi*p;
  32
  33         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  34         {
  35             fvScalarMatrix pEqn
  36             (
  37                 pDDtEqn - fvm::laplacian(rho*rAU, p)
  38              ==
  39                 parcels.Srho()
  40               + massSource.SuTot()
  41             );
  42
  43             if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
  44             {
  45                 pEqn.solve(mesh.solver("pFinal"));
  46             }
  47             else
  48             {
  49                 pEqn.solve();
  50             }
  51
  52             if (nonOrth == nNonOrthCorr)
  53             {
  54                 phi += pEqn.flux();
  55             }
  56         }
  57
  58         // Second part of thermodynamic density update
  59         thermo.rho() += psi*p;
  60     }
  61
  62     #include "rhoEqn.H"
  63     #include "compressibleContinuityErrs.H"
  64
  65     U -= rAU*fvc::grad(p);
  66     U.correctBoundaryConditions();
  67 }