From 5bd56ebdaf8271c92d442bfc79b4bbd94fe130e4 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 22 Feb 2011 11:55:00 +0000
Subject: [PATCH] ENH: grey absorption with multispecies solvers

---
 .../global/unitConversion/unitConversion.H    |  12 ++
 .../radiationModels/Make/options              |  17 ++-
 .../radiationModels/radiationModel/P1/P1.C    |   2 +-
 .../greyMeanAbsorptionEmission.C              | 117 +++++++++++++-----
 .../greyMeanAbsorptionEmission.H              |  17 +--
 .../constant/radiationProperties              |  56 +++++++++
 .../constant/radiationProperties              |  56 +++++++++
 7 files changed, 236 insertions(+), 41 deletions(-)

diff --git a/src/OpenFOAM/global/unitConversion/unitConversion.H b/src/OpenFOAM/global/unitConversion/unitConversion.H
index 3255a4d805..9913ed3647 100644
--- a/src/OpenFOAM/global/unitConversion/unitConversion.H
+++ b/src/OpenFOAM/global/unitConversion/unitConversion.H
@@ -53,6 +53,18 @@ inline scalar radToDeg(const scalar rad)
     return (rad * 180.0/constant::mathematical::pi);
 }
 
+//- Conversion from atm to Pa
+inline scalar atmToPa(const scalar atm)
+{
+    return (atm * 101035.0);
+}
+
+//- Conversion from atm to Pa
+inline scalar paToAtm(const scalar pa)
+{
+    return (pa / 101035.0);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/radiationModels/Make/options b/src/thermophysicalModels/radiationModels/Make/options
index 91f332c71f..f9a6356a6d 100644
--- a/src/thermophysicalModels/radiationModels/Make/options
+++ b/src/thermophysicalModels/radiationModels/Make/options
@@ -3,12 +3,25 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude
+
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lbasicThermophysicalModels \
     -lspecie \
     -lbasicSolidThermo \
-    -lmeshTools
+    -lmeshTools \
+    -lSLGThermo \
+    -lsolidMixtureProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lliquidProperties
+
 
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C b/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
index f59b5632f4..907ffc4566 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
+++ b/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
@@ -89,7 +89,7 @@ Foam::radiation::P1::P1(const volScalarField& T)
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 62c491b104..21dfa9786d 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -25,6 +25,8 @@ License
 
 #include "greyMeanAbsorptionEmission.H"
 #include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,13 +58,15 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
     coeffsDict_((dict.subDict(typeName + "Coeffs"))),
     speciesNames_(0),
     specieIndex_(0),
-    lookUpTable_
+    lookUpTablePtr_(),
+    thermo_
     (
-        fileName(coeffsDict_.lookup("lookUpTableFileName")),
-        mesh.time().constant(),
-        mesh
+        mesh,
+        const_cast<basicThermo&>
+        (
+            mesh.lookupObject<basicThermo>("thermophysicalProperties")
+        )
     ),
-    thermo_(mesh.lookupObject<basicThermo>("thermophysicalProperties")),
     EhrrCoeff_(readScalar(coeffsDict_.lookup("EhrrCoeff"))),
     Yj_(nSpecies_)
 {
@@ -83,17 +87,45 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
         nFunc++;
     }
 
+    if (coeffsDict_.found("lookUpTableFileName"))
+    {
+        const word name = coeffsDict_.lookup("lookUpTableFileName");
+        if (name != "none")
+        {
+            lookUpTablePtr_.set
+            (
+                new interpolationLookUpTable<scalar>
+                (
+                    fileName(coeffsDict_.lookup("lookUpTableFileName")),
+                    mesh.time().constant(),
+                    mesh
+                )
+            );
+
+            if (!mesh.foundObject<volScalarField>("ft"))
+            {
+                FatalErrorIn
+                (
+                    "Foam::radiation::greyMeanAbsorptionEmission(const"
+                    "dictionary& dict, const fvMesh& mesh)"
+                )   << "specie ft is not present to use with "
+                    << "lookUpTableFileName " << nl
+                    << exit(FatalError);
+            }
+        }
+    }
+
     // Check that all the species on the dictionary are present in the
     // look-up table and save the corresponding indices of the look-up table
 
     label j = 0;
     forAllConstIter(HashTable<label>, speciesNames_, iter)
     {
-        if (mesh.foundObject<volScalarField>("ft"))
+        if (!lookUpTablePtr_.empty())
         {
-            if (lookUpTable_.found(iter.key()))
+            if (lookUpTablePtr_().found(iter.key()))
             {
-                label index = lookUpTable_.findFieldIndex(iter.key());
+                label index = lookUpTablePtr_().findFieldIndex(iter.key());
 
                 Info<< "specie: " << iter.key() << " found on look-up table "
                     << " with index: " << index << endl;
@@ -120,18 +152,32 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
                     "dictionary& dict, const fvMesh& mesh)"
                 )   << "specie: " << iter.key()
                     << " is neither in look-up table: "
-                    << lookUpTable_.tableName()
+                    << lookUpTablePtr_().tableName()
                     << " nor is being solved" << nl
                     << exit(FatalError);
             }
         }
+        else if (mesh.foundObject<volScalarField>(iter.key()))
+        {
+            volScalarField& Y =
+                const_cast<volScalarField&>
+                (
+                    mesh.lookupObject<volScalarField>(iter.key())
+                );
+
+            Yj_.set(j, &Y);
+            specieIndex_[iter()] = 0;
+            j++;
+        }
         else
         {
             FatalErrorIn
             (
                 "Foam::radiation::greyMeanAbsorptionEmission(const"
                 "dictionary& dict, const fvMesh& mesh)"
-            )   << "specie ft is not present " << nl
+            )   << " there is not lookup table and the specie" << nl
+                << iter.key() << nl
+                << " is not found " << nl
                 << exit(FatalError);
 
         }
@@ -149,11 +195,10 @@ Foam::radiation::greyMeanAbsorptionEmission::~greyMeanAbsorptionEmission()
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
 {
-    const volScalarField& T = thermo_.T();
-    const volScalarField& p = thermo_.p();
-    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
+    const volScalarField& T = thermo_.thermo().T();
+    const volScalarField& p = thermo_.thermo().p();
 
-    label nSpecies = speciesNames_.size();
+    const basicMultiComponentMixture& mixture = thermo_.carrier();
 
     tmp<volScalarField> ta
     (
@@ -168,48 +213,60 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
                 IOobject::NO_WRITE
             ),
             mesh(),
-            dimensionedScalar("a", dimless/dimLength, 0.0)
+            dimensionedScalar("a", dimless/dimLength, 0.0),
+            zeroGradientFvPatchVectorField::typeName
         )
     );
 
     scalarField& a = ta().internalField();
 
-    forAll(a, i)
+    forAll(a, cellI)
     {
-        const List<scalar>& species = lookUpTable_.lookUp(ft[i]);
-
-        for (label n=0; n<nSpecies; n++)
+        forAllConstIter(HashTable<label>, speciesNames_, iter)
         {
-            label l = 0;
-            scalar Yipi = 0;
+            label n = iter();
+            scalar Xipi = 0.0;
             if (specieIndex_[n] != 0)
             {
+                //Specie found in the lookUpTable.
+                const volScalarField& ft =
+                    mesh_.lookupObject<volScalarField>("ft");
+
+                const List<scalar>& Ynft = lookUpTablePtr_().lookUp(ft[cellI]);
                 //moles x pressure [atm]
-                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+                Xipi = Ynft[specieIndex_[n]]*paToAtm(p[cellI]);
             }
             else
             {
-                // mass fraction
-                Yipi = Yj_[l][i];
-                l++;
+                scalar invWt = 0.0;
+                forAll (mixture.Y(), s)
+                {
+                    invWt += mixture.Y(s)[cellI]/mixture.W(s);
+                }
+
+                label index = mixture.species()[iter.key()];
+                scalar Xk = mixture.Y(index)[cellI]/(mixture.W(index)*invWt);
+
+                Xipi = Xk*paToAtm(p[cellI]);
             }
 
-            const absorptionCoeffs::coeffArray& b = coeffs_[n].coeffs(T[i]);
+            const absorptionCoeffs::coeffArray& b = coeffs_[n].coeffs(T[cellI]);
 
-            scalar Ti = T[i];
+            scalar Ti = T[cellI];
             // negative temperature exponents
             if (coeffs_[n].invTemp())
             {
-                Ti = 1./T[i];
+                Ti = 1.0/T[cellI];
             }
-            a[i] +=
-                Yipi
+            a[cellI] +=
+                Xipi
                *(
                     ((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti
                   + b[0]
                 );
         }
     }
+    ta().correctBoundaryConditions();
     return ta;
 }
 
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
index 6e8fbe8c45..a41d2adf09 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
@@ -89,7 +89,8 @@ SourceFiles
 #include "absorptionEmissionModel.H"
 #include "HashTable.H"
 #include "absorptionCoeffs.H"
-#include "basicThermo.H"
+#include "SLGThermo.H"
+//#include "basicThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -110,10 +111,10 @@ public:
 
     // Public data
 
-        // - Maximum number of species considered for absorptivity
+        // Maximum number of species considered for absorptivity
         static const int nSpecies_ = 5;
 
-        //  Absorption Coefficients
+        // Absorption Coefficients
         absorptionCoeffs coeffs_[nSpecies_];
 
 
@@ -127,14 +128,14 @@ private:
         //- Hash table of species names
         HashTable<label> speciesNames_;
 
-        // Indices of species in the look-up table
+        //- Indices of species in the look-up table
         FixedList<label, nSpecies_> specieIndex_;
 
-        // Look-up table of species related to ft
-        mutable interpolationLookUpTable<scalar> lookUpTable_;
+        //- Look-up table of species related to ft
+        mutable autoPtr<interpolationLookUpTable<scalar> > lookUpTablePtr_;
 
-        // Thermo package
-        const basicThermo& thermo_;
+        //- SLG thermo package
+        SLGThermo thermo_;
 
         //- Emission constant coefficient
         const scalar EhrrCoeff_;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
index 37b3eefe66..c030050fc4 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
@@ -135,6 +135,62 @@ greyMeanAbsorptionEmissionCoeffs
         );
     }
 
