rhoSimplecFoam: further developments

not entirely satisfactory

applications/solvers/compressible/rhoSimplecFoam/UEqn.H

-    // Solve the Momentum equation
-
-    tmp<fvVectorMatrix> UEqn
-    (
-        fvm::div(phi, U)
-      + turbulence->divDevRhoReff(U)
-    );
-
-    UEqn().relax();
-
-    solve(UEqn() == -fvc::grad(p));

applications/solvers/compressible/rhoSimplecFoam/hEqn.H

-{
-    fvScalarMatrix hEqn
-    (
-        fvm::div(phi, h)
-      - fvm::Sp(fvc::div(phi), h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
-      - p*fvc::div(phi/fvc::interpolate(rho))
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    rho = thermo.rho();
-
-    if (!transonic)
-    {
-        rho.relax();
-    }
-
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
-}

applications/solvers/compressible/rhoSimplecFoam/pEqn.H

+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
 volScalarField p0 = p;
 
 volScalarField AU = UEqn().A();
@@ -24,8 +29,8 @@ if (transonic)
           + phid*(fvc::interpolate(p) - fvc::interpolate(p, "UD"))
         );
 
-        refCast<mixedFvPatchScalarField>(p.boundaryField()[1]).refValue()
-            = p.boundaryField()[1];
+        //refCast<mixedFvPatchScalarField>(p.boundaryField()[1]).refValue()
+        //    = p.boundaryField()[1];
 
         fvScalarMatrix pEqn
         (
@@ -35,6 +40,7 @@ if (transonic)
           + fvc::div(phid)*p
           - fvm::laplacian(rho/AtU, p)
         );
+        //pEqn.relax();
 
         pEqn.setReference(pRefCell, pRefValue);
 

applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C

@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
         #include "readSIMPLEControls.H"
 
         p.storePrevIter();
+        rho.storePrevIter();
 
         if (!transonic)
         {