deepmodeling · OpenFOAMFans · May 16, 2023 · May 12, 2023
diff --git a/applications/solvers/df0DFoam/createFields.H b/applications/solvers/df0DFoam/createFields.H
@@ -2,7 +2,8 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 fluidThermo& thermo = *pThermo;
 
 thermo.validate(args.executable(), "ha");

diff --git a/applications/solvers/df0DFoam/df0DFoam.C b/applications/solvers/df0DFoam/df0DFoam.C
@@ -27,7 +27,8 @@ Description
 \*---------------------------------------------------------------------------*/
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 
 #ifdef USE_PYTORCH
 #include <pybind11/embed.h>

diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H
@@ -2,8 +2,10 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
-psiThermo& thermo = *pThermo;
+// psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// psiThermo& thermo = *pThermo;
+rhoThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
+rhoThermo& thermo = *pThermo;
 
 //move from creatFieldRefs.H to createFields.H
 //p needed to be created before e

diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C
@@ -32,7 +32,8 @@ Description
 
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 
 #ifdef USE_PYTORCH
 #include <pybind11/embed.h>
@@ -47,7 +48,8 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+// #include "psiThermo.H"
+#include "rhoThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"

diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H
@@ -2,7 +2,8 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 fluidThermo& thermo = *pThermo;
 thermo.validate(args.executable(), "ha");
 

diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C
@@ -35,7 +35,8 @@ Description
 
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 
 #ifdef USE_PYTORCH
 #include <pybind11/embed.h>

diff --git a/applications/solvers/dfSprayFoam/createFields.H b/applications/solvers/dfSprayFoam/createFields.H
@@ -1,7 +1,8 @@
 #include "readGravitationalAcceleration.H"
 
 Info<< "Reading thermophysical properties\n" << endl;
-fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 fluidThermo& thermo = *pThermo;
 thermo.validate(args.executable(), "ha");
 

diff --git a/applications/solvers/dfSprayFoam/dfSprayFoam.C b/applications/solvers/dfSprayFoam/dfSprayFoam.C
@@ -31,7 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 #include "turbulentFluidThermoModel.H"
 
 #ifdef USE_PYTORCH

diff --git a/applications/solvers/dfSprayFoam/pEqn.H b/applications/solvers/dfSprayFoam/pEqn.H
@@ -3,6 +3,10 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -93,6 +97,7 @@ else
 p.relax();
 
 // Recalculate density from the relaxed pressure
+thermo.correctRho(psi*p - psip0);
 rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);

diff --git a/src/dfCanteraMixture/CanteraMixture.H b/src/dfCanteraMixture/CanteraMixture.H
@@ -121,6 +121,15 @@ public:
             return CanteraGas_->density()/p;
         }
 
+        scalar rho
+        (
+            const scalar& p,
+            const scalar& T
+        ) const
+        {
+            return CanteraGas_->density();
+        }
+
         scalar mu
         (
             const scalar& p,

diff --git a/src/dfCanteraMixture/makeThermos.C b/src/dfCanteraMixture/makeThermos.C
@@ -24,9 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "basicThermo.H"
-#include "psiThermo.H"
+// #include "psiThermo.H"
+#include "rhoThermo.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -71,7 +73,8 @@ namespace Foam
     addThermoPhysicsThermo(fluidThermo, CThermo##Mixture##ThermoPhys)
 
 
-makeThermoPhysicsThermos(hePsiThermo, CanteraMixture, psiThermo);
+// makeThermoPhysicsThermos(hePsiThermo, CanteraMixture, psiThermo);
+makeThermoPhysicsThermos(heRhoThermo, CanteraMixture, rhoThermo);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

diff --git a/src/dfChemistryModel/dfChemistryModel.C b/src/dfChemistryModel/dfChemistryModel.C
@@ -67,9 +67,12 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     alpha_(const_cast<volScalarField&>(thermo.alpha())),
     T_(thermo.T()),
     p_(thermo.p()),
-    rho_(mesh_.objectRegistry::lookupObject<volScalarField>("rho")),
-    mu_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).mu()())),
-    psi_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).psi())),
+    // rho_(mesh_.objectRegistry::lookupObject<volScalarField>("rho")),
+    rho_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).rho())),
+    // mu_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).mu()())),
+    // psi_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).psi())),
+    mu_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).mu()())),
+    psi_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).psi())),
     Qdot_
     (
         IOobject
@@ -385,6 +388,8 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
 
         psi_[celli] = mixture_.psi(p_[celli],T_[celli]);
 
+        rho_[celli] = mixture_.rho(p_[celli],T_[celli]);
+
         mu_[celli] = mixture_.CanteraTransport()->viscosity(); // Pa-s
 
         alpha_[celli] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass()); // kg/(m*s)
