fix energy equation and change convection term in dfHighSpeedFoam, add new example for dfHighSpeedFoam (#287)

* fix energy equation and change convection term in dfHighSpeedFoam, add new example for dfHighSpeedFoam

* Update dfLowMachFoam.C

* fix dfLowMachFoam and dfSprayFoam

* fix df0DFoam and dfSprayFoam

* fix df0DFoam and dfSprayFoam

* Update CanteraMixture.C, fix name

* fix dfLowMachFoam, df0DFoam and dfSprayFoam

Author: pkuLmq

applications/solvers/df0DFoam/createFields.H

@@ -6,7 +6,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 fluidThermo& thermo = *pThermo;
 
-thermo.validate(args.executable(), "ha");
+// thermo.validate(args.executable(), "ha");
 
 volScalarField& p = thermo.p();
 
@@ -94,6 +94,8 @@ dfChemistryModel<basicThermo> chemistry(thermo);
 PtrList<volScalarField>& Y = chemistry.Y();
 const word inertSpecie(chemistry.lookup("inertSpecie"));
 const label inertIndex(chemistry.species()[inertSpecie]);
+chemistry.setEnergyName("ha");
+chemistry.updateEnergy();
 
 forAll(Y, i)
 {

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -29,7 +29,7 @@ dictionary thermoDict
     )
 );
 
-bool inviscid(thermoDict.lookup("inviscid"));
+bool inviscid(thermoDict.lookupOrDefault("inviscid",false));
 
 Info<< "Reading field U\n" << endl;
 volVectorField U
@@ -58,9 +58,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField e = thermo.he() - p/rho;
-volScalarField& ha = thermo.he();
-
 volVectorField rhoU
 (
     IOobject
@@ -74,19 +71,6 @@ volVectorField rhoU
     rho*U
 );
 
-volScalarField rhoE
-(
-    IOobject
-    (
-        "rhoE",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    rho*(e + 0.5*magSqr(U))
-);
-
 surfaceScalarField pos
 (
     IOobject
@@ -140,6 +124,22 @@ const label inertIndex(chemistry->species()[inertSpecie]);
 chemistry->setEnergyName("ea");
 chemistry->updateEnergy();
 
+volScalarField& ea = thermo.he();
+volScalarField ha = ea + p/rho;
+
+volScalarField rhoE
+(
+    IOobject
+    (
+        "rhoE",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    rho*(ea + 0.5*magSqr(U))
+);
+
 chemistry->correctThermo();
 Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
 
@@ -149,6 +149,28 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
+const label nspecies(chemistry->species().size());
+PtrList<volScalarField> rhoYi(nspecies);
+forAll(rhoYi,i)
+{
+    rhoYi.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoYi" + Y[i].name(),
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            rho*Y[i]
+        )
+    );
+}
+
 const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
 volScalarField diffAlphaD
 (

applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C

@@ -132,15 +132,57 @@ int main(int argc, char *argv[])
         surfaceScalarField rho_pos(interpolate(rho, pos));
         surfaceScalarField rho_neg(interpolate(rho, neg));
 
+        PtrList<surfaceScalarField> rhoYi_pos(nspecies);
+        PtrList<surfaceScalarField> rhoYi_neg(nspecies);
+        forAll(rhoYi_pos,i)
+        {
+            rhoYi_pos.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "rhoYi_pos" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    interpolate(rhoYi[i], pos,"Yi")
+                )
+            );
+        }
+
+        forAll(rhoYi_neg,i)
+        {
+            rhoYi_neg.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "rhoYi_neg" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    interpolate(rhoYi[i], neg,"Yi")
+                )
+            );
+        }
+
         surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
         surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
 
         volScalarField rPsi("rPsi", 1.0/psi);
         surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
         surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
 
-        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
-        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
+        surfaceScalarField ea_pos(interpolate(ea, pos, T.name()));
+        surfaceScalarField ea_neg(interpolate(ea, neg, T.name()));
 
         surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
         surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
@@ -217,6 +259,27 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
+        PtrList<surfaceScalarField> phiYi(nspecies);
+        forAll(phiYi,i)
+        {
+            phiYi.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "phiYi_" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i]
+                )
+            );
+        }
+
         surfaceVectorField phiUp
         (
             (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
@@ -226,8 +289,8 @@ int main(int argc, char *argv[])
         surfaceScalarField phiEp
         (
             "phiEp",
-            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
-          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+            aphiv_pos*(rho_pos*(ea_pos + 0.5*magSqr(U_pos)) + p_pos)
+          + aphiv_neg*(rho_neg*(ea_neg + 0.5*magSqr(U_neg)) + p_neg)
           + aSf*p_pos - aSf*p_neg
         );
 

applications/solvers/dfHighSpeedFoam/rhoEEqn.H

@@ -15,31 +15,31 @@ solve
     - fvc::div(sigmaDotU)
 );
 
-e = rhoE/rho - 0.5*magSqr(U);
-e.correctBoundaryConditions();
+ea = rhoE/rho - 0.5*magSqr(U);
+ea.correctBoundaryConditions();
 
-ha = e + p/rho;
+ha = ea + p/rho;
 chemistry->correctThermo(); // before this, we must update ha = e + p/rho
 
-rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField() + 0.5*magSqr(U.boundaryField()));
+rhoE.boundaryFieldRef() == rho.boundaryField()*(ea.boundaryField() + 0.5*magSqr(U.boundaryField()));
 
 if (!inviscid)
 {
     fvScalarMatrix eEqn
     (
-           fvm::ddt(rho, e) - fvc::ddt(rho, e)
+           fvm::ddt(rho, ea) - fvc::ddt(rho, ea)
         // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha
-         - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), e)
+         - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
          + diffAlphaD
       ==
            fvc::div(hDiffCorrFlux)
     );
 
