Compressible and reacting flow solvers: Changed the internal energy pressure work for consistency with enthalpy

The pressure work term for total internal energy is div(U p) which can be discretised is various ways, given a mass flux field phi it seems logical to implement it in the form div(phi/interpolate(rho), p) but this is not exactly consistent with the relationship between enthalpy and internal energy (h = e + p/rho) and the transport of enthalpy, it would be more consistent to implement it in the form div(phi, p/rho). A further improvement in consistency can be gained by using the same convection scheme for this work term and the convection term div(phi, e) and for reacting solvers this is easily achieved by using the multi-variate limiter mvConvection provided for energy and specie convection. This more consistent total internal energy work term has now been implemented in all the compressible and reacting flow solvers and provides more accurate solutions when running with internal energy, particularly for variable density mixing cases with small pressure variation. For non-reacting compressible solvers this improvement requires a change to the corresponding divScheme in fvSchemes: div(phiv,p) -> div(phi,(p|rho)) and all the tutorials have been updated accordingly.
2021-06-11 11:34:38 +01:00
parent 881bbfa591
commit 90831fbb55
32 changed files with 162 additions and 187 deletions
--- a/applications/solvers/combustion/PDRFoam/EaEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EaEqn.H
@ -7,12 +7,7 @@
      + betav*fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            hea.name() == "ea"
-          ? fvc::div
-            (
-                phi/fvc::interpolate(rho),
-                p,
-                "div(phiv,p)"
-            )
+          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
          : -betav*dpdt
        )
      - fvm::laplacian(Db, hea)
--- a/applications/solvers/combustion/PDRFoam/EauEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EauEqn.H
@ -8,12 +8,8 @@ if (ign.ignited())
      + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
      + (
            heau.name() == "eau"
-          ? fvc::div
-            (
-                phi/fvc::interpolate(rho),
-                p,
-                "div(phiv,p)"
-            )*rho/thermo.rhou()
+          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+           *rho/thermo.rhou()
          : -betav*dpdt*rho/thermo.rhou()
        )
      - fvm::laplacian(Db, heau)
--- a/applications/solvers/combustion/XiFoam/EaEqn.H
+++ b/applications/solvers/combustion/XiFoam/EaEqn.H
@ -7,12 +7,7 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            hea.name() == "ea"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      - fvm::laplacian(thermophysicalTransport->alphaEff(), hea)
--- a/applications/solvers/combustion/XiFoam/EauEqn.H
+++ b/applications/solvers/combustion/XiFoam/EauEqn.H
@ -8,12 +8,8 @@ if (ign.ignited())
      + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
      + (
            heau.name() == "eau"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )*rho/thermo.rhou()
+          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+           *rho/thermo.rhou()
          : -dpdt*rho/thermo.rhou()
        )
      - fvm::laplacian(thermophysicalTransport->alphaEff(), heau)
--- a/applications/solvers/combustion/reactingFoam/EEqn.H
+++ b/applications/solvers/combustion/reactingFoam/EEqn.H
@ -7,12 +7,7 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      + thermophysicalTransport->divq(he)
--- a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
@ -7,12 +7,7 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      + thermophysicalTransport->divq(he)
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H
@ -7,12 +7,7 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      + thermophysicalTransport->divq(he)
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
@ -12,12 +12,9 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? mvConvection.valid()
+            ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
+            : fvc::div(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      + thermophysicalTransport.divq(he)
--- a/applications/solvers/heatTransfer/thermoFoam/EEqn.H
+++ b/applications/solvers/heatTransfer/thermoFoam/EEqn.H
@ -7,12 +7,7 @@
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
+          ? fvc::div(fvc::absolute(phi, rho, U), p/rho)
          : -dpdt
        )
      + thermophysicalTransport->divq(he)
--- a/tutorials/combustion/XiEngineFoam/kivaTest/system/fvSchemes
+++ b/tutorials/combustion/XiEngineFoam/kivaTest/system/fvSchemes
@ -32,7 +32,6 @@ divSchemes
    div(phi,K)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(meshPhi,p)  Gauss upwind;
