applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

   1 {
   2     rho = thermo.rho();
   3     rho.relax();
   4
   5     volScalarField rAU(1.0/UEqn().A());
   6     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
   7
   8     U = rAU*UEqn().H();
   9     UEqn.clear();
  10
  11     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
  12     bool closedVolume = adjustPhi(phi, U, p_rgh);
  13
  14     surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
  15     phi -= buoyancyPhi;
  16
  17     for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
  18     {
  19         fvScalarMatrix p_rghEqn
  20         (
  21             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
  22         );
  23
  24         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
  25         p_rghEqn.solve();
  26
  27         if (nonOrth == simple.nNonOrthCorr())
  28         {
  29             // Calculate the conservative fluxes
  30             phi -= p_rghEqn.flux();
  31
  32             // Explicitly relax pressure for momentum corrector
  33             p_rgh.relax();
  34
  35             // Correct the momentum source with the pressure gradient flux
  36             // calculated from the relaxed pressure
  37             U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
  38             U.correctBoundaryConditions();
  39         }
  40     }
  41
  42     #include "continuityErrs.H"
  43
  44     p = p_rgh + rho*gh;
  45
  46     // For closed-volume cases adjust the pressure level
  47     // to obey overall mass continuity
  48     if (closedVolume)
  49     {
  50         p += (initialMass - fvc::domainIntegrate(psi*p))
  51             /fvc::domainIntegrate(psi);
  52         p_rgh = p - rho*gh;
  53     }
  54
  55     rho = thermo.rho();
  56     rho.relax();
  57     Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
  58         << endl;
  59 }