Rationalized pEqn in sprayFoam solvers

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -62,7 +62,6 @@ else
 
     while (pimple.correctNonOrthogonal())
     {
-        // Pressure corrector
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)

applications/solvers/lagrangian/sprayFoam/UEqn.H

@@ -1,4 +1,6 @@
-    fvVectorMatrix UEqn
+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
 (
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
@@ -9,13 +11,13 @@
   + fvOptions(rho, U)
 );
 
-    UEqn.relax();
+UEqn().relax();
 
-    fvOptions.constrain(UEqn);
+fvOptions.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-        solve(UEqn == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -1,10 +1,18 @@
 rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
 
-volScalarField rAU(1.0/UEqn.A());
+volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
 
 if (pimple.transonic())
 {
@@ -13,9 +21,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-        )/fvc::interpolate(rho)
+        )
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -77,6 +85,17 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H

@@ -1,98 +0,0 @@
-rho = thermo.rho();
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
-
-if (pimple.transonic())
-{
-    surfaceScalarField phid
-    (
-        "phid",
-        fvc::interpolate(psi)
-       *(
-            (
-                (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
-            )/fvc::interpolate(rho)
-        )
-    );
-
-    fvc::makeRelative(phid, psi, U);
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            parcels.Srho()
-          + fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi == pEqn.flux();
-        }
-    }
-}
-else
-{
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
-        )
-    );
-
-    fvc::makeRelative(phiHbyA, rho, U);
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            parcels.Srho()
-          + fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
-    }
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-U = HbyA - rAU*fvc::grad(p);
-U.correctBoundaryConditions();
-fvOptions.correct(U);
-K = 0.5*magSqr(U);
-
-{
-    rhoUf = fvc::interpolate(rho*U);
-    surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
-}
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);
-}