applications/solvers/compressible/rhoSimpleFoam/pEqn.H

   1 rho = thermo.rho();
   2 rho = max(rho, rhoMin);
   3 rho = min(rho, rhoMax);
   4 rho.relax();
   5
   6 volScalarField rUA = 1.0/UEqn().A();
   7 U = rUA*UEqn().H();
   8 UEqn.clear();
   9
  10 bool closedVolume = false;
  11
  12 if (transonic)
  13 {
  14     surfaceScalarField phid
  15     (
  16         "phid",
  17         fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
  18     );
  19
  20     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  21     {
  22         fvScalarMatrix pEqn
  23         (
  24             fvm::div(phid, p)
  25           - fvm::laplacian(rho*rUA, p)
  26         );
  27
  28         // Relax the pressure equation to ensure diagonal-dominance
  29         pEqn.relax(mesh.relaxationFactor("pEqn"));
  30
  31         pEqn.setReference(pRefCell, pRefValue);
  32
  33         // retain the residual from the first iteration
  34         if (nonOrth == 0)
  35         {
  36             eqnResidual = pEqn.solve().initialResidual();
  37             maxResidual = max(eqnResidual, maxResidual);
  38         }
  39         else
  40         {
  41             pEqn.solve();
  42         }
  43
  44         if (nonOrth == nNonOrthCorr)
  45         {
  46             phi == pEqn.flux();
  47         }
  48     }
  49 }
  50 else
  51 {
  52     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
  53     closedVolume = adjustPhi(phi, U, p);
  54
  55     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  56     {
  57         fvScalarMatrix pEqn
  58         (
  59             fvm::laplacian(rho*rUA, p) == fvc::div(phi)
  60         );
  61
  62         pEqn.setReference(pRefCell, pRefValue);
  63
  64         // Retain the residual from the first iteration
  65         if (nonOrth == 0)
  66         {
  67             eqnResidual = pEqn.solve().initialResidual();
  68             maxResidual = max(eqnResidual, maxResidual);
  69         }
  70         else
  71         {
  72             pEqn.solve();
  73         }
  74
  75         if (nonOrth == nNonOrthCorr)
  76         {
  77             phi -= pEqn.flux();
  78         }
  79     }
  80 }
  81
  82
  83 #include "incompressible/continuityErrs.H"
  84
  85 // Explicitly relax pressure for momentum corrector
  86 p.relax();
  87
  88 U -= rUA*fvc::grad(p);
  89 U.correctBoundaryConditions();
  90
  91 // For closed-volume cases adjust the pressure and density levels
  92 // to obey overall mass continuity
  93 if (closedVolume)
  94 {
  95     p += (initialMass - fvc::domainIntegrate(psi*p))
  96         /fvc::domainIntegrate(psi);
  97 }
  98
  99 rho = thermo.rho();
 100 rho = max(rho, rhoMin);
 101 rho = min(rho, rhoMax);
 102 rho.relax();
 103 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;