ENH: Updated LTSReactingParcelFoam solver

applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
+        #include "readPIMPLEControls.H"
         #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
@@ -86,17 +86,27 @@ int main(int argc, char *argv[])
         #include "timeScales.H"
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
-        #include "YEqn.H"
-        #include "hsEqn.H"
 
-        // --- PISO loop
+        // --- PIMPLE loop
-        for (int corr=0; corr<nCorr; corr++)
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
-            #include "pEqn.H"
+            bool finalIter = oCorr == nOuterCorr-1;
-        }
+            #include "addFinalIter.H"
 
-        turbulence->correct();
+            turbulence->correct();
+
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "hsEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            #include "removeFinalIter.H"
+        }
 
         if (runTime.write())
         {

applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H

@@ -1,4 +1,3 @@
-
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
@@ -10,7 +9,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     )
 );
 
-
 if (solveSpecies)
 {
     label inertIndex = -1;

applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H

@@ -73,12 +73,12 @@
 
     dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("PISO").lookup("rhoMax")
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax")
     );
 
     dimensionedScalar rhoMin
     (
-        mesh.solutionDict().subDict("PISO").lookup("rhoMin")
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
@@ -116,11 +116,11 @@
         dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );
 
-    volScalarField invTauFlow
+    volScalarField invTau
     (
         IOobject
         (
-            "invTauFlow",
+            "invTau",
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -40,7 +40,20 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    (
+                        finalIter
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H

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

applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H

@@ -1,10 +1,13 @@
 DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0);
 
-if (eWork)
+if (pressureWork)
 {
-    DpDt += -p*fvc::div(phi/fvc::interpolate(rho));
+    surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-}
+
-if (hWork)
+    DpDt += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-{
+
-    DpDt += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+    if (pressureWorkTimeDerivative)
+    {
+        DpDt += fvc::ddt(p);
+    }
 }

applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H

@@ -28,129 +28,40 @@ Info<< "Time scales min/max:" << endl;
 
 {
     // Cache old time scale field
-    tmp<volScalarField> tinvTauFlow0
+    tmp<volScalarField> tinvTau0
     (
         new volScalarField
         (
             IOobject
             (
-                "invTauFlow0",
+                "invTau0",
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            invTauFlow
+            invTau
         )
     );
-    const volScalarField& invTauFlow0 = tinvTauFlow0();
+    const volScalarField& invTau0 = tinvTau0();
 
 
     // Flow time scale
     // ~~~~~~~~~~~~~~~
     {
-        invTauFlow =
+        invTau =
             fvc::surfaceSum
             (
                 mag(phi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
             )
            /rho;
 
-        invTauFlow.max(1.0/maxDeltaT);
+        invTau.max(1.0/maxDeltaT);
 
         Info<< "    Flow        = "
-            << gMin(1/invTauFlow.internalField()) << ", "
+            << gMin(1/invTau.internalField()) << ", "
-            << gMax(1/invTauFlow.internalField()) << endl;
+            << gMax(1/invTau.internalField()) << endl;
-    }
-
-
-    // Mass source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-        scalarField tau
-        (
-            runTime.deltaTValue()*mag(parcels.Srho() + massSource.SuTot())
-        );
-
-        tau = alphaTauRho*rho/(tau + ROOTVSMALL);
-
-        Info<< "    Density     = "
-            << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
-    }
-
-
-    // Momentum source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-/*
-        // Method 1 - mag(U) limit using 'small' nominal velocity
-        scalarField tau
-        (
-            runTime.deltaTValue()
-           *mag
-            (
-                rho.dimensionedInternalField()*g
-              + parcels.UTrans()/(mesh.V()*runTime.deltaT())
-              + momentumSource.Su()
-            )
-           /rho
-        );
-
-        const scalar nomMagU(dimensionedScalar("1", dimVelocity, 1));
-        tau = alphaTauU*(nomMagU + mag(U))/(tau + ROOTVSMALL);
-*/
-/*
-        // Method 2 - based on fluxes and Co-like limit
-        volVectorField UEqnRhs
-        (
-            IOobject
-            (
-                "UEqnRhs",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedVector("zero", dimDensity*dimAcceleration, vector::zero),
-            zeroGradientFvPatchVectorField::typeName
-        );
-
-        UEqnRhs.internalField() =
-            rho.dimensionedInternalField()*g
-          + parcels.UTrans()/(mesh.V()*runTime.deltaT())
-          + momentumSource.Su();
-
-        surfaceScalarField phiSU
-        (
-            "phiSU",
-            fvc::interpolate(runTime.deltaT()*UEqnRhs) & mesh.Sf()
-        );
-
-        scalarField tau
-        (
-            alphaTauU*rho
-           /fvc::surfaceSum
-            (
-                mag(phi + phiSU)*mesh.deltaCoeffs()/mesh.magSf()
-              + dimensionedScalar("SMALL", dimDensity/dimTime, ROOTVSMALL)
-            )
-        );
-
-*/
-/*
-        Info<< "    Momentum    = " << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
-*/
     }
 
 
@@ -176,42 +87,7 @@ Info<< "Time scales min/max:" << endl;
         Info<< "    Temperature = " << min(maxDeltaT, gMin(tau)) << ", "
             << min(maxDeltaT, gMax(tau)) << endl;
 
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
+        invTau.internalField() = max(invTau.internalField(), 1/tau);
-    }
-
-
-    // Specie source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-        scalarField tau(mesh.nCells(), ROOTVGREAT);
-        forAll(Y, fieldI)
-        {
-            const volScalarField& Yi = Y[fieldI];
-            const scalarField deltaYi
-            (
-                runTime.deltaTValue()
-               *mag
-                (
-                    kappa*chemistry.RR(fieldI)()
-                  + massSource.Su(fieldI)
-                  + parcels.Srho(fieldI)
-                )
-               /rho
-            );
-            tau =
-                min
-                (
-                    tau,
-                    alphaTauSpecie
-                   /(deltaYi/(Yi + specieMaxUnbound) + ROOTVSMALL)
-                );
-        }
-
-        Info<< "    Specie      = " << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
     }
 
 
@@ -219,21 +95,21 @@ Info<< "Time scales min/max:" << endl;
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - reduce as much as required for flow, but limit source contributions
     const dimensionedScalar deltaTRamp("deltaTRamp", dimless, 1/(1 + 0.2));
-    invTauFlow = max(invTauFlow, invTauFlow0*deltaTRamp);
+    invTau = max(invTau, invTau0*deltaTRamp);
-    tinvTauFlow0.clear();
+    tinvTau0.clear();
 
 
     // Limit the largest time scale
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    invTauFlow.max(1/maxDeltaT);
+    invTau.max(1/maxDeltaT);
 
 
     // Spatially smooth the time scale field
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    fvc::smooth(invTauFlow, alphaTauSmooth);
+    fvc::smooth(invTau, alphaTauSmooth);
 
-    Info<< "    Overall     = " << min(1/invTauFlow).value()
+    Info<< "    Overall     = " << min(1/invTau).value()
-        << ", " << max(1/invTauFlow).value() << nl << endl;
+        << ", " << max(1/invTau).value() << nl << endl;
 }
 
 

src/finiteVolume/cfdTools/general/include/addFinalIter.H

@@ -0,0 +1,5 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+

src/finiteVolume/cfdTools/general/include/removeFinalIter.H

@@ -0,0 +1,4 @@
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}