From d9a9362ca607ddc7d364cab9e038d60d9ad29501 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sun, 30 Sep 2018 21:18:38 +0900
Subject: [PATCH] slfm with energy, working.

---
 .gitignore           |  7 +++++++
 EEqn.H               | 42 ++++++++++++++++++++++++++++++++++++++++++
 Make/files           |  2 +-
 Make/options         |  2 ++
 MixturefractionVar.H | 12 ++++++++++--
 SLFMFoam.C           |  9 ++++-----
 createFields.H       |  4 +++-
 pEqn.H               |  1 +
 setSDR.H             |  2 +-
 updateYi.H           |  5 +++++
 10 files changed, 76 insertions(+), 10 deletions(-)
 create mode 100644 .gitignore
 create mode 100644 EEqn.H

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e53befb
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+[._]*.s[a-v][a-z]
+[._]*.sw[a-p]
+[._]s[a-v][a-z]
+[._]sw[a-p]
+Make
+!Make/files
+!Make/options
diff --git a/EEqn.H b/EEqn.H
new file mode 100644
index 0000000..69e4026
--- /dev/null
+++ b/EEqn.H
@@ -0,0 +1,42 @@
+{
+    tmp<fv::convectionScheme<scalar> > mvConvection
+    (
+        fv::convectionScheme<scalar>::New
+        (
+            mesh,
+            fields,
+            phi,
+            mesh.divScheme("div(phi,Yi_h)")
+        )
+    );
+
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        mvConvection->fvmDiv(phi, he)
+      + (
+            he.name() == "ea"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        rho*(U&g)
+      + radiation->Sh(thermo)
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/Make/files b/Make/files
index 3349e8d..59a5f45 100644
--- a/Make/files
+++ b/Make/files
@@ -14,4 +14,4 @@ FlameStructure/FlameStructureIO.C
 
 SLFMFoam.C
 
-EXE = $(FOAM_APPBIN)/SLFMFoam
+EXE = $(FOAM_APPBIN)/slfmEnergyFoam
diff --git a/Make/options b/Make/options
index 3e21a2b..95b8843 100644
--- a/Make/options
+++ b/Make/options
@@ -13,6 +13,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
@@ -30,5 +31,6 @@ EXE_LIBS = \
     -lchemistryModel \
     -lODE \
     -lcombustionModels \
+    -lradiationModels \
     -lfvOptions \
     -lsampling
diff --git a/MixturefractionVar.H b/MixturefractionVar.H
index bcd78f1..bdc172b 100644
--- a/MixturefractionVar.H
+++ b/MixturefractionVar.H
@@ -1,6 +1,10 @@
 //mixture fraction variance equation
 //cf) mfVar should not be ZERO
-SDR = turbulence->epsilon() * mfVar / turbulence->k();
+
+Info<<"Calculating epsilon/k"<<endl;
+volScalarField epsk(turbulence->epsilon() / turbulence->k());
+
+Info<<"Calculating mixture fraction gradient"<<endl;
 volVectorField Gradmf = fvc::grad(mf);
 
 fvScalarMatrix mfVarEqn
@@ -9,9 +13,13 @@ fvScalarMatrix mfVarEqn
   - fvm::laplacian((1/Sc) * turbulence->mut(), mfVar)
  ==
     2 * ((1/Sc) * turbulence->mut()) * (Gradmf & Gradmf)
-  - 2 * rho * SDR
+  - fvm::Sp(2 * rho * epsk, mfVar)
 );
 
 mfVarEqn.relax();
 
 mfVarEqn.solve();
+
+mfVar.writeMinMax(Info);
+
+mfVar.max(0.0);
diff --git a/SLFMFoam.C b/SLFMFoam.C
index 8d1d489..a14bc2e 100644
--- a/SLFMFoam.C
+++ b/SLFMFoam.C
@@ -43,6 +43,7 @@ Contact
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "rhoCombustionModel.H"
+#include "radiationModel.H"
 #include "fvOptions.H"
 #include "simpleControl.H"
 #include "FlameStructure/FlameStructure.H"
@@ -93,16 +94,14 @@ int main(int argc, char *argv[])
             #include "Mixturefraction.H" //mean mixture fraction
             #include "MixturefractionVar.H" //mean mixture fraction variance
             #include "setSDR.H" //calculate scalar dissipation rate
-            #include "updateT_RHO.H" //update temperature and density
+            #include "updateYi.H" //update species
+            #include "EEqn.H" //update temperature and density
             #include "pEqn.H"
         }
 
         turbulence->correct();
 
-        if (runTime.write())
-        {
-            #include "updateYi.H" //update species
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/createFields.H b/createFields.H
index ea1c5ac..6eb2f20 100644
--- a/createFields.H
+++ b/createFields.H
@@ -8,7 +8,7 @@ autoPtr<combustionModels::rhoCombustionModel> combustion
 );
 
 rhoReactionThermo& thermo = combustion->thermo();
-thermo.validate(args.executable(), "h", "e");
+thermo.validate(args.executable(), "ha", "ea");
 
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -371,4 +371,6 @@ for(label j=0 ; j<etamax ; j++)
     );
 }
 
+#include "createRadiationModel.H"
+
 #include "createMRF.H"
diff --git a/pEqn.H b/pEqn.H
index e5f3980..1598c3c 100644
--- a/pEqn.H
+++ b/pEqn.H
@@ -47,6 +47,7 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    rho = thermo.rho();
     rho = max(rho, rhoMin);
     rho = min(rho, rhoMax);
     rho.relax();
diff --git a/setSDR.H b/setSDR.H
index c2e4168..bf84630 100644
--- a/setSDR.H
+++ b/setSDR.H
@@ -72,7 +72,7 @@ forAll(mf, cellI)
     }
     else
     {
-        Neta[lowerN][cellI] = SDR[cellI]*Coeff*C1coeff; //CSDR at stoichiometic m.f.
+        Neta[lowerN][cellI] = mfVar[cellI]*epsk[cellI]*Coeff*C1coeff; //CSDR at stoichiometic m.f.
     }
 }
 
diff --git a/updateYi.H b/updateYi.H
index 09a38c7..d27caf6 100644
--- a/updateYi.H
+++ b/updateYi.H
@@ -1,6 +1,11 @@
 //force the saving of the old-time values
 //important!!
 
+forAll(Y, yi)
+{
+	Y[yi].oldTime();
+}
+
 forAll(rho, cellI)
 {
     scalar& Nst = Neta[lowerN][cellI];