solvers::isothermalFluid: Improved the consistency, accuracy and stability of the pressure-work term

for moving mesh cases, in both the internal energy and enthalpy equations

applications/solvers/combustion/XiFoam/EaEqn.H

@@ -7,7 +7,7 @@
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
             hea.name() == "ea"
-          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
           : -dpdt
         )
       + thermophysicalTransport.divq(hea)

applications/solvers/combustion/XiFoam/EauEqn.H

@@ -8,7 +8,7 @@ if (ign.ignited())
       + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
       + (
             heau.name() == "eau"
-          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
            *rho/thermo.rhou()
           : -dpdt*rho/thermo.rhou()
         )

applications/solvers/combustion/XiFoam/PDRFoam/EaEqn.H

@@ -7,7 +7,7 @@
       + betav*fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
             hea.name() == "ea"
-          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
           : -betav*dpdt
         )
       - fvm::laplacian(Db, hea)

applications/solvers/combustion/XiFoam/PDRFoam/EauEqn.H

@@ -8,7 +8,7 @@ if (ign.ignited())
       + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
       + (
             heau.name() == "eau"
-          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
            *rho/thermo.rhou()
           : -betav*dpdt*rho/thermo.rhou()
         )

applications/solvers/modules/fluid/thermophysicalPredictor.C

@@ -37,7 +37,7 @@ void Foam::solvers::fluid::thermophysicalPredictor()
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
             he.name() == "e"
-          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
+          ? addMotionWork(fvc::div(phi, p/rho))
           : -dpdt
         )
       + thermophysicalTransport->divq(he)

applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C

@@ -265,7 +265,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
 
         if (mesh.moving())
         {
-            dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
+            dpdt -= motionWork();
         }
     }
 }

applications/solvers/modules/isothermalFluid/correctPressure.C

@@ -241,7 +241,7 @@ void Foam::solvers::isothermalFluid::correctPressure()
 
         if (mesh.moving())
         {
-            dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
+            dpdt -= motionWork();
         }
     }
 }

applications/solvers/modules/isothermalFluid/isothermalFluid.C

@@ -54,6 +54,37 @@ void Foam::solvers::isothermalFluid::continuityErrors()
 }
 
 
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::solvers::isothermalFluid::motionWork() const
+{
+    return fvc::div
+    (
+        fvc::interpolate(rho)*fvc::meshPhi(rho, U),
+        p/rho,
+        "div(phi,(p|rho))"
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::solvers::isothermalFluid::addMotionWork
+(
+    const tmp<volScalarField>& work
+) const
+{
+    if (mesh.moving())
+    {
+        return work + motionWork();
+    }
+    else
+    {
+        return work;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solvers::isothermalFluid::isothermalFluid

applications/solvers/modules/isothermalFluid/isothermalFluid.H

@@ -193,6 +193,17 @@ private:
         void correctBuoyantPressure();
 
 
+protected:
+
+        //- Return the mesh-motion work
+        //  Used to correct dpdt of the enthalpy equation
+        //  or added to the div(Up) term of the internal energy equation
+        tmp<volScalarField> motionWork() const;
+
+        //- Adds the mesh-motion work to the work term provided
+        tmp<volScalarField> addMotionWork(const tmp<volScalarField>&) const;
+
+
 public:
 
     //- Runtime type information

applications/solvers/modules/multicomponentFluid/thermophysicalPredictor.C

@@ -51,9 +51,9 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
-                + mvConvection->fvmDiv(phi, Yi)
-                + thermophysicalTransport->divj(Yi)
-                ==
+              + mvConvection->fvmDiv(phi, Yi)
+              + thermophysicalTransport->divj(Yi)
+             ==
                 reaction->R(Yi)
               + fvModels().source(rho, Yi)
             );
@@ -79,7 +79,7 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
             he.name() == "e"
-          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+          ? addMotionWork(mvConvection->fvcDiv(phi, p/rho))
           : -dpdt
         )
       + thermophysicalTransport->divq(he)

tutorials/modules/multicomponentFluid/lockExchange/system/fvSchemes

@@ -35,7 +35,8 @@ divSchemes
         sludge          vanLeer;
         ".*"            upwind;
     };
-    div(phi,K)      Gauss limitedLinear 1;
+    div(phi,K)      Gauss upwind;
+    div(phi,(p|rho))    Gauss upwind;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }