From 7d218b5e0221b938f124357e01b9eeb4fc5affef Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 24 Nov 2011 15:26:59 +0000
Subject: [PATCH] porousExplicitSourceReactingParcelFoam: Add total energy
 source term to enthalpy equation Pressure time derivative term is run-time
 switchable

---
 .../UEqn.H                                    |  3 +
 .../YEqn.H                                    |  2 +
 .../createFields.H                            |  6 ++
 .../hsEqn.H                                   | 62 ++++++-------------
 .../pEqn.H                                    |  3 +
 .../readAdditionalSolutionControls.H          |  3 +-
 .../backwardDdtScheme/backwardDdtScheme.C     | 10 ++-
 .../filter/constant/polyMesh/boundary         | 22 +++----
 .../filter/system/fvSolution                  |  3 +-
 .../parcelInBox/system/fvSchemes              |  2 +-
 .../parcelInBox/system/fvSolution             |  7 +--
 .../verticalChannel/system/fvSchemes          |  2 +-
 .../verticalChannel/system/fvSolution         |  5 +-
 13 files changed, 57 insertions(+), 73 deletions(-)

diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index 9bb0de54b1c..28f6f81d9f0 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/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);
     }
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
index e333a14712a..60f6e1a55da 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/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"));
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index ab6b01f9cea..a656d4ea8c5 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/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));
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index 27a9f6c8b3d..fc88ba785e4 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -1,56 +1,32 @@
 {
-    tmp<volScalarField> pWork
+    fvScalarMatrix hsEqn
     (
-        new volScalarField
-        (
-            IOobject
-            (
-                "pWork",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
-        )
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
+      + parcels.Sh(hs)
+      + radiation->Shs(thermo)
+      + combustion->Sh()
+      + sources(rho, hs)
     );
 
-    if (pressureWork)
+    if (pressureWorkTimeDerivative)
     {
-        surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-
-        pWork() += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-
-        if (pressureWorkTimeDerivative)
-        {
-            pWork() += fvc::ddt(p);
-        }
+        hsEqn -= dpdt;
     }
 
-    {
-        fvScalarMatrix hsEqn
-        (
-            fvm::ddt(rho, hs)
-          + mvConvection->fvmDiv(phi, hs)
-          - fvm::laplacian(turbulence->alphaEff(), hs)
-         ==
-            pWork()
-          + parcels.Sh(hs)
-          + radiation->Shs(thermo)
-          + combustion->Sh()
-          + sources(rho, hs)
-        );
+    hsEqn.relax();
 
-        sources.constrain(hsEqn);
+    sources.constrain(hsEqn);
 
-        hsEqn.solve();
+    hsEqn.solve();
 
-        thermo.correct();
+    thermo.correct();
 
-        radiation->correct();
+    radiation->correct();
 
-        Info<< "T gas min/max   = " << min(T).value() << ", "
-            << max(T).value() << endl;
-    }
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 23645a1b7be..080faeccedd 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/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);
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H
index 8136af1ad3e..050c24fe7ad 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H
+++ b/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);
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index cb77f67e760..857dd7637cb 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/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()
-                              - coefft00*U.oldTime().oldTime()
+                                coefft0*rho.oldTime()*U.oldTime()
+                              - coefft00*rho.oldTime().oldTime()*U.oldTime().oldTime()
                             )
                         ) & mesh().Sf()
                     )
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary
index 2016cf671ba..c9826291111 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary
+++ b/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;
-        startFace       3546;
+        nFaces          20;
+        startFace       3506;
         matchTolerance  0.0001;
         neighbourPatch  cycLeft_half1;
     }
     cycLeft_half1
     {
         type            cyclic;
-        nFaces          0;
-        startFace       3546;
+        nFaces          20;
+        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;
-        startFace       3546;
+        nFaces          20;
+        startFace       3566;
         matchTolerance  0.0001;
         neighbourPatch  cycRight_half0;
     }
@@ -71,7 +71,7 @@ FoamFile
     {
         type            empty;
         nFaces          3440;
-        startFace       3546;
+        startFace       3586;
     }
 )
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
index 574d2ac6a41..672f614cf5e 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
@@ -85,8 +85,7 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative true;
+    pressureWorkTimeDerivative false;
 }
 
 relaxationFactors
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
index 9b46cd8f1cc..a74c78c1fcd 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
+++ b/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;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
index 2927bf0ea22..ff633dd7d5e 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
@@ -83,19 +83,18 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative true;
+    pressureWorkTimeDerivative false;
 }
 
 relaxationFactors
 {
     fields
     {
-        ".*Final"       1;
+        ".*"       1;
     }
     equations
     {
-        ".*Final"       1;
+        ".*"       1;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
index 7a96d901447..15fb830b801 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
+++ b/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;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
index f5ac901b11c..3aeefebadd7 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/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;
     }
 }
 
-- 
GitLab