From 74499d57d9f41e9c98f8d099bb767b9a9f2e5197 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Thu, 28 Oct 2010 18:14:57 +0100
Subject: [PATCH] ENH: Updated steady parcel solver

---
 .../steadyReactingParcelFoam/chemistry.H      | 23 ++++++++++---------
 .../steadyReactingParcelFoam/pEqn.H           |  5 ++++
 .../readChemistryProperties.H                 |  2 +-
 .../steadyReactingParcelFoam/rhoEqn.H         |  3 +++
 4 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
index 50fc7f575b5..ec92091de37 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
+++ b/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
     {
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
index 6ae2712ea42..e7d7c55dbbe 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
+++ b/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;
 
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
index f0bcf7597fc..e742e9fea78 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
@@ -1,4 +1,4 @@
-Info<< "Reading chemistry properties\n" << endl;
+// Info<< "Reading chemistry properties\n" << endl;
 
 IOdictionary chemistryProperties
 (
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
index 33b7943cadd..e42db723993 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
@@ -40,6 +40,9 @@ Description
     );
 
     rhoEqn.solve();
+
+    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+        << endl;
 }
 
 // ************************************************************************* //
-- 
GitLab