Compressible solver robustness improvements

applications/solvers/compressible/sonicDyMFoam/UEqn.H

@@ -7,7 +7,10 @@
 
     UEqn.relax
     (
-        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()))
+        mesh.solutionDict().equationRelaxationFactor
+        (
+            U.select(pimple.finalIter())
+        )
     );
 
     solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -27,8 +27,6 @@
             fvm::ddt(psis, p)
           + fvm::div(phid, p)
             // Convective flux relaxation terms
-          + fvm::SuSp(-divPhid, p)
-          + divPhid*p
           + fvc::div(phid2)
           - fvm::laplacian(rho*rUA, p)
         );

applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H

@@ -8,7 +8,10 @@
 
     UEqn.relax
     (
-        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()))
+        mesh.solutionDict().equationRelaxationFactor
+        (
+            U.select(pimple.finalIter())
+        )
     );
 
     solve(UEqn == -fvc::grad(p));

applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H

@@ -12,4 +12,8 @@
     hEqn.solve();
 
     thermo.correct();
+
+    // Recalculate density
+    rho = thermo.rho();
+    rho.correctBoundaryConditions();
 }

applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H

@@ -1,6 +1,4 @@
 {
-    rho = thermo.rho();
-
     rUA = 1.0/UEqn.A();
     U = rUA*UEqn.H();
 
@@ -56,12 +54,10 @@
         p.relax();
     }
 
-#   include "rhoEqn.H"
 #   include "compressibleContinuityErrs.H"
 
     U -= rUA*fvc::grad(p);
     U.correctBoundaryConditions();
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-    dpdt = fvc::ddt(p);
 }

applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C

@@ -79,14 +79,15 @@ int main(int argc, char *argv[])
 
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
-        // Make flux absolute
-        phi += meshFlux;
+        // Make the fluxes absolute (using the ddt(rho, U) scheme)
+        phi += fvc::interpolate(rho)*fvc::meshPhi(rho, U);
 
         bool meshChanged = mesh.update();
 
 #       include "volContinuity.H"
 
-        mesh.setBoundaryVelocity(U);
+        // Make the fluxes relative (using the ddt(rho, U) scheme)
+        phi -= fvc::interpolate(rho)*fvc::meshPhi(rho, U);
 
         if (meshChanged)
         {
@@ -95,14 +96,7 @@ int main(int argc, char *argv[])
             rho.correctBoundaryConditions();
         }
 
-        meshFlux = fvc::interpolate(rho)*fvc::meshPhi(rho, U);
-
-        phi = fvc::interpolate(rho)
-            *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
-
-        DpDt = dpdt + fvc::div(phi/fvc::interpolate(rho), p)
-            - fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(U))*p;
-
+        if (meshChanged)
         {
 #           include "compressibleCourantNo.H"
         }
@@ -111,22 +105,20 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
 #           include "rhoEqn.H"
+#           include "hEqn.H"
 #           include "UEqn.H"
 
             // --- PISO loop
             while (pimple.correct())
             {
 #               include "pEqn.H"
-#               include "hEqn.H"
             }
+
+            turbulence->correct();
         }
 
-        turbulence->correct();
-
 #       include "logSummary.H"
 
-        rho = thermo.rho();
-
         runTime.write();
 
 #       include "infoDataOutput.H"