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

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