From f752bb2a503af37aacc11d7ce1182ddb1f50de82 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Thu, 1 Feb 2018 23:10:09 +0900
Subject: [PATCH] replaced hard coded constants

---
 .../MossBrookesSoot/MossBrookesSoot.C         | 71 ++++++++-----------
 .../MossBrookesSoot/MossBrookesSoot.H         |  9 +--
 .../khanGreeveSoot/khanGreeveSoot.C           |  7 +-
 3 files changed, 36 insertions(+), 51 deletions(-)

diff --git a/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C b/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C
index 33abb5d9a..256a5bc64 100644
--- a/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C
+++ b/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C
@@ -26,6 +26,7 @@ License
 #include "MossBrookesSoot.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvm.H"
+#include "constants.H"
 #include "basicSpecieMixture.H"
 
 
@@ -269,6 +270,7 @@ Foam::radiation::MossBrookesSoot::MossBrookesSoot
         mesh_,
         dimensionedScalar("SootMeanDiameter", dimLength, 0.0)
     ),
+    Nnorm("Nnorm", dimless, 1e15),
     PrtSoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
     rhoSoot(readScalar(coeffsDict_.lookup("SootDensity"))),
     Calpha(readScalar(coeffsDict_.lookup("Calpha"))),
@@ -279,9 +281,7 @@ Foam::radiation::MossBrookesSoot::MossBrookesSoot
     etaColl(readScalar(coeffsDict_.lookup("etaColl"))),
     Talpha(readScalar(coeffsDict_.lookup("Talpha"))),
     Tgamma(readScalar(coeffsDict_.lookup("Tgamma"))),
-    Navog(readScalar(coeffsDict_.lookup("AvogadroNumber"))),
     Cbeta(readScalar(coeffsDict_.lookup("Cbeta"))),
-    Boltzmann(readScalar(coeffsDict_.lookup("BoltzmannConst"))),
     radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect"))),
     cliqueNuc(readScalar(coeffsDict_.lookup("NucExpCorrectConstant"))),
     cliqueMassGrow(readScalar(coeffsDict_.lookup("MassGrowCorrectConstant")))
@@ -381,34 +381,13 @@ void Foam::radiation::MossBrookesSoot::correct()
 
     NucConcEqn.solve();
 
-    scalar minusindicator1 = 0.0;
-    scalar minusindicator2 = 0.0;
 
-    forAll(rho_, celli)
-    {
-        if (soot_[celli] < 0.0) {
-            soot_[celli] = 0.0;
-            minusindicator1 = 1.0;
-        }
-        if (NucConc_[celli] < 0.0) {
-            NucConc_[celli] =0.0;
-            minusindicator2 = 1.0;
-        }
+    soot_.max(0.0);
+    NucConc_.max(0.0);
 
-        NumDen_[celli] = rho_[celli]*(NucConc_[celli]*1e+15);
+    NumDen_ = rho_*NucConc_*Nnorm;
 
-        sootVF_[celli] = soot_[celli]*rho_[celli]/rhoSoot;
-
-    }
-    if (minusindicator1 == 1.0 and minusindicator2 == 0.0) {
-        Info<< "Minus soot value! corrected to zero" <<endl;
-    }
-    else if (minusindicator1 == 0.0 and minusindicator2 ==1.0) {
-        Info<< "Minus NucConc value! corrected to zero" <<endl;
-    }
-    else if (minusindicator1 == 1.0 and minusindicator2 == 1.0) {
-        Info<< "Minus soot and NucConc value! corrected to zero" << endl;
-    }
+    sootVF_ = soot_*rho_/dimensionedScalar("rhoSoot", dimDensity, rhoSoot);
 }
 
 
@@ -430,6 +409,12 @@ void Foam::radiation::MossBrookesSoot::calcSource()
     const volScalarField XSgs_ (mixture_.W()*Ygrow_/mixture_.W(growI));
     const volScalarField XOxid_(mixture_.W()*Yoxid_/mixture_.W(oxidI));
 
+    const dimensionedScalar NA (Foam::constant::physicoChemical::NA);
+    const dimensionedScalar kB (Foam::constant::physicoChemical::k);
+    const dimensionedScalar R (Foam::constant::physicoChemical::R);
+
+    const scalar pi314 = Foam::constant::mathematical::pi;
+
     forAll(rho_, celli)
     {
         const scalar XPreci = XPrec_[celli];
@@ -441,36 +426,38 @@ void Foam::radiation::MossBrookesSoot::calcSource()
         const scalar sooti = soot_[celli];
         const scalar NumDeni = NumDen_[celli];
 
-        const scalar MolarConcPreci = (XPreci*pi)/(8314.5*Ti);
-        const scalar MolarConcSgsi = (XSgsi*pi)/(8314.5*Ti);
-        const scalar MolarConcOxidi = (XOxidi*pi)/(8314.5*Ti);
+        // concentrations in [kmol/m3]
+        const scalar MolarConcPreci = (XPreci*pi)/(1e3*R.value()*Ti);
+        const scalar MolarConcSgsi  = (XSgsi *pi)/(1e3*R.value()*Ti);
+        const scalar MolarConcOxidi = (XOxidi*pi)/(1e3*R.value()*Ti);
 
-        scalar SootMeanDiameteri = 0.0;
-
-        if (NumDeni == 0.0) {
+        if (NumDeni == 0.0)
+        {
             sootMeanDiameter_[celli] = 0.0;
             Sf_[celli] = 0.0;
         }
-        else {
-            const scalar parameter = (6.0/3.1416)*(rhoi/rhoSoot)*(sooti/NumDeni);
-            sootMeanDiameter_[celli] = pow(parameter,0.33333);
-            SootMeanDiameteri = sootMeanDiameter_[celli];
-            Sf_[celli] = (3.1416)*pow(SootMeanDiameteri,2.0)*NumDeni;
+        else
+        {
+            const scalar parameter = (6.0/pi314)*(rhoi/rhoSoot)*(sooti/NumDeni);
+            sootMeanDiameter_[celli] = pow(parameter, 1./3.);
+            Sf_[celli] = pi314*pow(sootMeanDiameter_[celli],2.0)*NumDeni;
         }
 
+        const scalar sootMeanDiameteri = sootMeanDiameter_[celli];
+
         Snuc_[celli] = Mp*Calpha*exp(-cliqueNuc*Talpha/Ti)*MolarConcPreci;
 
         Smassgrow_[celli] = Cgamma*MolarConcSgsi*exp(-cliqueMassGrow*Tgamma/Ti)*Sf_[celli];
 
-        Soxid_[celli] = Coxid*Cw*etaColl*MolarConcOxidi*pow(Ti,0.5)*Sf_[celli];
+        Soxid_[celli] = Coxid*Cw*etaColl*MolarConcOxidi*sqrt(Ti)*Sf_[celli];
 
         Snet_[celli] = Snuc_[celli] + Smassgrow_[celli] - Soxid_[celli];
 
-        NSnuc_[celli] = (Calpha*Navog*MolarConcPreci*exp(-cliqueNuc*Talpha/Ti))*(1e-15);
+        NSnuc_[celli] = (Calpha*NA.value()*MolarConcPreci*exp(-cliqueNuc*Talpha/Ti))/Nnorm.value();
 
-        aggParam_[celli] = 24.0*Boltzmann*Ti/rhoSoot;
+        aggParam_[celli] = 24.0*kB.value()*Ti/rhoSoot;
 
-        NSagg_[celli] = (Cbeta*pow((aggParam_[celli]),0.5)*pow(SootMeanDiameteri,0.5)*pow(NumDeni,2.0))*(1e-15);
+        NSagg_[celli] = (Cbeta*sqrt(aggParam_[celli])*sqrt(sootMeanDiameteri)*pow(NumDeni,2.0))/Nnorm.value();
 
         NSnet_[celli] = NSnuc_[celli] - NSagg_[celli];
     }
diff --git a/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H b/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H
index b3589f2a6..35dfb3900 100644
--- a/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H
+++ b/src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H
@@ -159,6 +159,9 @@ class MossBrookesSoot
         //- Density of soot particle
         scalar rhoSoot;
 
+        //- soot number density normalization factor
+        const dimensionedScalar Nnorm;
+
         //- Model constant for soot inception rate
         scalar Calpha;
 
@@ -183,15 +186,9 @@ class MossBrookesSoot
         //- Activation temperature for soot growth rate
         scalar Tgamma;
 
-        //- Avogadro's number
-        scalar Navog;
-
         //- Model constant for coagulation rate
         scalar Cbeta;
 
-        //- Boltzmann constant
-        scalar Boltzmann;
-
         //- Soot radiative effect on/off flag
         bool radiationActive_;
 
diff --git a/src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C b/src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C
index ccb7e45da..e6f5e1185 100644
--- a/src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C
+++ b/src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C
@@ -348,12 +348,13 @@ void Foam::radiation::khanGreeveSoot::calcSource()
 
         S1_[celli] = A_*rhoi*sooti*epsiloni/ki;
 
-        if (sooti == 0.0 and Yfueli == 0.0 ){
+        if (sooti == 0.0 and Yfueli == 0.0 )
+        {
             S2_[celli] = 0.0;
         }
-        else {
+        else 
+        {
             S2_[celli] = A_*rhoi*(Yoxi/nuSoot_)*(sooti*nuSoot_/(sooti*nuSoot_+(Yfueli)*nuFuel_))*epsiloni/ki;
-
         }
 
         Scomb_[celli] = min(S1_[celli], S2_[celli]);