Understanding the Marshak boundary condition (radiation)

chriss85 · May 14, 2014, 08:18

I'm currently working on implementing a non-gray P1 model. This also requires me to modify the Marshak BC.

The problem is that I don't quite understand how it works.
I know that the radiative flux incident to the wall is:

With e=emissivity at the surface, sigma= Stefan Boltzmann constant and G the spectral intensity integrated over all angles.
And also:

(also valid at the wall _w)

The implementation uses a mixed boundary condition, that means it uses an equation of the form a*G+b*dG/dn = c. Does that mean that here this evaluates to:

-e_w/(2*(2-e_w))*G_w + q_rw = -e_w/(2*(2-e_w))*4*sigma*T_w^4
-->
-e_w/(2*(2-e_w))*G_w - Gamma * dG/dn = -e_w/(2*(2-e_w))*4*sigma*T_w^4

I don't quite understand the code of the implementation of the Marshak BC in OF.

The mixed BC fixes the value to this:

with x_p = refValue(), w=valueFraction(), dx/dn = refGrad() and x_c = value in adjacent cell.

In Marshak BC refValue() is 4*sigma*T^4, refGrad() is 0 and valueFraction() is
1.0/(1.0 + Gamma*Delta/(e_w/(2*(2 - e_w))).
I don't understand if and how this matches the mixed equation above. Can anyone shed some light?

In the non-gray model, G is exchanged with G_i, and sigma*T^4 with pi*B_i, where B_i is the integral over the black body spectrum over the frequency range of the band. Does that mean I can simply replace these values in the code and it will be fine? I suppose this should work, but I would like to understand how OpenFOAM implements this.

elia87 · May 23, 2014, 05:42

I think it's okay, because:

$q_{r,w}=-\frac{\epsilon_{w}}{2(2-\epsilon_{w})}(4\sigma T_{w}^4-G_w)$

Assume that

$\theta=\frac{\epsilon_{w}}{2(2-\epsilon_{w})}$

it becomes

$-\theta G_w+q_{r,w}=-\theta 4\sigma T_{w}^4$

you know that $q_r=-\Gamma\nabla G$

from the Fluent Guide, at the boundaries is $q_{r,w}=-\Gamma\frac{\partial G}{\partial n}$

so it becomes

$-\theta G_w+-\Gamma\frac{\partial G}{\partial n}=-\theta 4\sigma T_{w}^4$

$-\theta G_w+-\Gamma\frac{G_w-G_c}{\Delta}=-\theta 4\sigma T_{w}^4$

where $\Delta$ is the distance between face and cell center.

If you explicit

in the equation, you will find

$G_w=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}4\sigma T_{w}^4+\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}G_c$

Here you can find the explanation of the mixed BC and it's easy to recognize that in this case

$refValue=4\sigma T_{w}^4$
$f=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}$

because

$1-\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}=\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}$

the f value (called valueFraction in OpenFOAM) can also be written as

$\frac{1}{1+\frac{\Gamma}{\Delta\theta}}$

In the source code of MarhaskRadiationBC we can find

Code:

refGrad() = 0.0;

Code:

refValue() = 4.0*constant::physicoChemical::sigma.value()*pow4(Tp);

Code:

valueFraction() = 1.0/(1.0 + gamma*patch().deltaCoeffs()/Ep);

where

is what I called $\theta$

Code:

const scalarField Ep(temissivity/(2.0*(2.0 - temissivity)));

and where, according to the post below

$patch().deltaCoeffs()=\frac{1}{\Delta}$

So I think it's ok!

I have some doubts regarding the emissivity I choose in the BC, especially in case of conjugated heat transfer (chtMultiRegion solver) in OpenFOAM 2.3.

If I set (in the coupled wall patch, fluid side obviously)

Code:

emissivityMethod lookup;

I have to set the emissivity value in the same dictionary

Code:

emissivity uniform 0.9;

If I set

Code:

emissivityMethod solidRadiation;

OpenFOAM looks for emissivity value in the radiationProperties file of the coupled solid.
In the AbsorptionEmissionModel dictionary to be precise.

Code:

radiation       on;


radiationModel  opaqueSolid;


absorptionEmissionModel constantAbsorptionEmission;


constantAbsorptionEmissionCoeffs
{
    absorptivity    absorptivity [ 0 -1 0 0 0 0 0 ] 0.0;  //opaque
    emissivity      emissivity [ 0 -1 0 0 0 0 0 ] 0.1;
    E               E [ 1 -1 -3 0 0 0 0 ] 0;
}


scatterModel    none;


sootModel none;

But here I find an emissivity in 1/m, I thought that at the wall the emissivity should have been dimensionless from 0 to 1, like in the formula

$E(T)=\epsilon E_n(T)=\epsilon\sigma T^4$

Anyone can shed more light?

Elia

chriss85 · May 26, 2014, 05:25

Thanks, I think it's clear now how the formula translates to the implementation.
I would also agree with you that the emissivity at the surface should be in [0 .. 1].
Since the opaqueSolid radiation model doesn't do anything, I suspect it's just a different (somewhat unlucky) place to store the value in.

kmefun · December 13, 2014, 18:29

Quote:

Originally Posted by elia87

I think it's okay, because:

$q_{r,w}=-\frac{\epsilon_{w}}{2(2-\epsilon_{w})}(4\sigma T_{w}^4-G_w)$

Assume that

$\theta=\frac{\epsilon_{w}}{2(2-\epsilon_{w})}$

it becomes

$-\theta G_w+q_{r,w}=-\theta 4\sigma T_{w}^4$

you know that $q_r=-\Gamma\nabla G$

from the Fluent Guide, at the boundaries is $q_{r,w}=-\Gamma\frac{\partial G}{\partial n}$

so it becomes

$-\theta G_w+-\Gamma\frac{\partial G}{\partial n}=-\theta 4\sigma T_{w}^4$

$-\theta G_w+-\Gamma\frac{G_w-G_c}{\Delta}=-\theta 4\sigma T_{w}^4$

where $\Delta$ is the distance between face and cell center.

If you explicit

in the equation, you will find

$G_w=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}4\sigma T_{w}^4+\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}G_c$

Here you can find the explanation of the mixed BC and it's easy to recognize that in this case

$refValue=4\sigma T_{w}^4$
$f=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}$

because

$1-\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}=\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}$

the f value (called valueFraction in OpenFOAM) can also be written as

$\frac{1}{1+\frac{\Gamma}{\Delta\theta}}$

In the source code of MarhaskRadiationBC we can find

Code:

refGrad() = 0.0;

Code:

refValue() = 4.0*constant::physicoChemical::sigma.value()*pow4(Tp);

Code:

valueFraction() = 1.0/(1.0 + gamma*patch().deltaCoeffs()/Ep);

where

is what I called $\theta$

Code:

const scalarField Ep(temissivity/(2.0*(2.0 - temissivity)));

and where, according to the post below

$patch().deltaCoeffs()=\frac{1}{\Delta}$

So I think it's ok!

I have some doubts regarding the emissivity I choose in the BC, especially in case of conjugated heat transfer (chtMultiRegion solver) in OpenFOAM 2.3.

If I set (in the coupled wall patch, fluid side obviously)

Code:

emissivityMethod lookup;

I have to set the emissivity value in the same dictionary

Code:

emissivity uniform 0.9;

If I set

Code:

emissivityMethod solidRadiation;

OpenFOAM looks for emissivity value in the radiationProperties file of the coupled solid.
In the AbsorptionEmissionModel dictionary to be precise.

Code:

radiation       on;


radiationModel  opaqueSolid;


absorptionEmissionModel constantAbsorptionEmission;


constantAbsorptionEmissionCoeffs
{
    absorptivity    absorptivity [ 0 -1 0 0 0 0 0 ] 0.0;  //opaque
    emissivity      emissivity [ 0 -1 0 0 0 0 0 ] 0.1;
    E               E [ 1 -1 -3 0 0 0 0 ] 0;
}


scatterModel    none;


sootModel none;

But here I find an emissivity in 1/m, I thought that at the wall the emissivity should have been dimensionless from 0 to 1, like in the formula

$E(T)=\epsilon E_n(T)=\epsilon\sigma T^4$

Anyone can shed more light?

Elia

Generally speaking, the emissivity can change depending on the wavelength. That's why there is a per unit wavelength for its unit.

chriss85 · December 16, 2014, 09:38

But a wavelength-dependent emissivity still is dimensionless, as it's not a density like a spectral intensity for example.

Also, the emissivity in radiationProperties is used in two different contexts in the models, which is somewhat confusing.

Germilly · June 19, 2018, 14:08

Hello,

I was searching in the forum and this thread seems to be the right place to post my question.

I'm modelling radiation transfer inside a pipe which contain a participating media. My problem is about the boundary condition for the incident radiation G (I'm using the P1 radiationModel). In that pipe, I have 3 boundary patches:

1) "inlet" of the pipe: Which is in direct contact with the environment
1) "outlet" of the pipe: also in direct contact with the environment
3) "wall" (cylindrical surface): Bounded by a solid diffuse surface

For the wall, I'm using in the 0/G file the MarshakRadiation boundary condition:

Code:

wall
{
        type            MarshakRadiation;
        T               Ts; // Ts is the name of my temperature field
        value           uniform 0;;
 }

I think I have understood the parameters in the constant/boundaryRadiationProperties for the wall:

Code:

wall
{
        type              boundaryRadiation;
        mode              lookup;
        emissivity        uniform 0.8; // The emissivity of the wall
        value             uniform 0;
}

I also think I have understood the parameters in constant/radiationProperties which characterize the participating media.

But, for the inlet and outlet, I dont know what should I do.
I cannot define an emissivity for the inlet and outlet because there are no solid surface, there are only the participating media.

Do you know what should I define for these two boundary conditions (inlet and outlet)?

Thank you

Best regards,
Germilly Barreto

Joanne · June 27, 2018, 07:29

Hi Germilly,

I am interested in a similar problem as yours, I am unsure which BC to use for flow inlets/outlets when implementing a radiation model.

Please let me know if you make any discoveries, and I will do likewise.

Kind regards,
Joanne

Germilly · June 27, 2018, 08:16

Hello Joanne,

See my last post in the following thread:

Radiation boundary conditions for flow through boundaries in openFoam

I hope it can be helpful.

Regards,
Germilly Barreto

chriss85 · June 27, 2018, 10:06

Hello,

I replied to the other thread, I hope this helps somewhat.

May 14, 2014, 08:18	Understanding the Marshak boundary condition (radiation)	#1
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	I'm currently working on implementing a non-gray P1 model. This also requires me to modify the Marshak BC. The problem is that I don't quite understand how it works. I know that the radiative flux incident to the wall is: With e=emissivity at the surface, sigma= Stefan Boltzmann constant and G the spectral intensity integrated over all angles. And also: (also valid at the wall _w) The implementation uses a mixed boundary condition, that means it uses an equation of the form aG+bdG/dn = c. Does that mean that here this evaluates to: -e_w/(2(2-e_w))G_w + q_rw = -e_w/(2(2-e_w))4sigmaT_w^4 --> -e_w/(2(2-e_w))G_w - Gamma * dG/dn = -e_w/(2(2-e_w))4sigmaT_w^4 I don't quite understand the code of the implementation of the Marshak BC in OF. The mixed BC fixes the value to this: with x_p = refValue(), w=valueFraction(), dx/dn = refGrad() and x_c = value in adjacent cell. In Marshak BC refValue() is 4sigmaT^4, refGrad() is 0 and valueFraction() is 1.0/(1.0 + GammaDelta/(e_w/(2(2 - e_w))). I don't understand if and how this matches the mixed equation above. Can anyone shed some light? In the non-gray model, G is exchanged with G_i, and sigmaT^4 with piB_i, where B_i is the integral over the black body spectrum over the frequency range of the band. Does that mean I can simply replace these values in the code and it will be fine? I suppose this should work, but I would like to understand how OpenFOAM implements this. atulkjoy likes this.

May 23, 2014, 05:42		#2
elia87 New Member Elia Agnani Join Date: Oct 2011 Location: Modena, Italy Posts: 5 Rep Power: 15	I think it's okay, because: $q_{r,w}=-\frac{\epsilon_{w}}{2(2-\epsilon_{w})}(4\sigma T_{w}^4-G_w)$ Assume that $\theta=\frac{\epsilon_{w}}{2(2-\epsilon_{w})}$ it becomes $-\theta G_w+q_{r,w}=-\theta 4\sigma T_{w}^4$ you know that $q_r=-\Gamma\nabla G$ from the Fluent Guide, at the boundaries is $q_{r,w}=-\Gamma\frac{\partial G}{\partial n}$ so it becomes $-\theta G_w+-\Gamma\frac{\partial G}{\partial n}=-\theta 4\sigma T_{w}^4$ $-\theta G_w+-\Gamma\frac{G_w-G_c}{\Delta}=-\theta 4\sigma T_{w}^4$ where $\Delta$ is the distance between face and cell center. If you explicit in the equation, you will find $G_w=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}4\sigma T_{w}^4+\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}G_c$ Here you can find the explanation of the mixed BC and it's easy to recognize that in this case $refValue=4\sigma T_{w}^4$ $f=\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}$ because $1-\frac{\theta}{\theta+\frac{\Gamma}{\Delta}}=\frac{\frac{\Gamma}{\Delta}}{\theta+\frac{\Gamma}{\Delta}}$ the f value (called valueFraction in OpenFOAM) can also be written as $\frac{1}{1+\frac{\Gamma}{\Delta\theta}}$ In the source code of MarhaskRadiationBC we can find Code: refGrad() = 0.0; Code: refValue() = 4.0constant::physicoChemical::sigma.value()pow4(Tp); Code: valueFraction() = 1.0/(1.0 + gammapatch().deltaCoeffs()/Ep); where is what I called $\theta$ Code: const scalarField Ep(temissivity/(2.0(2.0 - temissivity))); and where, according to the post below $patch().deltaCoeffs()=\frac{1}{\Delta}$ So I think it's ok! I have some doubts regarding the emissivity I choose in the BC, especially in case of conjugated heat transfer (chtMultiRegion solver) in OpenFOAM 2.3. If I set (in the coupled wall patch, fluid side obviously) Code: emissivityMethod lookup; I have to set the emissivity value in the same dictionary Code: emissivity uniform 0.9; If I set Code: emissivityMethod solidRadiation; OpenFOAM looks for emissivity value in the radiationProperties file of the coupled solid. In the AbsorptionEmissionModel dictionary to be precise. Code: radiation on; radiationModel opaqueSolid; absorptionEmissionModel constantAbsorptionEmission; constantAbsorptionEmissionCoeffs { absorptivity absorptivity [ 0 -1 0 0 0 0 0 ] 0.0; //opaque emissivity emissivity [ 0 -1 0 0 0 0 0 ] 0.1; E E [ 1 -1 -3 0 0 0 0 ] 0; } scatterModel none; sootModel none; But here I find an emissivity in 1/m, I thought that at the wall the emissivity should have been dimensionless from 0 to 1, like in the formula $E(T)=\epsilon E_n(T)=\epsilon\sigma T^4$ Anyone can shed more light? Elia francescomarra, jherb, Tobi and 12 others like this. __________________ SnappyWiki

June 19, 2018, 14:08		#6
Germilly New Member Germilly Barreto Join Date: Jul 2016 Location: Portugal Posts: 25 Rep Power: 10	Hello, I was searching in the forum and this thread seems to be the right place to post my question. I'm modelling radiation transfer inside a pipe which contain a participating media. My problem is about the boundary condition for the incident radiation G (I'm using the P1 radiationModel). In that pipe, I have 3 boundary patches: 1) "inlet" of the pipe: Which is in direct contact with the environment 1) "outlet" of the pipe: also in direct contact with the environment 3) "wall" (cylindrical surface): Bounded by a solid diffuse surface For the wall, I'm using in the 0/G file the MarshakRadiation boundary condition: Code: wall { type MarshakRadiation; T Ts; // Ts is the name of my temperature field value uniform 0;; } I think I have understood the parameters in the constant/boundaryRadiationProperties for the wall: Code: wall { type boundaryRadiation; mode lookup; emissivity uniform 0.8; // The emissivity of the wall value uniform 0; } I also think I have understood the parameters in constant/radiationProperties which characterize the participating media. But, for the inlet and outlet, I dont know what should I do. I cannot define an emissivity for the inlet and outlet because there are no solid surface, there are only the participating media. Do you know what should I define for these two boundary conditions (inlet and outlet)? Thank you Best regards, Germilly Barreto atulkjoy, Joanne, altinel and 1 others like this. Last edited by Germilly; June 21, 2018 at 14:31.

June 27, 2018, 07:29		#7
Joanne New Member Joanne Join Date: Aug 2017 Location: Ireland Posts: 6 Rep Power: 9	Hi Germilly, I am interested in a similar problem as yours, I am unsure which BC to use for flow inlets/outlets when implementing a radiation model. Please let me know if you make any discoveries, and I will do likewise. Kind regards, Joanne Last edited by Joanne; June 27, 2018 at 07:30. Reason: grammar

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Low Mixing time Problem	Mavier	CFX	5	April 29, 2013 01:00
Velocity profile boundary condition	Tuca	FLOW-3D	1	April 23, 2013 13:02
External Radiation Boundary Condition (Two sided wall), Grid Interface	CFD XUE	FLUENT	0	July 8, 2010 07:49
[Commercial meshers] Trimmed cell and embedded refinement mesh conversion issues	michele	OpenFOAM Meshing & Mesh Conversion	2	July 15, 2005 05:15
The Boundary Condition about the Flat Plate	boing	Main CFD Forum	1	January 6, 2002 17:53

May 26, 2014, 05:25		#3
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	Thanks, I think it's clear now how the formula translates to the implementation. I would also agree with you that the emissivity at the surface should be in [0 .. 1]. Since the opaqueSolid radiation model doesn't do anything, I suspect it's just a different (somewhat unlucky) place to store the value in.

December 16, 2014, 09:38		#5
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	But a wavelength-dependent emissivity still is dimensionless, as it's not a density like a spectral intensity for example. Also, the emissivity in radiationProperties is used in two different contexts in the models, which is somewhat confusing.

June 27, 2018, 08:16		#8
Germilly New Member Germilly Barreto Join Date: Jul 2016 Location: Portugal Posts: 25 Rep Power: 10	Hello Joanne, See my last post in the following thread: Radiation boundary conditions for flow through boundaries in openFoam I hope it can be helpful. Regards, Germilly Barreto

June 27, 2018, 10:06		#9
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	Hello, I replied to the other thread, I hope this helps somewhat.