CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Understanding the Marshak boundary condition (radiation)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree16Likes
  • 1 Post By chriss85
  • 13 Post By elia87
  • 2 Post By Germilly

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 14, 2014, 07:18
Default Understanding the Marshak boundary condition (radiation)
  #1
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
atulkjoy likes this.
chriss85 is offline   Reply With Quote

Old   May 23, 2014, 04:42
Default
  #2
New Member
 
Elia Agnani
Join Date: Oct 2011
Location: Modena, Italy
Posts: 5
Rep Power: 11
elia87 is on a distinguished road
Send a message via Skype™ to 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 G_w 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

refGrad=0
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 E_p 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
__________________
SnappyWiki
elia87 is offline   Reply With Quote

Old   May 26, 2014, 04:25
Default
  #3
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
chriss85 is offline   Reply With Quote

Old   December 13, 2014, 17:29
Default
  #4
Member
 
Kaufman
Join Date: Jul 2013
Posts: 55
Rep Power: 10
kmefun is on a distinguished road
Quote:
Originally Posted by elia87 View Post
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 G_w 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

refGrad=0
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 E_p 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.
kmefun is offline   Reply With Quote

Old   December 16, 2014, 08:38
Default
  #5
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
chriss85 is offline   Reply With Quote

Old   June 19, 2018, 13:08
Default
  #6
New Member
 
Germilly Barreto
Join Date: Jul 2016
Location: Portugal
Posts: 24
Rep Power: 7
Germilly is on a distinguished road
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 and Joanne like this.

Last edited by Germilly; June 21, 2018 at 13:31.
Germilly is offline   Reply With Quote

Old   June 27, 2018, 06:29
Default
  #7
New Member
 
Joanne
Join Date: Aug 2017
Location: Ireland
Posts: 6
Rep Power: 6
Joanne is on a distinguished road
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 06:30. Reason: grammar
Joanne is offline   Reply With Quote

Old   June 27, 2018, 07:16
Default
  #8
New Member
 
Germilly Barreto
Join Date: Jul 2016
Location: Portugal
Posts: 24
Rep Power: 7
Germilly is on a distinguished road
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
Germilly is offline   Reply With Quote

Old   June 27, 2018, 09:06
Default
  #9
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
Hello,

I replied to the other thread, I hope this helps somewhat.
chriss85 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 15:22.