applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

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