From 454d0c8cf0830c6ba243a5ef7f58cd82d7ed9853 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Thu, 17 Aug 2017 09:40:17 +0900
Subject: [PATCH] LagrangianCMCFoam solves absolute enthalpy equation

---
 .../constant/thermophysicalProperties         |  2 +-
 .../FlameD_LagrangianCMC/system/fvSchemes     |  2 +-
 .../FlameD_LagrangianCMC/system/fvSolution    |  4 +-
 .../LagrangianCMCFoam/CMCintegration.H        | 59 ++++++-------------
 solvers_post/LagrangianCMCFoam/EEqn.H         |  5 +-
 .../LagrangianCMCFoam/LagrangianCMCFoam.C     |  1 -
 solvers_post/LagrangianCMCFoam/Make/options   |  1 +
 solvers_post/LagrangianCMCFoam/createFields.H |  2 +-
 8 files changed, 26 insertions(+), 50 deletions(-)

diff --git a/examples/FlameD_LagrangianCMC/constant/thermophysicalProperties b/examples/FlameD_LagrangianCMC/constant/thermophysicalProperties
index 0a8e826..79245b9 100755
--- a/examples/FlameD_LagrangianCMC/constant/thermophysicalProperties
+++ b/examples/FlameD_LagrangianCMC/constant/thermophysicalProperties
@@ -21,7 +21,7 @@ thermoType
     mixture         reactingMixture;
     transport       sutherland;
     thermo          janaf;
-    energy          sensibleEnthalpy;
+    energy          absoluteEnthalpy;
     equationOfState perfectGas;
     specie          specie;
 }
diff --git a/examples/FlameD_LagrangianCMC/system/fvSchemes b/examples/FlameD_LagrangianCMC/system/fvSchemes
index 79ce9ac..8b73dcb 100755
--- a/examples/FlameD_LagrangianCMC/system/fvSchemes
+++ b/examples/FlameD_LagrangianCMC/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss linearUpwindV cellLimited Gauss linear 1;
     div(phi,Yi_h)   Gauss multivariateSelection
     {
-	h linearUpwind cellLimited Gauss linear 1;
+	ha linearUpwind cellLimited Gauss linear 1;
 	H2 linearUpwind cellLimited Gauss linear 1;
 	H linearUpwind cellLimited Gauss linear 1;
 	O linearUpwind cellLimited Gauss linear 1;
diff --git a/examples/FlameD_LagrangianCMC/system/fvSolution b/examples/FlameD_LagrangianCMC/system/fvSolution
index 11fc7d5..0f40fac 100755
--- a/examples/FlameD_LagrangianCMC/system/fvSolution
+++ b/examples/FlameD_LagrangianCMC/system/fvSolution
@@ -46,14 +46,14 @@ solvers
         relTol           0.0;
     }
 
-    "(U|Yi|h|k|epsilon|mf|mfVar|Fk)"
+    "(U|Yi|ha|k|epsilon|mf|mfVar|Fk)"
     {
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-08;
         relTol          0;
     }
