compressible solvers: Name pressure diffusivity Dp

applications/solvers/compressible/rhoPimpleFoam/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
+wmake rhoPimplecFoam
 wmake rhoPorousMRFPimpleFoam
 wmake rhoPorousMRFLTSPimpleFoam
 

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -4,13 +4,12 @@ tmp<fvVectorMatrix> UEqn
 (
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
+  - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U)
   + turbulence->divDevRhoReff(U)
 );
 
 UEqn().relax();
 
-volScalarField rAU(1.0/UEqn().A());
-
 if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));

applications/solvers/compressible/rhoPimpleFoam/hEqn.H

@@ -3,10 +3,14 @@
     (
         fvm::ddt(rho, h)
       + fvm::div(phi, h)
+      - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
         dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
+      - (
+            fvc::ddt(rho, K) + fvc::div(phi, K)
+          - (fvc::ddt(rho) + fvc::div(phi))*K
+        )
     );
 
     hEqn.relax();

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -3,6 +3,7 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
+volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
@@ -23,13 +24,15 @@ if (pimple.transonic())
         )
     );
 
+    volScalarField Dp("Dp", rho*rAU);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -52,6 +55,8 @@ else
         )
     );
 
+    volScalarField Dp("Dp", rho*rAU);
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -59,7 +64,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -66,7 +66,10 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "rhoEqn.H"
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C

@@ -74,7 +74,10 @@ int main(int argc, char *argv[])
 
         #include "setrDeltaT.H"
 
-        #include "rhoEqn.H"
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H

@@ -13,8 +13,6 @@ UEqn().relax();
 mrfZones.addCoriolis(rho, UEqn());
 pZones.addResistance(UEqn());
 
-volScalarField rAU(1.0/UEqn().A());
-
 if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H

@@ -25,13 +25,15 @@ if (pimple.transonic())
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
+    volScalarField Dp("Dp", rho*rAU);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -56,6 +58,8 @@ else
 
     mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
+    volScalarField Dp("Dp", rho*rAU);
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -63,7 +67,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -69,7 +69,10 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "rhoEqn.H"
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())

applications/solvers/compressible/rhoSimpleFoam/hEqn.H

@@ -11,7 +11,6 @@
     );
 
     hEqn.relax();
-
     hEqn.solve();
 
     thermo.correct();

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H

@@ -1,61 +0,0 @@
-    Info<< "Reading thermophysical properties\n" << endl;
-
-    autoPtr<basicPsiThermo> pThermo
-    (
-        basicPsiThermo::New(mesh)
-    );
-    basicPsiThermo& thermo = pThermo();
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rho()
-    );
-
-    volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo.psi();
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "compressibleCreatePhi.H"
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
-
-    dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));
-    dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));
-
-    Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
-    (
-        compressible::RASModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo
-        )
-    );
-
-    dimensionedScalar initialMass = fvc::domainIntegrate(rho);

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -3,11 +3,11 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField AU(UEqn().A());
+volScalarField rAU(1.0/UEqn().A());
-volScalarField AtU(AU - UEqn().H1());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = UEqn().H()/AU;
+HbyA = rAU*UEqn().H();
 
 UEqn.clear();
 
@@ -15,30 +15,32 @@ bool closedVolume = false;
 
 if (simple.transonic())
 {
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    surfaceScalarField phic
+    (
+        "phic",
+        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
+    );
+
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField Dp("Dp", rho*rAtU);
+
     while (simple.correctNonOrthogonal())
     {
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
-        );
-
-        surfaceScalarField phic
-        (
-            "phic",
-            fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf()
-        );
-
-        HbyA -= fvc::grad(p)*(1.0/AU - 1.0/AtU);
-
         fvScalarMatrix pEqn
         (
             fvm::div(phid, p)
           + fvc::div(phic)
-          - fvm::laplacian(rho/AtU, p)
+          - fvm::laplacian(Dp, p)
         );
 
-        // Relax the pressure equation to ensure diagonal-dominance
+        // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -53,25 +55,25 @@ if (simple.transonic())
 }
 else
 {
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    closedVolume = adjustPhi(phiHbyA, U, p);
+
+    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField Dp("Dp", rho*rAtU);
+
     while (simple.correctNonOrthogonal())
     {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            fvc::interpolate(rho)*fvc::interpolate(HbyA) & mesh.Sf()
-        );
-
-        closedVolume = adjustPhi(phiHbyA, U, p);
-
-        phiHbyA +=
-            fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
-
-        HbyA -= fvc::grad(p)*(1.0/AU - 1.0/AtU);
-
         fvScalarMatrix pEqn
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(rho/AtU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -85,14 +87,14 @@ else
     }
 }
 
