slfm with energy, working.

.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

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;
+}

Make/files

@@ -14,4 +14,4 @@ FlameStructure/FlameStructureIO.C
 
 SLFMFoam.C
 
-EXE = $(FOAM_APPBIN)/SLFMFoam
+EXE = $(FOAM_APPBIN)/slfmEnergyFoam

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

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);

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())
+        runTime.write();
-        {
-            #include "updateYi.H" //update species
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

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"

pEqn.H

@@ -47,6 +47,7 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    rho = thermo.rho();
     rho = max(rho, rhoMin);
     rho = min(rho, rhoMax);
     rho.relax();

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.
     }
 }
 

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];