compressibleInterFoam family: Corrected transonic option

Resolves bug-report https://bugs.openfoam.org/view.php?id=2785

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -31,27 +31,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         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)
+            pos(alpha1)
-          + correction
+           *(
-            (
+                (
-                psi1*fvm::ddt(p_rgh)
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+            pos(alpha2)
-          + correction
+           *(
-            (
+               (
-                psi2*fvm::ddt(p_rgh)
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
-            );
+               )/rho2
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H

@@ -28,27 +28,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         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)
+            pos(alpha1)
-          + correction
+           *(
-            (
+                (
-                psi1*fvm::ddt(p_rgh)
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+            pos(alpha2)
-          + correction
+           *(
-            (
+               (
-                psi2*fvm::ddt(p_rgh)
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
-            );
+               )/rho2
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -28,27 +28,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         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)
+            pos(alpha1)
-          + correction
+           *(
-            (
+                (
-                psi1*fvm::ddt(p_rgh)
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+            pos(alpha2)
-          + correction
+           *(
-            (
+               (
-                psi2*fvm::ddt(p_rgh)
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
-            );
+               )/rho2
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -331,10 +331,6 @@ while (pimple.correct())
                     ).ptr()
                 );
 
-                deleteDemandDrivenData
-                (
-                    pEqnComps[phasei].faceFluxCorrectionPtr()
-                );
                 pEqnComps[phasei].relax();
             }
             else

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -255,7 +255,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
         else
@@ -297,7 +296,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
             pEqnComp2.ref().relax();
         }
         else

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -241,7 +241,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
         else
@@ -283,7 +282,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
             pEqnComp2.ref().relax();
         }
         else

applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H

@@ -262,7 +262,6 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
         pEqnComp1.ref().relax();
 
         pEqnComp2 =
@@ -278,7 +277,6 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
         pEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H

@@ -243,7 +243,6 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
         pEqnComp1.ref().relax();
 
         pEqnComp2 =
@@ -259,7 +258,6 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
         pEqnComp2.ref().relax();
     }
     else

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -1384,14 +1384,9 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::correction
 {
     tmp<Foam::fvMatrix<Type>> tAcorr = A - (A & A.psi());
 
-    if
+    // Delete the faceFluxCorrection from the correction matrix
-    (
+    // as it does not have a clear meaning or purpose
-        (A.hasUpper() || A.hasLower())
+    deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
-     && A.psi().mesh().fluxRequired(A.psi().name())
-    )
-    {
-        tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
-    }
 
     return tAcorr;
 }
@@ -1405,17 +1400,9 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::correction
 {
     tmp<Foam::fvMatrix<Type>> tAcorr = tA - (tA() & tA().psi());
 
-    // Note the matrix coefficients are still that of matrix A
+    // Delete the faceFluxCorrection from the correction matrix
-    const fvMatrix<Type>& A = tAcorr();
+    // as it does not have a clear meaning or purpose
-
+    deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
-    if
-    (
-        (A.hasUpper() || A.hasLower())
-     && A.psi().mesh().fluxRequired(A.psi().name())
-    )
-    {
-        tAcorr.ref().faceFluxCorrectionPtr() = (-A.flux()).ptr();
-    }
 
     return tAcorr;
 }