applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H

   1 volScalarField rUA = 1.0/UEqn().A();
   2 U = rUA*UEqn().H();
   3 UEqn.clear();
   4 phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
   5 bool closedVolume = adjustPhi(phi, U, p);
   6 phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
   7
   8 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
   9 {
  10     fvScalarMatrix pdEqn
  11     (
  12         fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
  13     );
  14
  15     pdEqn.setReference(pdRefCell, pdRefValue);
  16     // retain the residual from the first iteration
  17     if (nonOrth == 0)
  18     {
  19         eqnResidual = pdEqn.solve().initialResidual();
  20         maxResidual = max(eqnResidual, maxResidual);
  21     }
  22     else
  23     {
  24         pdEqn.solve();
  25     }
  26
  27     if (nonOrth == nNonOrthCorr)
  28     {
  29         phi -= pdEqn.flux();
  30     }
  31 }
  32
  33 #include "continuityErrs.H"
  34
  35 // Explicitly relax pressure for momentum corrector
  36 pd.relax();
  37
  38 p = pd + rho*gh + pRef;
  39
  40 U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
  41 U.correctBoundaryConditions();
  42
  43 // For closed-volume cases adjust the pressure and density levels
  44 // to obey overall mass continuity
  45 if (closedVolume)
  46 {
  47     p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
  48         /fvc::domainIntegrate(thermo->psi());
  49     pd == p - (rho*gh + pRef);
  50 }
  51
  52 rho = thermo->rho();
  53 rho.relax();
  54 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;