applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

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