updates to ../../applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -9,35 +9,10 @@
       + parcels.SU()
     );
 
-    UEqn.relax();
+    pZones.addResistance(UEqn);
 
-    tmp<volScalarField> trAU;
+    if (momentumPredictor)
-    tmp<volTensorField> trTU;
-
-    if (pressureImplicitPorosity)
     {
-        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
+        solve(UEqn == -fvc::grad(p));
-        pZones.addResistance(UEqn, tTU());
-        trTU = inv(tTU());
-        trTU().rename("rAU");
-
-        volVectorField gradp = fvc::grad(p);
-
-        for (int UCorr=0; UCorr<nUCorr; UCorr++)
-        {
-            U = trTU() & (UEqn.H() - gradp);
-        }
-        U.correctBoundaryConditions();
     }
-    else
-    {
-        pZones.addResistance(UEqn);
 
-        trAU = 1.0/UEqn.A();
-        trAU().rename("rAU");
-
-        solve
-        (
-            UEqn == -fvc::grad(p)
-        );
-    }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -43,5 +43,4 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
-
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -73,13 +73,7 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
-    volScalarField DpDt
-    (
-        "DpDt",
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
-
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll (Y, i)

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H

@@ -1,21 +1,3 @@
     Info<< "Creating porous zones" << nl << endl;
 
     porousZones pZones(mesh);
-    Switch pressureImplicitPorosity(false);
-
-    label nUCorr = 0;
-    if (pZones.size())
-    {
-        // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
-        {
-            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
-                >> nUCorr;
-        }
-
-        if (nUCorr > 0)
-        {
-            pressureImplicitPorosity = true;
-        }
-    }
-

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H

@@ -1,20 +1,51 @@
 {
-    fvScalarMatrix hEqn
+    tmp<volScalarField> pWork
     (
-        fvm::ddt(rho, h)
+        new volScalarField
-      + mvConvection->fvmDiv(phi, h)
+        (
-      - fvm::laplacian(turbulence->alphaEff(), h)
+            IOobject
-     ==
+            (
-        DpDt
+                "pWork",
-     +  parcels.Sh()
+                runTime.timeName(),
-     +  radiation->Sh(thermo)
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
+        )
     );
 
-    hEqn.relax();
+    if (dpdt)
+    {
+        pWork() += fvc::ddt(p);
+    }
+    if (eWork)
+    {
+        pWork() = -p*fvc::div(phi/fvc::interpolate(rho));
+    }
+    if (hWork)
+    {
+        pWork() += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+    }
 
-    hEqn.solve();
+    {
+        solve
+        (
+            fvm::ddt(rho, h)
+          + mvConvection->fvmDiv(phi, h)
+          - fvm::laplacian(turbulence->alphaEff(), h)
+         ==
+            pWork()
+          + parcels.Sh()
+          + radiation->Sh(thermo)
+        );
 
-    thermo.correct();
+        thermo.correct();
 
-    radiation->correct();
+        radiation->correct();
+
+        Info<< "T gas min/max   = " << min(T).value() << ", "
+            << max(T).value() << endl;
+    }
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -5,65 +5,13 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    if (pressureImplicitPorosity)
+    volScalarField rAU = 1.0/UEqn.A();
+    U = rAU*UEqn.H();
+
+    if (pZones.size() > 0)
     {
-        U = trTU()&UEqn.H();
+        // ddtPhiCorr not well defined for cases with porosity
-    }
+        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    else
-    {
-        U = trAU()*UEqn.H();
-    }
-
-    if (transonic)
-    {
-        surfaceScalarField phiv = fvc::interpolate(U) & mesh.Sf();
-
-        phi = fvc::interpolate(rho)*phiv;
-
-        surfaceScalarField phid("phid", fvc::interpolate(psi)*phiv);
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            tmp<fvScalarMatrix> lapTerm;
-
-            if (pressureImplicitPorosity)
-            {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
-            fvScalarMatrix pEqn
-            (
-                fvc::ddt(rho) + fvc::div(phi)
-              + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-              - lapTerm()
-            ==
-                parcels.Srho()
-              + pointMassSources.Su()
-            );
-
-            if
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solver("pFinal"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi == pEqn.flux();
-            }
-        }
     }
     else
     {
@@ -71,50 +19,34 @@
             fvc::interpolate(rho)
            *(
                 (fvc::interpolate(U) & mesh.Sf())
-//            + fvc::ddtPhiCorr(trAU(), rho, U, phi)
+              + fvc::ddtPhiCorr(rAU, rho, U, phi)
             );
+    }
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rAU, p)
+        ==
+            parcels.Srho()
+          + pointMassSources.Su()
+        );
+
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            tmp<fvScalarMatrix> lapTerm;
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
 
-            if (pressureImplicitPorosity)
+        if (nonOrth == nNonOrthCorr)
-            {
+        {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
+            phi += pEqn.flux();
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
-            fvScalarMatrix pEqn
-            (
-                fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-              + fvc::div(phi)
-              - lapTerm()
-            ==
-                parcels.Srho()
-              + pointMassSources.Su()
-            );
-
-            if
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solver("pFinal"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi += pEqn.flux();
-            }
         }
     }
 
@@ -124,15 +56,6 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    if (pressureImplicitPorosity)
+    U -= rAU*fvc::grad(p);
-    {
-        U -= trTU()&fvc::grad(p);
-    }
-    else
-    {
-        U -= trAU()*fvc::grad(p);
-    }
     U.correctBoundaryConditions();
-
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -23,15 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    trackedReactingParcelFoam
+    porousExplicitSourceReactingParcelFoam
 
 Description
-    - reacting parcel cloud tracking
+    - reacting parcel cloud
     - porous media
     - point mass sources
     - polynomial based, incompressible thermodynamics (f(T))
 
-    Note: ddtPhiCorr not used here - not well defined for porous calcs
+    Note: ddtPhiCorr not used here when porous zones are active
+    - not well defined for porous calcs
 
 \*---------------------------------------------------------------------------*/
 
@@ -74,6 +75,7 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "readPISOControls.H"
+        #include "readAdditionalSolutionControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -88,23 +90,16 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
+        #include "YEqn.H"
+        #include "hEqn.H"
 
-        // --- PIMPLE loop
+        // --- PISO loop
-        for (int oCorr=1; oCorr<=nOuterCorr; oCorr++)
+        for (int corr=0; corr<nCorr; corr++)
         {
-            #include "YEqn.H"
+            #include "pEqn.H"
-            #include "hEqn.H"
-
-            // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
-            {
-                #include "pEqn.H"
-            }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
+
         turbulence->correct();
 
         rho = thermo.rho();

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -0,0 +1,20 @@
+dictionary additional = mesh.solutionDict().subDict("additional");
+
+bool dpdt = true;
+if (additional.found("dpdt"))
+{
+    additional.lookup("dpdt") >> dpdt;
+}
+
+bool eWork = true;
+if (additional.found("eWork"))
+{
+    additional.lookup("eWork") >> eWork;
+}
+
+bool hWork = true;
+if (additional.found("hWork"))
+{
+    additional.lookup("hWork") >> hWork;
+}
+