updates

applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H

     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+//        fvm::ddt(rho, U)
+        pZones.ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
      ==
@@ -32,13 +33,13 @@
     {
         pZones.addResistance(UEqn);
 
+        trAU = 1.0/UEqn.A();
+        trAU().rename("rAU");
+
         solve
         (
             UEqn == -fvc::grad(p)
         );
-
-        trAU = 1.0/UEqn.A();
-        trAU().rename("rAU");
     }
 
     if (momentumPredictor)

applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H

@@ -28,7 +28,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 parcels.Srho(i)
               + kappa*chemistry.RR(i)().dimensionedInternalField()
-              + pointMassSources.Su(i)
+              + pointMassSources.Su(i),
+                mesh.solver("Yi")
             );
 
             Yi.max(0.0);

applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H

 Info<< "\nConstructing reacting cloud" << endl;
-BasicTrackedReactingCloud<gasThermoPhysics> parcels
+BasicTrackedReactingCloud<icoPoly8ThermoPhysics> parcels
 (
     "reactingCloud1",
     rho,

applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H

-rho = thermo.rho();
-
-if (pressureImplicitPorosity)
-{
-    U = trTU()&UEqn.H();
-}
-else
 {
-    U = trAU()*UEqn.H();
-}
+    rho = thermo.rho();
 
-if (transonic)
-{
-    surfaceScalarField phid
-    (
-        "phid",
-        fvc::interpolate(psi)
-       *(
-            (fvc::interpolate(U) & mesh.Sf())
-//          + fvc::ddtPhiCorr(rUA, rho, U, phi)
-        )
-    );
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    if (pressureImplicitPorosity)
+    {
+        U = trTU()&UEqn.H();
+    }
+    else
     {
-        tmp<fvScalarMatrix> lapTerm;
+        U = trAU()*UEqn.H();
+    }
 
-        if (pressureImplicitPorosity)
-        {
-            lapTerm = fvm::laplacian(rho*trTU(), p);
-        }
-        else
-        {
-            lapTerm = fvm::laplacian(rho*trAU(), p);
-        }
+    if (transonic)
+    {
+        surfaceScalarField phiv = fvc::interpolate(U) & mesh.Sf();
 
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - lapTerm()
-         ==
-            parcels.Srho()
-          + pointMassSources.Su()
-        );
+        phi = fvc::interpolate(rho)*phiv;
 
-        pEqn.solve();
+        surfaceScalarField phid("phid", fvc::interpolate(psi)*phiv);
 
-        if (nonOrth == nNonOrthCorr)
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
-            phi == pEqn.flux();
+            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
-{
-    phi =
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(trAU(), rho, U, phi)
-        );
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    else
     {
-        tmp<fvScalarMatrix> lapTerm;
-
-        if (pressureImplicitPorosity)
-        {
-            lapTerm = fvm::laplacian(rho*trTU(), p);
-        }
-        else
+        phi =
+            fvc::interpolate(rho)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+//            + fvc::ddtPhiCorr(trAU(), rho, U, phi)
+            );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
-            lapTerm = fvm::laplacian(rho*trAU(), p);
+            tmp<fvScalarMatrix> lapTerm;
+
+            if (pressureImplicitPorosity)
+            {
+                lapTerm = fvm::laplacian(rho*trTU(), p);
+            }
+            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();
+            }
         }
+    }
 
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phi)
-          - lapTerm()
-         ==
-            parcels.Srho()
-          + pointMassSources.Su()
-        );
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
 
-        pEqn.solve();
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
 
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi += pEqn.flux();
-        }
+    if (pressureImplicitPorosity)
+    {
+        U -= trTU()&fvc::grad(p);
     }
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
+    else
+    {
+        U -= trAU()*fvc::grad(p);
+    }
+    U.correctBoundaryConditions();
 
-if (pressureImplicitPorosity)
-{
-    U -= trTU()&fvc::grad(p);
+    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 }
-else
-{
-    U -= trAU()*fvc::grad(p);
-}
-U.correctBoundaryConditions();
-
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-

applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve
+    fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
@@ -39,6 +39,10 @@ Description
         parcels.Srho()
       + pointMassSources.Su()
     );
+
+    rhoEqn.relax();
+
+    rhoEqn.solve();
 }
 
 // ************************************************************************* //

applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C

@@ -23,16 +23,23 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
+    trackedReactingParcelFoam
 
 Description
+    - reacting parcel cloud tracking
+    - porous media
+    - point mass sources
+    - polynomial based, incompressible thermodynamics (f(T))
+
+    Note: ddtPhiCorr not used here - not well defined for porous calcs
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
+#include "hReactionThermo.H"
 #include "turbulenceModel.H"
 #include "BasicTrackedReactingCloud.H"
-#include "psiChemistryModel.H"
+#include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
 #include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
@@ -81,17 +88,17 @@ int main(int argc, char *argv[])
 
         #include "chemistry.H"
         #include "rhoEqn.H"
+        #include "UEqn.H"
 
         // --- PIMPLE loop
-        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        for (int oCorr=1; oCorr<=nOuterCorr; oCorr++)
         {
-            #include "UEqn.H"
             #include "YEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
                 #include "pEqn.H"
             }