-    div(phiv,p)     Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,R)      Gauss upwind;
--- a/tutorials/combustion/coldEngineFoam/freePiston/system/fvSchemes
+++ b/tutorials/combustion/coldEngineFoam/freePiston/system/fvSchemes
@ -32,7 +32,7 @@ divSchemes
    div(phi,U)          Gauss linearUpwind grad(U);
    div(phi,e)          Gauss linearUpwind grad(e);

-    div(phiv,p)     Gauss linear;
+    div(phi,(p|rho))    Gauss linear;
    div(phi,K)          Gauss linear;
    div(meshPhi,p)      Gauss linear;

--- a/tutorials/combustion/reactingFoam/Lagrangian/filter/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/Lagrangian/filter/system/fvSchemes
@ -32,7 +32,6 @@ divSchemes
    div(phi,U)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(phi,K)      Gauss linear;
-    div(phiv,p)     Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(U)          Gauss linear;
--- a/tutorials/combustion/reactingFoam/Lagrangian/parcelInBox/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/Lagrangian/parcelInBox/system/fvSchemes
@ -32,7 +32,6 @@ divSchemes
    div(phi,U)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(phi,K)      Gauss upwind;
-    div(phiv,p)     Gauss upwind;
    div(phi,e)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
--- a/tutorials/combustion/reactingFoam/Lagrangian/verticalChannel/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/Lagrangian/verticalChannel/system/fvSchemes
@ -31,7 +31,6 @@ divSchemes
    div(phi,U)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(phi,K)      Gauss upwind;
-    div(phiv,p)     Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega) Gauss upwind;
--- a/tutorials/combustion/reactingFoam/Lagrangian/verticalChannelLTS/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/Lagrangian/verticalChannelLTS/system/fvSchemes
@ -31,7 +31,6 @@ divSchemes
    div(phi,U)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(phi,K)      Gauss upwind;
-    div(phiv,p)     Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega) Gauss upwind;
--- a/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSchemes
@ -28,14 +28,16 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss LUST grad(U);
    div(phi,e)          Gauss LUST grad(e);
    div(phi,K)          Gauss linear;
-    div(phiv,p)     Gauss linear;
+    div(phi,(p|rho))    Gauss linear;
    div(phi,k)          Gauss limitedLinear 1;
    div(phi,B)          Gauss limitedLinear 1;
    div(phi,muTilda)    Gauss limitedLinear 1;
    div(B)              Gauss linear;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/RAS/aerofoilNACA0012/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/aerofoilNACA0012/system/fvSchemes
@ -44,8 +44,8 @@ divSchemes
    div(phi,e)          $energy;
    div(phi,K)          $energy;
    div(phi,Ekp)        $energy;
+    div(phi,(p|rho))    $energy;

-    div(phiv,p)     Gauss upwind;
    div((phi|interpolate(rho)),p) Gauss upwind;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes
@ -28,18 +28,21 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          bounded Gauss upwind;
    div(phid,p)         bounded Gauss upwind;
-    div(phiv,p)     bounded Gauss linear;
+    div(phi,(p|rho))    bounded Gauss upwind;
    div(phi,K)          bounded Gauss linear;
    div(phi,e)          bounded Gauss upwind;
    div(phi,k)          bounded Gauss upwind;
    div(phi,epsilon)    bounded Gauss upwind;
    div(phi,R)          bounded Gauss upwind;
    div(phi,omega)      bounded Gauss upwind;
+
    div((rho*R))        Gauss linear;
    div(R)              Gauss linear;
    div(U)              Gauss linear;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/RAS/nacaAirfoil/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/nacaAirfoil/system/fvSchemes
@ -28,13 +28,15 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss limitedLinearV 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phid,p)         Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
    div(phi,k)          Gauss upwind;
    div(phi,epsilon)    Gauss upwind;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/RAS/prism/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/prism/system/fvSchemes
@ -28,13 +28,15 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss limitedLinearV 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phid,p)         Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
    div(phi,k)          Gauss upwind;
    div(phi,epsilon)    Gauss upwind;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
@ -36,7 +36,7 @@ divSchemes
    div(phi,epsilon)    Gauss linearUpwind limited;
    div(phi,k)          Gauss linearUpwind limited;
    div(phi,K)          Gauss linearUpwind limited;