-// The incompressibe for of the continuity error check is appropriate for
+// The incompressibe form of the continuity error check is appropriate for
 // steady-state compressible also.
 #include "incompressible/continuityErrs.H"
 
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - fvc::grad(p)/AtU;
+U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 // For closed-volume cases adjust the pressure and density levels
@@ -103,6 +105,7 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
+// Recalculate density from the relaxed pressure
 rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -14,13 +14,16 @@ surfaceScalarField phid
     )
 );
 
+
+volScalarField Dp("Dp", rho*rAU);
+
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix pEqn
     (
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
-      - fvm::laplacian(rho*rAU, p)
+      - fvm::laplacian(Dp, p)
     );
 
     pEqn.solve();

applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H

@@ -14,13 +14,15 @@ surfaceScalarField phid
     )
 );
 
+volScalarField Dp("Dp", rho*rAU);
+
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix pEqn
     (
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
-      - fvm::laplacian(rho*rAU, p)
+      - fvm::laplacian(Dp, p)
     );
 
     pEqn.solve();

tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes

@@ -42,7 +42,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DBEff,B) Gauss linear corrected;

tutorials/compressible/rhoPimpleFoam/ras/angledDuct/0/U

@@ -48,7 +48,7 @@ boundaryField
     }
     outlet
     {
-        type            fluxCorrectedVelocity; //inletOutlet;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
         inletValue      uniform (0 0 0);
     }

tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes

@@ -52,7 +52,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian(DomegaEff,omega) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
 }
 

tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution

@@ -52,6 +52,7 @@ solvers
 PIMPLE
 {
     momentumPredictor yes;
+    transonic       no;
     nOuterCorrectors 50;
     nCorrectors     1;
     nNonOrthogonalCorrectors 0;
@@ -66,6 +67,8 @@ PIMPLE
             tolerance       0.0001;
         }
     }
+
+    turbOnFinalIterOnly off;
 }
 
 relaxationFactors

tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes

@@ -52,7 +52,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear orthogonal;
     laplacian(DREff,R) Gauss linear orthogonal;
     laplacian(DomegaEff,omega) Gauss linear orthogonal;
-    laplacian((rho*(1|A(U))),p) Gauss linear orthogonal;
+    laplacian(Dp,p) Gauss linear orthogonal;
     laplacian(alphaEff,h) Gauss linear orthogonal;
 }
 

tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes

@@ -52,7 +52,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian(DomegaEff,omega) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
 }
 

tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes

@@ -40,7 +40,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;

tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes

@@ -39,7 +39,6 @@ laplacianSchemes
 {
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian((rho|A(U)),p) Gauss linear corrected;
     laplacian((rho*rAU),p) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;

tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary

@@ -8,7 +8,7 @@
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;

tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes

@@ -40,7 +40,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(alphaEff,e) Gauss linear corrected;
 }

tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes

@@ -40,7 +40,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian((rho*(1|A(U))),p) Gauss linear orthogonal;
+    laplacian(Dp,p) Gauss linear orthogonal;
     laplacian(muEff,U) Gauss linear orthogonal;
     laplacian(alphaEff,e) Gauss linear orthogonal;
 }

tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes

@@ -46,7 +46,7 @@ laplacianSchemes
     laplacian(DkEff,k) Gauss linear limited 0.5;
     laplacian(DREff,R) Gauss linear limited 0.5;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.5;
+    laplacian(Dp,p) Gauss linear limited 0.5;
     laplacian(alphaEff,e) Gauss linear limited 0.5;
 }
 

tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes

@@ -46,7 +46,7 @@ laplacianSchemes
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
     laplacian(alphaEff,e) Gauss linear corrected;
 }
 

tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSchemes

@@ -37,7 +37,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(mu,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(Dp,p) Gauss linear corrected;
 }
 
 interpolationSchemes