+    N2
+    {
+        Tcommon         300;
+        invTemp         false;
+        Tlow            200; 
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0.01
+            0   
+            0   
+            0   
+            0   
+            0   
+        );
+        hiTcoeffs
+        (
+            0.01
+            0
+            0
+            0
+            0
+            0
+        );
+
+    }
+
+    O2
+    {
+        Tcommon         300;
+        invTemp         false;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0.01 
+            0       
+            0         
+            0         
+            0         
+            0         
+        );
+        hiTcoeffs     
+        (
+            0.01
+            0
+            0
+            0
+            0
+            0
+        );
+
+    }
+
 }
 
 scatterModel    constantScatter;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
index 500834ae4a..aaee75d4a1 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
@@ -135,6 +135,62 @@ greyMeanAbsorptionEmissionCoeffs
             0
         );
     }
+    
+    N2
+    {
+        Tcommon         300;
+        invTemp         false;
+        Tlow            200; 
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0.01
+            0   
+            0   
+            0   
+            0   
+            0   
+        );
+        hiTcoeffs
+        (
+            0.01
+            0
+            0
+            0
+            0
+            0
+        );
+
+    }
+
+    O2
+    {
+        Tcommon         300;
+        invTemp         false;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0.01 
+            0       
+            0         
+            0         
+            0         
+            0         
+        );
+        hiTcoeffs     
+        (
+            0.01
+            0
+            0
+            0
+            0
+            0
+        );
+
+    }
 
 }