-    div(phiv,p)         Gauss linearUpwind limited;
+    div(phi,(p|rho))    Gauss linearUpwind limited;

    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
 }
--- a/tutorials/compressible/rhoPimpleFoam/laminar/blockedChannel/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/blockedChannel/system/fvSchemes
@ -29,11 +29,13 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss linearUpwind limited;
    div(phi,e)          Gauss linearUpwind limited;
    div(phi,tracer)     Gauss linearUpwind limited;
    div(phi,K)          Gauss linear;
-    div(phiv,p)         Gauss linear;
+    div(phi,(p|rho))    Gauss linear;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/laminar/decompressionTank/decompressionTank/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/decompressionTank/decompressionTank/system/fvSchemes
@ -33,7 +33,7 @@ divSchemes
    div(phid,p)         Gauss limitedLinear 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
--- a/tutorials/compressible/rhoPimpleFoam/laminar/forwardStep/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/forwardStep/system/fvSchemes
@ -32,7 +32,8 @@ divSchemes
    div(phid,p)         Gauss limitedLinear 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSchemes
@ -28,10 +28,12 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss limitedLinearV 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/compressible/rhoPimpleFoam/laminar/shockTube/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/shockTube/system/fvSchemes
@ -28,11 +28,13 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          Gauss upwind;
    div(phid,p)         Gauss limitedLinear 1;
    div(phi,e)          Gauss limitedLinear 1;
    div(phi,K)          Gauss limitedLinear 1;
-    div(phiv,p)     Gauss limitedLinear 1;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/heatTransfer/buoyantPimpleFoam/BernardCells/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/BernardCells/system/fvSchemes
@ -36,7 +36,7 @@ divSchemes
    div(phi,k)          $turbulence;
    div(phi,epsilon)    $turbulence;

-    div(phiv,p)     Gauss linear;
+    div(phi,(p|rho))    Gauss limitedLinear 1;
    div(phi,K)          Gauss linear;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomBoussinesq/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomBoussinesq/system/fvSchemes
@ -35,7 +35,7 @@ divSchemes
    div(phi,k)          Gauss upwind;
    div(phi,epsilon)    Gauss upwind;

-    div(phiv,p)     Gauss linear;
+    div(phi,(p|rho))    Gauss linear;
    div(phi,K)          Gauss linear;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
--- a/tutorials/heatTransfer/chtMultiRegionFoam/coolingCylinder2D/system/fluid/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/coolingCylinder2D/system/fluid/fvSchemes
@ -40,8 +40,8 @@ divSchemes
    div(phi,omega)      $turbulence;
    div(phi,e)          $turbulence;
    div(phi,K)          $turbulence;
+    div(phi,(p|rho))    $turbulence;

-    div(phiv,p)     Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
 }

--- a/tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/air/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/air/fvSchemes
@ -29,13 +29,15 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          bounded Gauss upwind;
    div(phi,h)          bounded Gauss upwind;
    div(phi,e)          bounded Gauss upwind;
    div(phi,K)          bounded Gauss upwind;
-    div(phiv,p)     Gauss upwind;
+    div(phi,(p|rho))    Gauss upwind;
    div(phi,k)          bounded Gauss upwind;
    div(phi,epsilon)    bounded Gauss upwind;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/porous/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/porous/fvSchemes
@ -27,13 +27,15 @@ gradSchemes
 divSchemes
 {
    default             none;
+
    div(phi,U)          bounded Gauss upwind;
    div(phi,h)          bounded Gauss upwind;
    div(phi,e)          bounded Gauss upwind;
-    div(phiv,p)     Gauss upwind;
+    div(phi,(p|rho))    bounded Gauss upwind;
    div(phi,K)          bounded Gauss upwind;
    div(phi,k)          bounded Gauss upwind;
    div(phi,epsilon)    bounded Gauss upwind;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }

--- a/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/system/gas/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/system/gas/fvSchemes
@ -34,7 +34,8 @@ divSchemes
    div(phi,k)          Gauss linearUpwind limitedGrad;
    div(phi,epsilon)    Gauss linearUpwind limitedGrad;
    div(phi,K)          Gauss linear;
-    div(phiv,p)     Gauss linear;
+    div(phi,(p|rho))    Gauss linear;
+
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }