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 rAU(1.0/UEqn().A());
   7 U = rAU*UEqn().H();
   8 UEqn.clear();
   9
  10 bool closedVolume = false;
  11
  12 if (simple.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<=simple.nNonOrthCorr(); nonOrth++)
  21     {
  22         fvScalarMatrix pEqn
  23         (
  24             fvm::div(phid, p)
  25           - fvm::laplacian(rho*rAU, 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         pEqn.solve();
  34
  35         if (nonOrth == simple.nNonOrthCorr())
  36         {
  37             phi == pEqn.flux();
  38         }
  39     }
  40 }
  41 else
  42 {
  43     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
  44     closedVolume = adjustPhi(phi, U, p);
  45
  46     for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
  47     {
  48         fvScalarMatrix pEqn
  49         (
  50             fvm::laplacian(rho*rAU, p) == fvc::div(phi)
  51         );
  52
  53         pEqn.setReference(pRefCell, pRefValue);
  54
  55         pEqn.solve();
  56
  57         if (nonOrth == simple.nNonOrthCorr())
  58         {
  59             phi -= pEqn.flux();
  60         }
  61     }
  62 }
  63
  64
  65 #include "incompressible/continuityErrs.H"
  66
  67 // Explicitly relax pressure for momentum corrector
  68 p.relax();
  69
  70 U -= rAU*fvc::grad(p);
  71 U.correctBoundaryConditions();
  72
  73 // For closed-volume cases adjust the pressure and density levels
  74 // to obey overall mass continuity
  75 if (closedVolume)
  76 {
  77     p += (initialMass - fvc::domainIntegrate(psi*p))
  78         /fvc::domainIntegrate(psi);
  79 }
  80
  81 rho = thermo.rho();
  82 rho = max(rho, rhoMin);
  83 rho = min(rho, rhoMax);
  84 rho.relax();
  85 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;