-    eEqn.solve("e");
+    eEqn.solve("ea");
 
-    ha = e + p/rho;
+    ha = ea + p/rho;
     chemistry->correctThermo();
-    rhoE = rho*(e + 0.5*magSqr(U));
+    rhoE = rho*(ea + 0.5*magSqr(U));
 }
 
 Info<< "min/max(T) = "

applications/solvers/dfHighSpeedFoam/rhoYEqn.H

@@ -37,32 +37,63 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     {
         volScalarField& Yi = Y[i];
 
+        if (!inviscid)
+        {
+            hDiffCorrFlux += chemistry->hai(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError);
+            diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hai(i), Yi);
+        }
+
         if (i != inertIndex)
         {          
-            fvScalarMatrix YiEqn
+            // original convection term
+            // fvScalarMatrix YiEqn
+            // (
+            //      fvm::ddt(rho, Yi)
+            //    + mvConvection->fvmDiv(phi, Yi)
+            //   ==
+            //     combustion->R(Yi)
+            // );
+
+            solve
             (
-                 fvm::ddt(rho, Yi)
-               + mvConvection->fvmDiv(phi, Yi)
-              ==
-                combustion->R(Yi)
+                  fvm::ddt(rhoYi[i])
+                + fvc::div(phiYi[i])
+               == 
+                  chemistry->RR(i)
             );
 
+            Yi=rhoYi[i]/rho;
+            Yi.max(0.0);
+
             if (!inviscid)
             {
                 const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
-
-                hDiffCorrFlux += chemistry->hai(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError);
-                diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hai(i), Yi);
                 tmp<volScalarField> DEff = chemistry->rhoD(i) + turbulence->mut()/Sct;
 
-                YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi);
-            }
+                // original term
+                // YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi);
 
-            YiEqn.relax();
+                fvScalarMatrix YiEqn
+                (
+                    fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi)
+                  - fvm::laplacian(DEff(), Yi)
+                  + mvConvection->fvmDiv(phiUc, Yi)
+                );
 
-            YiEqn.solve("Yi");
+                YiEqn.relax();
 
-            Yi.max(0.0);
+                YiEqn.solve("Yi");
+
+                Yi.max(0.0);
+
+            }
+
+            // original term
+            // YiEqn.relax();
+            // YiEqn.solve("Yi");
+            // Yi.max(0.0);
+
+            rhoYi[i] = rho*Yi;
             Yt += Yi;
         }
     }

applications/solvers/dfLowMachFoam/EEqn.H

@@ -25,7 +25,7 @@
 
         EEqn.relax();
 
-        EEqn.solve();
+        EEqn.solve("ha");
 
 
 }

applications/solvers/dfLowMachFoam/createFields.H

@@ -5,7 +5,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 // 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");
+// thermo.validate(args.executable(), "ha");
 
 const volScalarField& psi = thermo.psi();
 volScalarField& p = thermo.p();
@@ -89,6 +89,8 @@ dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
 PtrList<volScalarField>& Y = chemistry->Y();
 const word inertSpecie(chemistry->lookup("inertSpecie"));
 const label inertIndex(chemistry->species()[inertSpecie]);
+chemistry->setEnergyName("ha");
+chemistry->updateEnergy();
 
 
 chemistry->correctThermo();

applications/solvers/dfLowMachFoam/dfLowMachFoam.C

@@ -196,7 +196,7 @@ int main(int argc, char *argv[])
         rho = thermo.rho();
 
         runTime.write();
-
+      
         Info<< "========Time Spent in diffenet parts========"<< endl;
         Info<< "Chemical sources           = " << time_monitor_chem << " s" << endl;
         Info<< "Species Equations          = " << time_monitor_Y << " s" << endl;

applications/solvers/dfSprayFoam/EEqn.H

@@ -48,7 +48,7 @@
 
     //fvOptions.constrain(EEqn);
 
-    EEqn.solve();
+    EEqn.solve("ha");
 
     //fvOptions.correct(he);
 

applications/solvers/dfSprayFoam/createFields.H

@@ -4,7 +4,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 // 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");
+// thermo.validate(args.executable(), "ha");
 
 SLGThermo slgThermo(mesh, thermo);
 
@@ -88,6 +88,8 @@ dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
 PtrList<volScalarField>& Y = chemistry->Y();
 const word inertSpecie(chemistry->lookup("inertSpecie"));
 const label inertIndex(chemistry->species()[inertSpecie]);
+chemistry->setEnergyName("ha");
+chemistry->updateEnergy();
 
 
 chemistry->correctThermo();