applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

   1 {
   2     bool closedVolume = p.needReference();
   3
   4     rho = thermo.rho();
   5
   6     // Thermodynamic density needs to be updated by psi*d(p) after the
   7     // pressure solution - done in 2 parts. Part 1:
   8     thermo.rho() -= psi*p;
   9
  10     volScalarField rUA = 1.0/UEqn.A();
  11     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
  12
  13     U = rUA*UEqn.H();
  14
  15     surfaceScalarField phiU
  16     (
  17         fvc::interpolate(rho)
  18        *(
  19             (fvc::interpolate(U) & mesh.Sf())
  20           + fvc::ddtPhiCorr(rUA, rho, U, phi)
  21         )
  22     );
  23
  24     phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
  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(rhorUAf, p)
  33         );
  34
  35         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
  36         {
  37             pEqn.solve(mesh.solver(p.name() + "Final"));
  38         }
  39         else
  40         {
  41             pEqn.solve(mesh.solver(p.name()));
  42         }
  43
  44         if (nonOrth == nNonOrthCorr)
  45         {
  46             phi += pEqn.flux();
  47         }
  48     }
  49
  50     // Second part of thermodynamic density update
  51     thermo.rho() += psi*p;
  52
  53     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
  54     U.correctBoundaryConditions();
  55
  56     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
  57
  58     #include "rhoEqn.H"
  59     #include "compressibleContinuityErrs.H"
  60
  61     // For closed-volume cases adjust the pressure and density levels
  62     // to obey overall mass continuity
  63     if (closedVolume)
  64     {
  65         p +=
  66             (initialMass - fvc::domainIntegrate(psi*p))
  67             /fvc::domainIntegrate(psi);
  68         thermo.rho() = psi*p;
  69         rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
  70     }
  71 }