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

   1 {
   2     bool closedVolume = p_rgh.needReference();
   3     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
   4     bool compressible = (compressibility.value() > SMALL);
   5
   6     rho = thermo.rho();
   7
   8     volScalarField rAU(1.0/UEqn().A());
   9     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
  10
  11     U = rAU*UEqn().H();
  12
  13     surfaceScalarField phiU
  14     (
  15         fvc::interpolate(rho)
  16        *(
  17             (fvc::interpolate(U) & mesh.Sf())
  18           + fvc::ddtPhiCorr(rAU, rho, U, phi)
  19         )
  20     );
  21
  22     phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
  23
  24     {
  25         fvScalarMatrix p_rghDDtEqn
  26         (
  27             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
  28           + fvc::div(phi)
  29         );
  30
  31         // Thermodynamic density needs to be updated by psi*d(p) after the
  32         // pressure solution - done in 2 parts. Part 1:
  33         thermo.rho() -= psi*p_rgh;
  34
  35         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  36         {
  37             fvScalarMatrix p_rghEqn
  38             (
  39                 p_rghDDtEqn
  40               - fvm::laplacian(rhorAUf, p_rgh)
  41             );
  42
  43             p_rghEqn.solve
  44             (
  45                 mesh.solver
  46                 (
  47                     p_rgh.select
  48                     (
  49                         (
  50                            oCorr == nOuterCorr-1
  51                         && corr == nCorr-1
  52                         && nonOrth == nNonOrthCorr
  53                         )
  54                     )
  55                 )
  56             );
  57
  58             if (nonOrth == nNonOrthCorr)
  59             {
  60                 phi += p_rghEqn.flux();
  61             }
  62         }
  63
  64         // Second part of thermodynamic density update
  65         thermo.rho() += psi*p_rgh;
  66     }
  67
  68     // Correct velocity field
  69     U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
  70     U.correctBoundaryConditions();
  71
  72     p = p_rgh + rho*gh;
  73
  74     // Update pressure substantive derivative
  75     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
  76
  77     // Solve continuity
  78     #include "rhoEqn.H"
  79
  80     // Update continuity errors
  81     #include "compressibleContinuityErrors.H"
  82
  83     // For closed-volume cases adjust the pressure and density levels
  84     // to obey overall mass continuity
  85     if (closedVolume && compressible)
  86     {
  87         p += (initialMass - fvc::domainIntegrate(thermo.rho()))
  88             /compressibility;
  89         rho = thermo.rho();
  90         p_rgh = p - rho*gh;
  91     }
  92
  93     // Update thermal conductivity
  94     K = thermoFluid[i].Cp()*turb.alphaEff();
  95 }