applications/solvers/compressible/rhoSimpleFoam/pEqn.H

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