-    "(U|Yi|h|k|epsilon|mf|mfVar|Fk)Final"
+    "(U|Yi|ha|k|epsilon|mf|mfVar|Fk)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/solvers_post/LagrangianCMCFoam/CMCintegration.H b/solvers_post/LagrangianCMCFoam/CMCintegration.H
index 636a991..cc12438 100644
--- a/solvers_post/LagrangianCMCFoam/CMCintegration.H
+++ b/solvers_post/LagrangianCMCFoam/CMCintegration.H
@@ -1,50 +1,27 @@
-                
+
         forAll(rho, celli)
         {
-                if(F_total[celli] > 0)
+            if(F_total[celli] > 0)
+            {
+                scalarField Ytemp(Ysize, 0.0);
+
+                for(label i=0 ; i<Ysize ; i++) //mass fraction integration over an eta-space for species i
                 {
-			scalarField Ytemp(Ysize, 0.0);
-			
-                        for(label i=0 ; i < Ysize ; i++) //reaction rate integration over an eta-space for species i
+                    for(label k=0 ; k<=group ; k++)
+                    {
+                        for(label j=0 ; j<=etamax ; j++)
                         {
-        			for(label k=0 ; k<=group ; k++)
-        			{
-					for(label j=0 ; j<= etamax ; j++)
-					{
-						pdf[j] = Peta[j][celli];
-						f[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
-					}
-					Ytemp[i] += Ffrac[k][celli]*integration(deltaftn[celli], MFcut, Neta, pdf, f);
-                		}
-        		}
-
-                        if( Ytemp[o2index] > Y[o2index][celli] || Ytemp[fuelindex] > Y[fuelindex][celli])
-                        {
-                                SUMSh[celli] = 0;
-                                //keep Yi value
-                        }
-                        else
-                        {
-                                scalarField ReactionRate(Ysize, 0.0);
-                                scalar tShtemp(0);
-
-                                for(label i = 0 ; i<Ysize ; i++)
-                                {
-                                        ReactionRate[i] = rho[celli] * (Ytemp[i] - Y[i][celli]) / runTime.deltaT().value();
-                                        tShtemp -= chemistry->calculateShCMC( ReactionRate , i);
-                                }
-                                SUMSh[celli] = tShtemp;
-
-                                for(label i = 0 ; i <Ysize ; i++)
-                                {
-                                        Y[i][celli] = Ytemp[i];
-                                }
-
+                            pdf[j] = Peta[j][celli];
+                            f[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
                         }
+                        Ytemp[i] += Ffrac[k][celli]*integration(deltaftn[celli], MFcut, Neta, pdf, f);
+                    }
                 }
-                else
+
+                for(label i=0 ; i<Ysize ; i++)
                 {
-                        SUMSh[celli] = 0;
+                    Y[i][celli] = Ytemp[i];
                 }
+            }
         }
- 
+
diff --git a/solvers_post/LagrangianCMCFoam/EEqn.H b/solvers_post/LagrangianCMCFoam/EEqn.H
index 750e6a9..b1d6ea2 100644
--- a/solvers_post/LagrangianCMCFoam/EEqn.H
+++ b/solvers_post/LagrangianCMCFoam/EEqn.H
@@ -6,7 +6,7 @@
         fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
-            he.name() == "e"
+            he.name() == "ea"
           ? fvc::div
             (
                 fvc::absolute(phi/fvc::interpolate(rho), U),
@@ -17,8 +17,7 @@
         )
       - fvm::laplacian(1.47*turbulence->mut(), he)
      ==
-        SUMSh
-       +fvOptions(rho, he)
+        fvOptions(rho, he)
     );
 
     EEqn.relax();
diff --git a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
index bc833f1..8ed9abd 100644
--- a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
+++ b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
@@ -96,7 +96,6 @@ int main(int argc, char *argv[])
                 init_start_CMC = false;
         }
         #include "CMCequation.H" //Solve CMC equations
-        #include "YEqn.H"
         #include "CMCintegration.H" //Conditional field integration and get reaction rate
         #include "EEqn.H"
 
diff --git a/solvers_post/LagrangianCMCFoam/Make/options b/solvers_post/LagrangianCMCFoam/Make/options
index 3e56e5f..8da350c 100644
--- a/solvers_post/LagrangianCMCFoam/Make/options
+++ b/solvers_post/LagrangianCMCFoam/Make/options
@@ -18,6 +18,7 @@ EXE_LIBS = \
     -L$(FOAM_USER_LIBBIN) \
         -lchemistryModel_POSTECH \
         -lcombustionModels_POSTECH \
+        -lthermoAdd \
     -lfiniteVolume \
     -lfvOptions \
     -lmeshTools \
diff --git a/solvers_post/LagrangianCMCFoam/createFields.H b/solvers_post/LagrangianCMCFoam/createFields.H
index f8201a2..6fb50cd 100644
--- a/solvers_post/LagrangianCMCFoam/createFields.H
+++ b/solvers_post/LagrangianCMCFoam/createFields.H
@@ -10,7 +10,7 @@ autoPtr<combustionModels::psiChemistryCombustion> reaction
 autoPtr<psiChemistryModel> chemistry = reaction->chem();
 psiReactionThermo& thermo = chemistry->thermo();
 
-thermo.validate(args.executable(), "h", "e");
+thermo.validate(args.executable(), "ha");
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();