porousExplicitSourceReactingParcelFoam: Add total energy source term to enthalpy equation

Pressure time derivative term is run-time switchable

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -10,6 +10,8 @@
       + sources(rho, U)
     );
 
+    UEqn.relax();
+
     sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
@@ -17,5 +19,6 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -35,6 +35,8 @@ if (solveSpecies)
               + sources(rho, Yi)
             );
 
+            YiEqn.relax();
+
             sources.constrain(YiEqn);
 
             YiEqn.solve(mesh.solver("Yi"));

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -98,3 +98,9 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -1,47 +1,24 @@
 {
-    tmp<volScalarField> pWork
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "pWork",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
-        )
-    );
-
-    if (pressureWork)
-    {
-        surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-
-        pWork() += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-
-        if (pressureWorkTimeDerivative)
-        {
-            pWork() += fvc::ddt(p);
-        }
-    }
-
-    {
     fvScalarMatrix hsEqn
     (
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-            pWork()
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + parcels.Sh(hs)
       + radiation->Shs(thermo)
       + combustion->Sh()
       + sources(rho, hs)
     );
 
+    if (pressureWorkTimeDerivative)
+    {
+        hsEqn -= dpdt;
+    }
+
+    hsEqn.relax();
+
     sources.constrain(hsEqn);
 
     hsEqn.solve();
@@ -52,5 +29,4 @@
 
     Info<< "T gas min/max   = " << min(T).value() << ", "
         << max(T).value() << endl;
-    }
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -59,4 +59,7 @@
 
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    K = 0.5*magSqr(U);
+
+    dpdt = fvc::ddt(p);
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -1,9 +1,8 @@
 dictionary additional = mesh.solutionDict().subDict("additional");
 
 // pressure work term for enthalpy equation
-bool pressureWork = additional.lookupOrDefault("pressureWork", true);
 bool pressureWorkTimeDerivative =
-    additional.lookupOrDefault("pressureWorkTimeDerivative", true);
+    additional.lookupOrDefault("pressureWorkTimeDerivative", false);
 
 // flag to activate solve transport for each specie (Y vector)
 bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -645,20 +645,18 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
                     phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
                 )
                *(
-                    fvc::interpolate(rA*rho.oldTime())
+                    fvc::interpolate(rA)
                    *(
                        coefft0*phiAbs.oldTime()
-                      /fvc::interpolate(rho.oldTime())
                      - coefft00*phiAbs.oldTime().oldTime()
-                      /fvc::interpolate(rho.oldTime().oldTime())
                     )
                   - (
                         fvc::interpolate
                         (
-                            rA*rho.oldTime()*
+                            rA*
                             (
-                                coefft0*U.oldTime()
+                                coefft0*rho.oldTime()*U.oldTime()
-                              - coefft00*U.oldTime().oldTime()
+                              - coefft00*rho.oldTime().oldTime()*U.oldTime().oldTime()
                             )
                         ) & mesh().Sf()
                     )

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary

@@ -21,40 +21,40 @@ FoamFile
     {
         type            wall;
         nFaces          172;
-        startFace       3334;
+        startFace       3294;
     }
     inlet
     {
         type            patch;
         nFaces          20;
-        startFace       3506;
+        startFace       3466;
     }
     outlet
     {
         type            patch;
         nFaces          20;
-        startFace       3526;
+        startFace       3486;
     }
     cycLeft_half0
     {
         type            cyclic;
-        nFaces          0;
+        nFaces          20;
-        startFace       3546;
+        startFace       3506;
         matchTolerance  0.0001;
         neighbourPatch  cycLeft_half1;
     }
     cycLeft_half1
     {
         type            cyclic;
-        nFaces          0;
+        nFaces          20;
-        startFace       3546;
+        startFace       3526;
         matchTolerance  0.0001;
         neighbourPatch  cycLeft_half0;
     }
     cycRight_half0
     {
         type            cyclic;
-        nFaces          0;
+        nFaces          20;
         startFace       3546;
         matchTolerance  0.0001;
         neighbourPatch  cycRight_half1;
@@ -62,8 +62,8 @@ FoamFile
     cycRight_half1
     {
         type            cyclic;
-        nFaces          0;
+        nFaces          20;
-        startFace       3546;
+        startFace       3566;
         matchTolerance  0.0001;
         neighbourPatch  cycRight_half0;
     }
@@ -71,7 +71,7 @@ FoamFile
     {
         type            empty;
         nFaces          3440;
-        startFace       3546;
+        startFace       3586;
     }
 )
 

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution

@@ -85,8 +85,7 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
+    pressureWorkTimeDerivative false;
-    pressureWorkTimeDerivative true;
 }
 
 relaxationFactors

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes

@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phi,K)      Gauss linear;
+    div(phi,K)      Gauss upwind;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution

@@ -83,19 +83,18 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
+    pressureWorkTimeDerivative false;
-    pressureWorkTimeDerivative true;
 }
 
 relaxationFactors
 {
     fields
     {
-        ".*Final"       1;
+        ".*"       1;
     }
     equations
     {
-        ".*Final"       1;
+        ".*"       1;
     }
 }
 

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phi,K)      Gauss linear;
+    div(phi,K)      Gauss upwind;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution

@@ -89,7 +89,6 @@ potentialFlow
 
 additional
 {
-    pressureWork    false;
     pressureWorkTimeDerivative false;
 }
 
@@ -97,11 +96,11 @@ relaxationFactors
 {
     fields
     {
-        ".*Final"       1;
+        ".*"       1;
     }
     equations
     {
-        ".*Final"       1;
+        ".*"       1;
     }
 }