View Issue Details

IDProjectCategoryView StatusLast Update
0002881OpenFOAMBugpublic2018-04-26 05:43
Reporterkevnor Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
Fixed in Versiondev 
Summary0002881: wideBandAbsorptionEmission coefficient bug and inconsistency with
DescriptionwideBandAbsorptionEmission::aCont attempts to use the coefficients list with inverted ordering (i.e. coeffs_[species][band] instead of coeffs_[band][species]) which results in incorrect coefficients being used and segmentation faults. I'm guessing that no-one else is currently using this model?!

Additionally, the behaviour of wideBandAbsorptionEmission is inconsistent with greyMeanAbsorptionEmission in reading the look-up tables and using mole-fractions versus mass-fractions for solved-for species.

I've attached a patch (for 5.x) that I believe corrects the bug and makes the behaviour consistent.
Steps To ReproduceRun any case with wideBandAbsorpionEmission model. E.g. modify fireFoam oppositeBurningPanels tutorial case to use this model.
TagsNo tags attached.

Activities

kevnor

2018-03-15 08:23

reporter  

wideBandAbsorptionEmission.patch (9,385 bytes)   
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index aa0f3829a..4608a67fb 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
@@ -25,6 +25,8 @@ License
 
 #include "wideBandAbsorptionEmission.H"
 #include "addToRunTimeSelectionTable.H"
+#include "basicSpecieMixture.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,12 +58,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
     coeffsDict_((dict.optionalSubDict(typeName + "Coeffs"))),
     speciesNames_(0),
     specieIndex_(label(0)),
-    lookUpTable_
-    (
-        fileName(coeffsDict_.lookup("lookUpTableFileName")),
-        mesh.time().constant(),
-        mesh
-    ),
+    lookUpTablePtr_(),
     thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
     Yj_(nSpecies_),
     totalWaveLength_(0)
@@ -95,7 +92,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
                 if (!speciesNames_.found(key))
                 {
                     FatalErrorInFunction
-                        << "specie: " << key << "is not in all the bands"
+                        << "specie: " << key << " is not in all the bands"
                         << nl << exit(FatalError);
                 }
             }
@@ -106,36 +103,82 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
     }
     nBands_ = nBand;
 
+    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"))
+            {
+                FatalErrorInFunction
+                    << "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
+    // look-up table and save the corresponding indices of the look-up table
 
     label j = 0;
     forAllConstIter(HashTable<label>, speciesNames_, iter)
     {
-        if (lookUpTable_.found(iter.key()))
+        if (!lookUpTablePtr_.empty())
         {
-            label index = lookUpTable_.findFieldIndex(iter.key());
-            Info<< "specie: " << iter.key() << " found in look-up table "
-                << " with index: " << index << endl;
-            specieIndex_[iter()] = index;
+            if (lookUpTablePtr_().found(iter.key()))
+            {
+                label index = lookUpTablePtr_().findFieldIndex(iter.key());
+
+                Info<< "specie: " << iter.key() << " found on look-up table "
+                    << " with index: " << index << endl;
+
+                specieIndex_[iter()] = index;
+            }
+            else if (mesh.foundObject<volScalarField>(iter.key()))
+            {
+                Yj_.set(j, &mesh.lookupObjectRef<volScalarField>(iter.key()));
+                specieIndex_[iter()] = 0;
+                j++;
+                Info<< "specie: " << iter.key() << " is being solved" << endl;
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "specie: " << iter.key()
+                    << " is neither in look-up table: "
+                    << lookUpTablePtr_().tableName()
+                    << " nor is being solved" << nl
+                    << exit(FatalError);
+            }
         }
         else if (mesh.foundObject<volScalarField>(iter.key()))
         {
             Yj_.set(j, &mesh.lookupObjectRef<volScalarField>(iter.key()));
-
-            specieIndex_[iter()] = 0.0;
+            specieIndex_[iter()] = 0;
             j++;
-            Info<< "species: " << iter.key() << " is being solved" << endl;
         }
         else
         {
             FatalErrorInFunction
-                << "specie: " << iter.key()
-                << " is neither in look-up table : "
-                << lookUpTable_.tableName() << " nor is being solved"
+                << " there is no lookup table and the specie" << nl
+                << iter.key() << nl
+                << " is not found " << nl
                 << exit(FatalError);
+
         }
     }
+    
 }
 
 
@@ -151,11 +194,11 @@ Foam::radiation::wideBandAbsorptionEmission::~wideBandAbsorptionEmission()
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
 {
+    const basicSpecieMixture& mixture =
+        dynamic_cast<const basicSpecieMixture&>(thermo_);
+
     const volScalarField& T = thermo_.T();
     const volScalarField& p = thermo_.p();
-    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
-
-    label nSpecies = speciesNames_.size();
 
     tmp<volScalarField> ta
     (
@@ -176,38 +219,47 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
 
     scalarField& a = ta.ref().primitiveFieldRef();
 
-    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.0;
+            label n = iter();
+            scalar Xipi = 0.0;
             if (specieIndex_[n] != 0)
             {
-                // moles x pressure [atm]
-                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+	        const volScalarField& ft =
+		    mesh_.lookupObject<volScalarField>("ft");
+
+                const List<scalar>& Ynft = lookUpTablePtr_().lookUp(ft[celli]);
+                //moles x pressure [atm]
+                Xipi = Ynft[specieIndex_[n]]*paToAtm(p[celli]);
             }
             else
             {
-                // mass fraction from species being solved
-                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]);
             }
 
-            scalar Ti = T[i];
+            scalar Ti = T[celli];
 
             const absorptionCoeffs::coeffArray& b =
-                coeffs_[n][bandI].coeffs(T[i]);
+	        coeffs_[bandI][n].coeffs(T[celli]);
 
-            if (coeffs_[n][bandI].invTemp())
+            if (coeffs_[bandI][n].invTemp())
             {
-                Ti = 1.0/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]
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
index e1fb10cf0..880537b68 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
@@ -37,11 +37,8 @@ Description
 
     The emission constant proportionality is specified per band (EhrrCoeff).
 
-    The coefficients for the species in the lookup table have to be specified
-    for use in moles x P [atm].i.e. (k[i] = species[i]*p*9.869231e-6).
-
-    The coefficients for CO and soot or any other added are multiplied by the
-    respective mass fraction being solved.
+    The coefficients for the species have to be specified for use in 
+    moles x P [atm], i.e. (k[i] = species[i]*p*9.869231e-6).
 
     The look Up table file should be in the constant directory.
 
@@ -156,8 +153,8 @@ private:
         //- Proportion of the heat released rate emitted
         FixedList<scalar, maxBands_> iEhrrCoeffs_;
 
-        //- Lookup 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 fluidThermo& thermo_;

henry

2018-03-15 16:24

manager   ~0009427

Do you have a test case which demonstrates the need for this change?

kevnor

2018-03-16 21:20

reporter   ~0009429

Using the attached radiationProperties as a replacement within the oppositeBurningPanels tutorial case provides a demonstration of the issues. It should simply replicate the results from using the greyMeanAbsorptionEmission model since it simply contains two bands with identical properties.

Attempting to run this case will exit with an error message that a string is expected under the entry lookupTableFileName, i.e. that 'none' is not accepted, unlike the greyMeanAbsorptionEmission model.

Specifying the attached "wideBandTable" dummy lookup table file allows the solver to run further, but it then halts when looking for the mixture fraction "ft" field, which isn't defined for the present combustion model (and indeed shouldn't be required since all of the required fields are being explicitly solved for).

Hacking the solver to manually add a dummy "ft" field again allows the solver to continue. However, it then generates a slew of warning messages of the form:

    usinf absCoeffs out of temperature range:
    0.0532686334933 -> 0.0326473267652; T = 298.15

Closer look at the code shows that this is because it is trying to use unitialised coefficients due to the described interchanged indices in accessing coeffs_ within wideBandAbsorptionEmission::aCont()

Correcting these interchanged indices allows the solver to run to completion. However, comparison of the results with those obtained using the greyMeanAbsorptionEmissionModel in the same radiationProperties shows differences. This is due to the wideBand model using the mass fraction of the species as the pre-multiplier for the absorption coefficients, while the greyMean (more sensibly) uses the partial pressure.
radiationProperties (4,959 bytes)   
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      radiationProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

radiation       on;

radiationModel  fvDOM;

fvDOMCoeffs
{
    nPhi    3;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
    nTheta  6;          // polar angles in PI (from Z to X-Y plane)
    tolerance   0.05;   // convergence tolerance for radiation iteration
    maxIter 3;          // maximum number of iterations
}

// Number of flow iterations per radiation iteration
solverFreq 10;

absorptionEmissionModel wideBandAbsorptionEmission;
//absorptionEmissionModel greyMeanAbsorptionEmission;

wideBandAbsorptionEmissionCoeffs
{
  lookUpTableFileName none;
  // lookUpTableFileName "wideBandTable";
  
  band0
  {
     bandLimits (1e-20 1e-6);
     EhrrCoeff       0.2;
     species
     {
         O2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 C3H8
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 H2O
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 CO2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 N2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }	 
     }
  }

  band1
  {
     bandLimits (1e-6 1e10);
     EhrrCoeff       0.2;
     species
     {
         O2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 C3H8
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 H2O
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 CO2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 N2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }	 
     }
  }
  
}

constantAbsorptionEmissionCoeffs
{
    absorptivity    absorptivity    [ m^-1 ]       0.01;
    emissivity      emissivity      [ m^-1 ]       0.01;
    E               E               [ kg m^-1 s^-3 ]  0;
}

greyMeanAbsorptionEmissionCoeffs
{
    lookUpTableFileName     none;

    EhrrCoeff               0.2;

         O2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 C3H8
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 H2O
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 CO2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }
	 N2
	 {
	     Tcommon         200.;
	     Tlow            200.;
	     Thigh           2500.;
	     invTemp         false;
	     loTcoeffs (0 0 0 0 0 0) ;
	     hiTcoeffs (.1 0 0 0 0 0);
	 }	 

}

scatterModel    none;

sootModel       none;

// ************************************************************************* //
radiationProperties (4,959 bytes)   

kevnor

2018-03-16 21:21

reporter  

wideBandTable (888 bytes)   
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      radiationProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

fields (
 {
    N 1;
    min 0;
    max 1;
    name "dummy";
 }
);

output ();
values (
	(1)
);
wideBandTable (888 bytes)   

henry

2018-03-16 21:55

manager   ~0009430

Does your patch resolve all of the above issues or will some of the changes listed need to be made to reproduce your findings? If so should any further changes be committed to OpenFOAM-dev?

kevnor

2018-03-17 18:17

reporter   ~0009431

I believe the patch resolves all of the above issues. Much of the code has just been transferred from the greyMeanAbsorptionEmission class.

henry

2018-03-21 18:11

manager   ~0009436

Thanks for the patches and test-case. I have updated the code for consistency with the OpenFOAM style: https://openfoam.org/dev/coding-style-guide/

and committed it to OpenFOAM-dev: c677ba11856b116677253ff0575e2c0e17973233

Could you check these updates?

For future patches it would be useful if you could ensure they conform to the OpenFOAM style, in particular avoid tabs as they are difficult to remove reliably without upsetting the indentation requiring time consuming editing by hand.

kevnor

2018-03-27 19:42

reporter   ~0009452

The commit seems to do the trick for OpenFOAM-dev, thanks. Sorry I forgot to "untabify" the updates before creating the patches!

Issue History

Date Modified Username Field Change
2018-03-15 08:23 kevnor New Issue
2018-03-15 08:23 kevnor File Added: wideBandAbsorptionEmission.patch
2018-03-15 16:24 henry Note Added: 0009427
2018-03-16 21:20 kevnor File Added: radiationProperties
2018-03-16 21:20 kevnor Note Added: 0009429
2018-03-16 21:21 kevnor File Added: wideBandTable
2018-03-16 21:55 henry Note Added: 0009430
2018-03-17 18:17 kevnor Note Added: 0009431
2018-03-21 18:11 henry Note Added: 0009436
2018-03-27 19:42 kevnor Note Added: 0009452
2018-03-27 21:11 henry Assigned To => henry
2018-03-27 21:11 henry Status new => resolved
2018-03-27 21:11 henry Resolution open => fixed
2018-03-27 21:11 henry Fixed in Version => dev