@@ -418,7 +423,7 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
 
     const volScalarField::Boundary& pBf = p_.boundaryField();
 
-    const volScalarField::Boundary& rhoBf = rho_.boundaryField();
+    volScalarField::Boundary& rhoBf = rho_.boundaryFieldRef();
 
     volScalarField::Boundary& TBf = T_.boundaryFieldRef();
 
@@ -433,7 +438,7 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
     forAll(T_.boundaryField(), patchi)
     {
         const fvPatchScalarField& pp = pBf[patchi];
-        const fvPatchScalarField& prho = rhoBf[patchi];
+        fvPatchScalarField& prho = rhoBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& ph = hBf[patchi];
@@ -454,6 +459,8 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
 
                 ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
 
+                prho[facei] = mixture_.rho(pp[facei],pT[facei]);
+
                 pmu[facei] = mixture_.CanteraTransport()->viscosity();
 
                 palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
@@ -495,6 +502,8 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
 
                 ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
 
+                prho[facei] = mixture_.rho(pp[facei],pT[facei]);
+
                 pmu[facei] = mixture_.CanteraTransport()->viscosity();
 
                 palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());

diff --git a/src/dfChemistryModel/dfChemistryModel.H b/src/dfChemistryModel/dfChemistryModel.H
@@ -62,7 +62,8 @@ SourceFiles
 #include "scalarField.H"
 #include "volFields.H"
 #include "hashedWordList.H"
-#include "psiThermo.H"
+// #include "psiThermo.H"
+#include "rhoThermo.H"
 #include "physicoChemicalConstants.H" // for R
 #include "ChemistryProblem.H"
 #include "ChemistrySolution.H"
@@ -121,7 +122,7 @@ public IOdictionary
         volScalarField& alpha_;
         volScalarField& T_;
         const volScalarField& p_;
-        const volScalarField& rho_;
+        volScalarField& rho_;
         volScalarField& mu_;
         volScalarField& psi_;
         // heat release rate, [J/m^3/s]

diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGM.C b/src/dfCombustionModels/FGM/baseFGM/baseFGM.C
@@ -44,7 +44,8 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     combustion_(this->coeffs().lookupOrDefault("combustion", false)),
     solveEnthalpy_(this->coeffs().lookupOrDefault("solveEnthalpy", false)),
     flameletT_(this->coeffs().lookupOrDefault("flameletT", false)),
-    psi_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).psi())),
+    // psi_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).psi())),
+    psi_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).psi())),
     Wt_ 
     (
         IOobject

diff --git a/src/dfCombustionModels/PaSR/PaSR.C b/src/dfCombustionModels/PaSR/PaSR.C
@@ -37,7 +37,8 @@ Foam::combustionModels::PaSR<ReactionThermo>::PaSR
 )
 :
     laminar<ReactionThermo>(modelType, thermo, turb, combustionProperties),
-    mu_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).mu()())),
+    // mu_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).mu()())),
+    mu_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).mu()())),
     p_(this->thermo().p()),
     T_(this->thermo().T()),
     mixingScaleDict_(this->coeffs().subDict("mixingScale")),