From 2e799dbc774ff55a7c197a86fba18bc8b565686f Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Tue, 23 Jan 2018 04:54:04 +0900
Subject: [PATCH] implemented EMD combustionModel (added edm() to
 chemistryModel)

---
 src/combustionModels/EDM/EDM.C                | 45 ++++++++-
 .../basicChemistryModel/basicChemistryModel.H |  8 ++
 .../chemistryModel/chemistryModel.C           | 98 +++++++++++++++++++
 .../chemistryModel/chemistryModel.H           |  8 ++
 4 files changed, 156 insertions(+), 3 deletions(-)

diff --git a/src/combustionModels/EDM/EDM.C b/src/combustionModels/EDM/EDM.C
index d672d2c06..6ec62221f 100644
--- a/src/combustionModels/EDM/EDM.C
+++ b/src/combustionModels/EDM/EDM.C
@@ -77,12 +77,51 @@ void Foam::combustionModels::EDM<Type>::correct()
 {
     if (this->active())
     {
+        tmp<volScalarField> tk(this->turbulence().k());
+        const volScalarField& k = tk();
+
+        tmp<volScalarField> tepsilon(this->turbulence().epsilon());
+        const volScalarField& epsilon = tepsilon();
+
+        this->chemistryPtr_->edm(epsilon/k, A_, B_);
+
         if (finiteRate_)
         {
+            // Copy edm rate stored in chemistryPtr_->RR(i)
+            PtrList<DimensionedField<scalar, volMesh>> edmRates
+            (
+                this->thermo().composition().Y().size()
+            );
+            forAll(edmRates, i)
+            {
+                edmRates.set
+                (
+                    i,
+                    new DimensionedField<scalar, volMesh>
+                    (
+                        this->chemistryPtr_->RR(i)
+                    )
+                );
+            }
+
             this->chemistryPtr_->calculate();
-        }
-        else
-        {
+
+            forAll(edmRates, i)
+            {
+                DimensionedField<scalar, volMesh> &finiteRatei
+                (
+                    this->chemistryPtr_->RR(i)
+                );
+
+                forAll(finiteRatei, celli)
+                {
+                    finiteRatei[celli] = min(finiteRatei[celli], edmRates[i][celli]);
+                    finiteRatei[celli] =
+                        finiteRatei[celli]*edmRates[i][celli]
+                        / (finiteRatei[celli] + edmRates[i][celli] + SMALL);
+                }
+
+            }
         }
     }
 }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index eb927e71b..ade959215 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -159,6 +159,14 @@ public:
                 //- Calculates the reaction rates
                 virtual void calculate() = 0;
 
+                //- Calculates the eddy dissipation reaction rates
+                virtual void edm
+                (
+                    const scalarField& mixing,
+                    const scalar A,
+                    const scalar B
+                ) = 0;
+
                 //- Solve the reaction system for the given time step
                 //  and return the characteristic time
                 virtual scalar solve(const scalar deltaT) = 0;
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 9747e979c..91f2143d3 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -756,6 +756,104 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
 }
 
 
+template<class CompType, class ThermoType>
+void Foam::chemistryModel<CompType, ThermoType>::edm
+(
+    const scalarField& mixing,
+    const scalar A,
+    const scalar B
+)
+{
+    if (!this->chemistry_)
+    {
+        return;
+    }
+
+    const volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            this->time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->thermo().rho()
+    );
+
+    forAll(rho, celli)
+    {
+        const scalar rhoi = rho[celli];
+        const scalar mixingi = mixing[celli];
+
+        scalarField Yi(nSpecie_, 0.0);
+        for (label k=0; k<nSpecie_; k++)
+        {
+            Yi[k] = Y_[k][celli];
+        }
+
+        scalarField om(nSpecie_, 0.0);
+
+        scalar minYR = 1.0;
+        scalar YP = 0.0;
+        forAll(reactions_, m)
+        {
+            const Reaction<ThermoType>& R = reactions_[m];
+
+            scalar num = 0.0;
+            scalar den = 0.0;
+
+            forAll(R.lhs(), s)
+            {
+                const label k = R.lhs()[s].index;
+                const scalar vk = R.lhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                minYR = min(minYR, Yi[k]/Wk/vk);
+            }
+
+            forAll(R.rhs(), s)
+            {
+                const label k = R.rhs()[s].index;
+                const scalar vk = R.rhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                num += Yi[k];
+                den += vk*Wk;
+            }
+            YP = num/den;
+
+            scalar omegam = min (A*rhoi*mixingi*minYR, A*B*rhoi*mixingi*YP);
+
+            forAll(R.lhs(), s)
+            {
+                const label k = R.lhs()[s].index;
+                const scalar vk = R.lhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                om[k] -= vk*Wk*omegam;
+            }
+
+            forAll(R.rhs(), s)
+            {
+                const label k = R.rhs()[s].index;
+                const scalar vk = R.rhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                om[k] += vk*Wk*omegam;
+            }
+        }
+
+        for (label k=0; k<nSpecie_; k++)
+        {
+            RR_[k][celli] = om[k];
+        }
+    }
+}
+
+
 template<class CompType, class ThermoType>
 template<class DeltaTType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 720ad96b9..d21295eb4 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -191,6 +191,14 @@ public:
         //- Calculates the reaction rates
         virtual void calculate();
 
+        //- Calculates the eddy dissipation reaction rates
+        virtual void edm
+        (
+            const scalarField& mixing,
+            const scalar A,
+            const scalar B
+        );
+
 
         // Chemistry model functions (overriding abstract functions in
         // basicChemistryModel.H)