ENH: Updated steady parcel solver

applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H

@@ -1,22 +1,23 @@
+if (chemistry.chemistry())
 {
     Info<< "Solving chemistry" << endl;
 
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
+    // update reaction rates
+    chemistry.calculate();
 
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
-        DimensionedField<scalar, volMesh> tc =
-            chemistry.tc()().dimensionedInternalField();
+        typedef DimensionedField<scalar, volMesh> dsfType;
 
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+        const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL);
+
+        const dsfType tk =
+            Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0));
+
+        const dsfType tc = chemistry.tc()().dimensionedInternalField();
+
+        kappa = tc/(tc + tk);
     }
     else
     {

applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H

@@ -43,6 +43,11 @@
         }
     }
 
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+
     // Second part of thermodynamic density update
     thermo.rho() += psi*p;
 

applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H

@@ -1,4 +1,4 @@
-Info<< "Reading chemistry properties\n" << endl;
+// Info<< "Reading chemistry properties\n" << endl;
 
 IOdictionary chemistryProperties
 (

applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H

@@ -40,6 +40,9 @@ Description
     );
 
     rhoEqn.solve();
+
+    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+        << endl;
 }
 
 // ************************************************************************* //