CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Understanding the Marshak boundary condition (radiation) (https://www.cfd-online.com/Forums/openfoam-programming-development/135502-understanding-marshak-boundary-condition-radiation.html)

chriss85 May 14, 2014 07:18

Understanding the Marshak boundary condition (radiation)
 
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:
http://combust.hit.edu.cn:8080/fluen...ug/img2712.gif
With e=emissivity at the surface, sigma= Stefan Boltzmann constant and G the spectral intensity integrated over all angles.
And also:

http://combust.hit.edu.cn:8080/fluen...ug/img2683.gif (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:
http://foam.sourceforge.net/docs/cpp/form_22.png
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 04: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 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

chriss85 May 26, 2014 04: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 17:29

Quote:

Originally Posted by elia87 (Post 493757)
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.

chriss85 December 16, 2014 08: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 13: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 06: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 07:16

Hello Joanne,

See my last post in the following thread:

https://www.cfd-online.com/Forums/op...tml#post696795

I hope it can be helpful.

Regards,
Germilly Barreto

chriss85 June 27, 2018 09:06

Hello,

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


All times are GMT -4. The time now is 12:47.