applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H

   1 {
   2     volScalarField rUA = 1.0/UEqn.A();
   3
   4     surfaceScalarField psisf = fvc::interpolate(psis);
   5     surfaceScalarField rhof = fvc::interpolate(rho);
   6
   7     // Needs to be outside of loop since p is changing, but psi and rho are not.
   8     surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
   9
  10     for (int corr = 0; corr < nCorr; corr++)
  11     {
  12         U = rUA*UEqn.H();
  13
  14         // Calculate phi for boundary conditions
  15         phi = rhof*fvc::interpolate(U) & mesh.Sf();
  16
  17         surfaceScalarField phid2 = rhoReff/rhof*phi;
  18
  19         surfaceScalarField phid("phid", psisf/rhof*phi);
  20
  21         // Make fluxes relative within the MRF zone
  22         mrfZones.relativeFlux(rhoReff, phi);
  23         mrfZones.relativeFlux(psisf, phid);
  24         mrfZones.relativeFlux(rhoReff, phid2);
  25
  26         p.storePrevIter();
  27
  28         for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
  29         {
  30             fvScalarMatrix pEqn
  31             (
  32                 fvm::ddt(psis, p)
  33               + fvm::div(phid, p)
  34               + fvc::div(phid2)
  35               - fvm::laplacian(rho*rUA, p)
  36             );
  37
  38             // Retain the residual from the first pressure solution
  39             eqnResidual = pEqn.solve().initialResidual();
  40
  41             if (corr == 0 && nonOrth == 0)
  42             {
  43                 maxResidual = max(eqnResidual, maxResidual);
  44             }
  45
  46             // Calculate the flux
  47             if (nonOrth == nNonOrthCorr)
  48             {
  49                 phi = phid2 + pEqn.flux();
  50             }
  51         }
  52
  53 #       include "compressibleContinuityErrs.H"
  54
  55         // Relax the pressure
  56         p.relax();
  57
  58         U -= rUA*fvc::grad(p);
  59         U.correctBoundaryConditions();
  60     }
  61
  62     // Bound the pressure
  63     if (min(p) < pMin || max(p) > pMax)
  64     {
  65         p.max(pMin);
  66         p.min(pMax);
  67         p.correctBoundaryConditions();
  68     }
  69
  70     // Bound the velocity
  71     volScalarField magU = mag(U);
  72
  73     if (max(magU) > UMax)
  74     {
  75         volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU)
  76             + neg(magU - UMax);
  77         Ulimiter.max(scalar(0));
  78         Ulimiter.min(scalar(1));
  79
  80         U *= Ulimiter;
  81         U.correctBoundaryConditions();
  82     }
  83 }