compressibleInterFoam: Improve handling of compressibility and instate support for transonic flow

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -1,11 +1,11 @@
 {
-    solve
+    fvScalarMatrix TEqn
     (
         fvm::ddt(rho, T)
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
       + (
-            p*fvc::div(phi)
+            fvc::div(fvc::absolute(phi, U), p)
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(
@@ -14,5 +14,8 @@
         )
     );
 
+    TEqn.relax();
+    TEqn.solve();
+
     twoPhaseProperties.correct();
 }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -104,8 +104,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = alpha1*rho1 + alpha2*rho2;
-
         runTime.write();
 
         Info<< "ExecutionTime = "

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -38,9 +38,9 @@
 
     volScalarField& p = twoPhaseProperties.p();
     volScalarField& T = twoPhaseProperties.T();
-    const volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
+    volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
     const volScalarField& psi1 = twoPhaseProperties.thermo1().psi();
-    const volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
+    volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
     const volScalarField& psi2 = twoPhaseProperties.thermo2().psi();
 
     volScalarField rho
@@ -93,18 +93,5 @@
         compressible::turbulenceModel::New(rho, U, rhoPhi, twoPhaseProperties)
     );
 
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt
-    (
-        IOobject
-        (
-            "dpdt",
-            runTime.timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-    );
-
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -26,28 +26,44 @@
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
 
-    //if (transonic)
+    if (pimple.transonic())
-    //{
-    //}
-    //else
     {
         surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
         surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
 
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+          + correction
+            (
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
+        p_rghEqnComp1().relax();
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+          + correction
+            (
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
+        p_rghEqnComp2().relax();
+    }
+    else
+    {
         p_rghEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-          + fvc::div(phid1, p_rgh)
+          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
-          - fvc::Sp(fvc::div(phid1), p_rgh);
 
         p_rghEqnComp2 =
             fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-          + fvc::div(phid2, p_rgh)
+          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
-          - fvc::Sp(fvc::div(phid2), p_rgh);
     }
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // Cache p_rgh prior to solve for density update
-    // pressure solution - done in 2 parts. Part 1:
+    volScalarField p_rgh_0(p_rgh);
-    //twoPhaseProperties.rho() -= psi*p_rgh;
 
     while (pimple.correctNonOrthogonal())
     {
@@ -69,8 +85,8 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            // Second part of thermodynamic density update
+            //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            //twoPhaseProperties.rho() += psi*p_rgh;
+            //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
 
             dgdt =
             (
@@ -88,15 +104,14 @@
 
     p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
 
-    // twoPhaseProperties.correct();
+    // Update densities from change in p_rgh
+    rho1 += psi1*(p_rgh - p_rgh_0);
+    rho2 += psi2*(p_rgh - p_rgh_0);
 
-    Info<< "max(U) " << max(mag(U)).value() << endl;
+    rho = alpha1*rho1 + alpha2*rho2;
-    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
 
     K = 0.5*magSqr(U);
 
-    if (twoPhaseProperties.dpdt())
+    Info<< "max(U) " << max(mag(U)).value() << endl;
-    {
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
-        dpdt = fvc::ddt(p);
-    }
 }

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H

@@ -94,6 +94,16 @@ public:
             return thermo2_();
         }
 
+        rhoThermo& thermo1()
+        {
+            return thermo1_();
+        }
+
+        rhoThermo& thermo2()
+        {
+            return thermo2_();
+        }
+
         //- Update properties
         virtual void correct();
 

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes

@@ -34,6 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(rho*phi,K)  Gauss upwind;
+    div(phi,p)      Gauss upwind;
     div(phi,k)      Gauss vanLeer;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties

@@ -1,59 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      transportProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-phases (water air);
-
-water
-{
-    transportModel  Newtonian;
-    nu              1e-06;
-    rho             1000;
-    k               0; // 0.613;
-    Cv              4179;
-
-    equationOfState
-    {
-        type            perfectFluid;
-
-        rho0            1000;
-        R               3000;
-    }
-}
-
-air
-{
-    transportModel  Newtonian;
-    nu              1.589e-05;
-    rho             1;
-    k               0; // 2.63e-2;
-    Cv              721;
-
-    equationOfState
-    {
-        type            perfectFluid;
-
-        rho0            0;
-        R               287;
-    }
-}
-
-pMin            pMin [ 1 -1 -2 0 0 0 0 ] 10000;
-
-sigma           sigma [ 1 0 -2 0 0 0 0 ] 0.07;
-
-
-// ************************************************************************* //

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes

@@ -34,6 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(rho*phi,K)  Gauss upwind;
+    div(phi,p)      Gauss upwind;
     div(phi,k)      Gauss vanLeer;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }