View Issue Details

IDProjectCategoryView StatusLast Update
0002882OpenFOAMPatchpublic2018-03-27 21:11
ReporterkevnorAssigned Tohenry 
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Product Version5.x 
Fixed in Versiondev 
Summary0002882: Radiative energy source terms broken for fvDOM radiation model with wideBandAbsorptionEmission
DescriptionThe energy transfer terms for the continuous phase (fvDOM::Ru and fvDOM::Rp) appear to be completely incorrect for the fvDOM radiation model coupled with the wideBand absorption-emission model! The continuous phase absorption coefficient for the first band is currently applied for the entire incident radiation (and all the other absorption coefficients are hence ignored).

I believe the expression should rather be a sum (over all bands) of the absorption coefficient for the band multiplied by the incident radiation within that band, i.e.

  Ru = Sum_i (a_i G_i + E_i)

where G_i = Sum_n I_(n,i) omega_n is the incident radiation within the band

Similarly, the Rp term should be based on the sum of the absorption coefficient
multiplied by the proportion of the black body energy flux in the corresponding band, i.e.

  Rp = 4 \sigma Sum_i (a_i (f(\lambda_{i+1} T)-f(\lambda_i T)))

where f is the blackBodyEmission::fLambdaT function

I've compiled my suggested changes into the attached patches. Note that the computation of the proportion of the black-body energy flux within each band required the addition of an extra function in the blackBodyEmission class.

I've currently computed the continuous phase absorption coefficients in the same way as in radiativeIntensityRay, i.e. by subtracting the dispersed phase absorption coefficient from the (cached) total absorption. However, since the continuous phase absorption now gets used three times (as compared to once for the total absorption), it would be more efficient to cache the continuous phase absorption instead. I can update the patch to use this approach if you prefer?

(See also https://bugs.openfoam.org/view.php?id=2881 )
TagsNo tags attached.

Activities

kevnor

2018-03-15 08:39

reporter  

wideBandfvDOM.patch (6,093 bytes)
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
index 1b50a08cf..81547766e 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
@@ -40,7 +40,7 @@ void Foam::radiation::absorptionCoeffs::checkT(const scalar T) const
     if (T < Tlow_ || T > Thigh_)
     {
         WarningInFunction
-            << "usinf absCoeff out of temperature range:" << nl
+            << "using absorptionCoeffs out of temperature range:" << nl
             << "    " << Tlow_ << " -> " << Thigh_ << ";  T = " << T
             << nl << endl;
     }
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
index b2fdb85d1..30b3467a1 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
@@ -197,6 +197,51 @@ Foam::scalar Foam::radiation::blackBodyEmission::fLambdaT
 }
 
 
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::blackBodyEmission::DeltaLambdaT
+(
+    const volScalarField& T,
+    const Vector2D<scalar>& band
+) const
+{
+    tmp<volScalarField> DeltaLambdaT
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Eb",
+                T.mesh().time().timeName(),
+                T.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+	    T.mesh(),
+	    dimensionedScalar("", dimless, 1.0)
+        )
+    );
+
+
+    if (band == Vector2D<scalar>::one)
+    {
+        return DeltaLambdaT;
+    }
+    else
+    {
+        scalarField& DeltaLambdaTf = DeltaLambdaT.ref();
+
+        forAll(T, i)
+        {
+            scalar T1 = fLambdaT(band[1]*T[i]);
+            scalar T2 = fLambdaT(band[0]*T[i]);
+            DeltaLambdaTf[i] *= T1 - T2;
+        }
+
+        return DeltaLambdaT;
+    }
+}
+
+
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::blackBodyEmission::EbDeltaLambdaT
 (
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
index d7f175fc5..ae9cee868 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
@@ -130,6 +130,12 @@ public:
                 const Vector2D<scalar>& band
             ) const;
 
+            //- Proportion of total energy at T from lambda1 to lambda2
+            tmp<Foam::volScalarField> DeltaLambdaT
+            (
+                const volScalarField& T,
+                const Vector2D<scalar>& band
+            ) const;
 
     // Edit
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
index 94d365081..b20d70781 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
@@ -455,7 +455,8 @@ void Foam::radiation::fvDOM::calculate()
 
 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
 {
-    return tmp<volScalarField>
+    // Construct using contribution from first frequency band
+    tmp<volScalarField> tRp
     (
         new volScalarField
         (
@@ -468,28 +469,73 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
                 IOobject::NO_WRITE,
                 false
             ),
-            // Only include continuous phase emission
-            4*absorptionEmission_->aCont()*physicoChemical::sigma
+	    (
+	        4*physicoChemical::sigma
+	      * (aLambda_[0]-absorptionEmission_->aDisp(0)())
+	      * blackBody_.DeltaLambdaT(T_,absorptionEmission_->bands(0))
+	    )
         )
     );
+
+    volScalarField& Rp=tRp.ref();
+
+    // Add contributions over remaining frequency bands
+    for (label j=1; j < nLambda_; j++)
+    {
+        Rp +=
+	(
+	    4*physicoChemical::sigma
+	  * (aLambda_[j]-absorptionEmission_->aDisp(j)())
+	  * blackBody_.DeltaLambdaT(T_,absorptionEmission_->bands(j))
+	);
+    }
+    
+    return tRp;
 }
 
 
 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
 Foam::radiation::fvDOM::Ru() const
 {
+    tmp<DimensionedField<scalar, volMesh>> tRu
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "Ru",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+	    mesh_,
+	    dimensionedScalar("zero",dimensionSet(1,-1,-3,0,0), 0.0)
+        )
+    );
 
-    const volScalarField::Internal& G =
-        G_();
+    DimensionedField<scalar, volMesh>& Ru=tRu.ref();
 
-    const volScalarField::Internal E =
-        absorptionEmission_->ECont()()();
+    // Sum contributions over all frequency bands
+    for (label j=0; j < nLambda_; j++)
+    {
+	// Compute total incident radiation within frequency band
+	tmp<DimensionedField<scalar, volMesh>> Gj
+	(
+	    IRay_[0].ILambda(j)()*IRay_[0].omega()
+	);
 
-    // Only include continuous phase absorption
-    const volScalarField::Internal a =
-        absorptionEmission_->aCont()()();
+	for (label rayI=1; rayI < nRay_; rayI++)
+        {
+	    Gj.ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
+	}
+	
+	Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
+	    - absorptionEmission_->ECont(j)()();
+    }
 
-    return a*G - E;
+    return tRu;
 }
 
 
wideBandfvDOM.patch (6,093 bytes)

henry

2018-03-15 16:25

manager   ~0009428

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

kevnor

2018-03-17 22:19

reporter  

fvDOMTest.tar.gz (3,309 bytes)

kevnor

2018-03-17 22:19

reporter  

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-17 22:29

reporter   ~0009432

I've uploaded a simple test code (hacked together using fireFoam as a base) that computes the explicit energy transfer using radiationModel::Ru and ::Rp and compares it to the explicit energy transfer computed using the sum over all of the rays in the fvDOM model.

I've also attached a radiationProperties file that can be used as a replacement within the oppositeBurningPanels test case which creates two frequency bands with different opacities. Note that this will also need the corrections from https://bugs.openfoam.org/view.php?id=2881 ). Running fireFoam and then the above post-processing code without the patch shows significant energy conservation errors, while the supplied patch eliminates these.

henry

2018-03-21 18:11

manager   ~0009435

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:41

reporter   ~0009451

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:39 kevnor New Issue
2018-03-15 08:39 kevnor File Added: wideBandfvDOM.patch
2018-03-15 16:25 henry Note Added: 0009428
2018-03-17 22:19 kevnor File Added: fvDOMTest.tar.gz
2018-03-17 22:19 kevnor File Added: radiationProperties
2018-03-17 22:29 kevnor Note Added: 0009432
2018-03-21 18:11 henry Note Added: 0009435
2018-03-27 19:41 kevnor Note Added: 0009451
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