Implementation of The Rosseland and The P1 Radiation Models in The System of Navier-Stokes Equations With The Boundary Element Method
Implementation of The Rosseland and The P1 Radiation Models in The System of Navier-Stokes Equations With The Boundary Element Method
net/publication/315732794
CITATION READS
1 971
4 authors, including:
Matjaz Hriberšek
University of Maribor
105 PUBLICATIONS 678 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
A numerical model of translational and rotational momentum transfer of small non-spherical rigid particles in fluid dominated two-phase flows View project
All content following this page was uploaded by L. Škerget on 14 August 2017.
ABSTRACT
The objective of this article is to develop a boundary element numerical model to solve coupled prob-
lems involving heat energy diffusion, convection and radiation in a participating medium. In this
study, the contributions from radiant energy transfer are presented using two approaches for optical
thick fluids: the Rosseland diffusion approximation and the P1 approximation. The governing Navier–
Stokes equations are written in the velocity–vorticity formulation for the kinematics and kinetics of the
fluid motion. The approximate numerical solution algorithm is based on a boundary element numeri-
cal model in its macro-element formulation. Validity of the proposed implementation is tested on a
one-dimensional test case using a grey participating medium at radiative equilibrium between two
isothermal black surfaces.
Keywords: compressible fluid flow, radiation models, boundary element method
1 INTRODUCTION
The Navier–Stokes equations set is commonly used as a frame for the solution of transport
phenomena in a fluid flow. It provides a mathematical model of physical conservation laws
of mass, momentum and energy considering specific rheological models describing non-con-
vective fluxes of momentum and energy. In general, all three physically different mechanisms
of heat transport can occur, that is diffusion, convection and radiation. The energy radiation
phenomenon, which is a complex non-linear mode of heat transfer, gains importance at suf-
ficiently high temperature [1]. At temperatures which are high enough these processes are
essentially interdependent; energy transfer by one mechanism can influence heat exchange
by the other mechanism and vice versa. The objective of this article is to develop a boundary
element numerical simulation model to solve coupled problems involving heat energy diffu-
sion, convection and radiation in a participating viscous compressible fluid flow.
The governing equation for radiative heat transfer is the radiative transfer equation [2],
which is based on an energy balance for radiation passing through a differential volume in a
participating medium in local thermo-dynamic equilibrium.
The radiation impact on overall heat transfer is conveyed in the energy equation, where,
in the non-convective energy flux, besides the diffusion heat flux the radiative heat flux also
needs to be taken into consideration. This is done in such a way that we include the term for
the divergence of the radiative flux vector into the energy equation as the radiative energy
source [3, 4]. The radiative transfer equations (RTE) is an integro-differential equation
presenting a serious issue in computational fluid dynamics. Applying the chosen radiation
model means, under the given physical circumstances, a simplification of the radiative trans-
fer equation. In this study, the contributions from radiant energy transfer are presented using
two approaches for optical thick fluids, that is the Rosseland diffusion approximation and the
P1 approximation.
2 GOVERNING EQUATIONS
The analytical description of the motion of a continuous viscous compressible heat radiation
semi-transparent fluid is based on the conservation of mass, momentum and heat energy
with associated rheological models for the non-convective fluxes of the momentum and heat
energy and equations of state. The present development is focused on the laminar flow of
compressible isotropic radiation semi-transparent fluid in solution domain R = Ω × T , where
Ω stands for the two-dimensional plane domain bounded by boundary Γ defined by the
outward-pointing unit normal n, whilst T represents the time dimension of the transport
phenomenon.
The field functions of interest are the velocity vector field ui(rj, t), scalar pressure field p(rj, t),
temperature field T(rj, t) and the field of mass density ρ(rj,t), so that the mass, momentum and
energy equations are given by the following set of non-linear equations:
∂u j 1 Dr
=− = D, (1)
∂x j r Dt
Dui ∂tij ∂p
r =− − + r gi (2)
Dt ∂x j ∂x i
DT ∂q Dj ∂q Rj
c =− − (3)
Dt ∂x j ∂x j
in the Cartesian frame xi , where ρ and c denote changeable mass density and isobaric spe-
cific heat capacity per unit volume, c = cp ρ, t is the time, gi is the gravitational acceleration
vector and tij represents the tensor components of the momentum diffusion, whilst the vector
R
variables q D and q j are heat diffusion and radiation fluxes, respectively. The differential
j
operator D ( ⋅) / Dt stands for the Stokes material derivative.
The conservation eqns (2) and (3) contain two molecular diffusive fluxes, that is tij and q Dj ,
representing the diffusion of linear momentum and heat energy, respectively. The New-
ton linear momentum diffusion constitutive model for compressible viscous shear fluid is
considered, such as
2
tij = 2 heij − hD dij , (4)
3
where D = diυυ =∈ ii represents the divergence of the velocity field or local expansion field,
and η is a dynamic viscosity. For most heat transfer problems of practical importance, the
simplification known as the Fourier law of heat diffusion is accurate enough, namely
350 P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017)
∂T
qiD = −k D , (5)
∂xi
where k D is the thermal heat conductivity.
All bodies at absolute temperature T emit electromagnetic radiation continuously over a
wide range of wavelengths. At temperatures which are high enough, the simulation of the
heat transfer becomes very complex. The mechanism of the heat transfer plays an important
role in radiation which can present a great deal of the total heat flux. The governing equation
for radiative heat transfer is the radiative transfer equation (RTE) [2], which is based on an
energy balance for radiation passing through a differential volume in a participating medium
in local thermodynamic equilibrium (LTE). The change in spectral intensity iλ (r ) along a
path from r to r + dr , where the time dependence of the intensity is neglected, expresses the
quasi-steady form of the RTE
4p
dil s
dl
= ∇il ⋅ i = al il b − K l il + Sl
4p ∫
w =0
il fl d w (6)
The spectral extinction coefficient K l (r ) = al (r ) + sSl (r ) is defined as the sum of the
spectral absorption coefficient aλ and the spectral scattering coefficient σ Sλ , iλ b is the black
body isotropic spectral intensity and ϕl is the scattering phase function. The eqn (6) is a
first-order integro-differential equation for iλ in a fixed direction r . Due to the dependence
on three spatial coordinates, two local direction coordinates and wavelength, an analytical
solution is almost impossible for most engineering applications.
Thus, the eqn (6) has to be solved numerically using radiation transport models for spatial
and directional dependencies and spectral models for the spectral dependency. In this study,
we present an analysis of two common approximations for modelling radiative heat transfer
that occurs in optically thick fluids. The contributions from the radiative energy transfer are
presented using two approaches: the Rosseland diffusion approximation and the P1 approx-
imation.
4π ∂ib (T ) 16 n2σ T 3 ∂T ∂T
qiR (rj ) = − =− = −k R (7)
3K R ∂xi 3K R ∂xi ∂xi
where K R is the Rosseland mean extinction coefficient, n is the refractive index and σ is the
Stefan–Boltzmann constant. Although the Rosseland model provides a substantial simpli-
fication of the RTE
L
and is recommended for use in problems where the optical differential
thickness kl = ∫0 K l dr exceeds 10 [2], it is often used for simulation of the radiation processes
in many engineering applications. In analogy to eqn (5) the radiation heat conductivity k R is
introduced in eqn (7).
P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017) 351
For the Rosseland radiation model it is possible to specify an adiabatic boundary condition,
that is a zero-temperature gradient dT / dn = 0 at the solid wall or a specific wall temperature
as the Dirichlet boundary condition. However, near a boundary, the diffusion approximation
may not be accurate as the radiation is not isotropic [2]. To overcome this difficulty, the
boundary condition at the edge of the medium is modified by using the effective jump bound-
ary condition. Using the Deissler jump boundary condition concept for pure radiation [5],
Goldstein and Howell introduced a similar concept for combined conduction and radiation
[6]. In this model, the radiative heat flux at the wall boundary qwR is defined using the jump
coefficient
s Tw4 − T 4 ( x → 0 )
Ψw = (8)
qwR
where Tw is the wall temperature and T ( x → 0) is the extrapolated temperature of the medium
at the wall. The jump coefficient is a function only of the conduction-radiation parameter
N w = k D K R / 4sTw3 , which expresses a measure of the ratio of the energy transferred by
conduction and radiation. For large N w the jump effect can be neglected, as heat conduction
dominates over radiation effects near the wall. In general, the jump coefficient is approxi-
mated by a curve fit to the plot given in [3]
3 1/ 2 for N w ≤ 0.01,
2 x + 3 x 2 − 12 x + 7
Ψw = for 0.01 ≤ N w ≤ 10, (9)
54
0 for N w ≥ 10,
∂T s T 4 ( x → 0) − Tw4
qwR = − k R = (10)
∂n w Ψw
As shown in the derivation of Ψ w , the conditions for which the diffusion model
is valid lead to the temperature jump Tw − T ( x → 0 ) being small, so the difference
( )
Tw4 − T 4 ( x → 0 ) ≈ 4Tw3 Tw − T ( x → 0 ) can be linearized. The jump temperature follows
the relation
4 Ψ w ∂T
Tw = T ( x → 0) − . (11)
3K R ∂ n w
case if only four terms in the series are retained. The method is a generalization of the Milne–
Eddington equations analysed in [7, 2]. In engineering radiative transfer problems, the P1
model should typically be used for spectral optical thickness κ λ > 1 [3].
A medium with spectral extinction coefficient K l (r ) and isotropic scattering is consid-
ered. The spectral incident radiation at position r is defined as
4p
Gl (r ) = ∫ i (r , l )d w, (12)
0
l
where the integration takes over all solid angles [2]. Note that the spectral incident radiation
divided by the speed of light Gλ / c is the spectral radiative energy density at location r in the
radiation field. The P1 radiation model yields two spatial differential governing equations,
one for the gradient of the directionally averaged spectral intensity
1 1
qλR = − ∇Gλ = − ∇Gλ = −Γ λ ∇Gλ , (13)
3K λ 3 aλ + σ Sλ
where the parameter Γ λ = 1 3K λ , and another for the divergence of net radiative heat flux
density vector
∇ ⋅ qλR = aλ 4π iλ b − Gλ . (14)
Equations (13) and (14) can be combined to yield a second-order elliptic PDE for the incident
radiation.
∇ ⋅ qλR = −∇ ⋅ Γ λ ∇Gλ = aλ 4π iλ b − Gλ . (15)
This equation is simply a statement that the net radiative heat flux out of any region occupied
by the medium is the difference between that emitted and that absorbed in the volume under
consideration.
In the grey medium with constant absorption and extinction coefficients eqn (15) is simpli-
fied to a non-linear inhomogeneous modified Helmholtz equation:
∂2 G
− bG + b = 0 with b = 3aK and b = 4 bsT 4 , (16)
∂x j ∂x j
whilst the divergence of the radiation flux vector in eqn (3) can be expressed as the local
radiation source term S R:
∂q Rj
∂x j
( )
= − a G − 4sT 4 = − S R . (17).
If it is assumed that the walls are diffuse grey surfaces, the eqn (16) is solved using Marshak
boundary condition [3]
−Γ
∂G
∂n w
=
ew
2 ( 2 − ew )
( )
4sTw4 − Gw , (18)
P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017) 353
where ∈w is the emissivity of the wall and the subscript w denotes the value of the indicated
variable at the wall. Unlike the Rosseland diffusion approximation discussed above, there is
no ambiguity about boundary conditions for the P1 approximation.
For the coupling of the radiative heat transport with the fluid dynamics, LTE is assumed
and the time dependence of the radiative transfer equation is neglected. It follows from LTE
that the temperature of the fluid and the corresponding radiative temperature in the medium
are equal.
For the known vorticity and local expansion field functions, the corresponding velocity
vector can be determined by solving eqn (19), provided that appropriate boundary condi-
tions for the velocity are prescribed, that is the normal and tangential component of the
velocity vector. The kinetic aspect of the fluid motion is governed by the vorticity transport
equation.
∂wi ∂u j wi ∂ 2 wi ∂w j ui 1 ∂r gk 1 ∂f m
+ = no + + eijk + eijk k (20)
∂t ∂x j ∂x j ∂x j ∂x j ro ∂x j ro ∂x j
describing the redistribution of the vorticity in the fluid domain by different transport phe-
nomena, for example diffusion, convection, twisting and stretching, whilst the buoyancy,
compressibility, and the non-linear terms act as a source or strengthen terms. The vorticity
transport equation is a highly non-linear partial differential equation due to the products
of velocity and vorticity in convective and in stretching-twisting terms, and the velocity
field function is kinematically dependent on vorticity and local expansion. However, strong
coupling of the kinematics and kinetics can be clearly observed, even in the case of incom-
pressible fluid.
The energy conservation equations are
DT ∂ ∂T R
c = keff + S , (21)
Dt ∂x j ∂x j
∂2G
− β G + b = 0, (22)
∂x j ∂x j
where keff = k D + k R and S R = 0 for the Rosseland radiation model and keff = k D and
( )
S R = − a G − 4sT 4 for the P1 one, respectively.
354 P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017)
The pseudo body force term f m and pseudo heat source term STm were introduced into
the vorticity fransport eqn (20) and into energy eqn (21) respectively, capturing the variable
transport property effects, and given by expressions
4
f m = −∇ × ( hw
) + 2∇h × w + 2∇u ⋅∇h + ∇ ( hD ) − 2 D∇h − r a, (23)
3
while the pseudo heat source term is given by an expression
DT
STm = ∇ k∇T − c
( . (24) )
Dt
∂2T ∂T
L [T ] + b = ao − + b = 0, (25)
∂x j ∂x j ∂t
∫
+ TF −1u*F −1 d Ω.
Ω
The boundary integrals describe the total heat flux on the boundary due to molecular diffusion
and convection. The first domain integral gives the influence of the perturbated convection
and the nonlinear diffusion flux, the second domain integral includes the non-linear material
effects and radiation source, while the last domain integral is due to the initial temperature
distribution effect on the development of the temperature field in subsequent time interval.
The incident radiation equation is an elliptic modified Helmholtz equation and, therefore,
by employing the linear elliptic modified Helmholtz differential operator, we obtain the
following expression
∂2G
L[G ] + b = − β G + b = 0 (27)
∂x j ∂x j
P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017) 355
and the corresponding boundary-domain integral representation to eqn (22) can be stated as
∂G
∫
c (x ) G (x ) + Gq * d Γ = ∫ ∂n u * d Γ + ∫ 4 bsT u * d Ω (28)
4
Γ Γ Ω
where u∗ is now the modified Helmholtz fundamental solution given by
u* =
1
2p
K0 ( br ) and q* =
dn
2pr 2
b rK1 ( br ) , (29)
whilst Kα is a modified Bessel function of the second kind of order α .
6 VALIDATION
To check the validity of the implemented Rosseland and the P1 radiation models we inves-
tigated a one-dimensional case. We consider a grey participating medium at radiative
equilibrium between two isothermal black surfaces (∈= 1) at temperatures Th = 600 K and
Tc = 300 K . We consider radiation as being coupled to the energy equation via Rosseland
diffusion approximation with specified Dirichlet boundary condition and radiative heat
transfer as a source term using the P1 radiation model with specified Marshak boundary
condition.
The coefficient of the diffusion thermal heat conductivity is reduced to the value where
k D = 10 −4 W mK . The test example is analysed for a grey medium with optical thickness
kL = aL = 10 and κ L = 2. The influence of the natural convection was neglected. The exact
results of the RTE are available for these cases and can be found in [1]. Figures 1 and 2 show
the non-dimensional temperature T ∗ = T 4 − Tc4
( Th4 − Tc4 versus non-dimensional coor-
)( )
dinate x∗ = x L and compare the results of the P1 and the Rosseland radiation models with
the results of [1].
The results reveal good agreement between present numerical results and the exact solu-
tion of Modest [1]. Results also reveal a temperature discontinuity (a sharp temperature
profile) at the walls. In a limiting case of a transparent medium κ L → 0, the non-dimensional
temperature takes the value of 0.5. The temperature slip at the walls decreases as the optical
thickness increases and vanishes as κ L → ∞.
356 P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017)
Figure 1:
The comparison of pure radiation simulation results for the dimensionless
temperature versus non-dimensional coordinate for a grey medium between
two plates. Results of the P1 model are shown versus benchmark results of
Modest [1].
Figure 2:
The comparison of pure radiation simulation results for the dimensionless
temperature versus non-dimensional coordinate for a grey medium between two
plates. Results of the Rosseland model are shown versus benchmark results of
Modest [1].
P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017) 357
While the Rosseland model reveals the temperature discontinuity at the walls, the P1 model,
due to low heat conductivity, results in a more physically realistic temperature profile at the wall.
7 CONCLUSIONS
We have implemented two approaches to solve for radiation energy transport within an optically
thick fluid. The choice of any one of the solution methods will depend upon the computational
effort needed for the solution as well as for the accuracy of the solution. Based upon the results
of this study, there are some clear recommendations. The Rosseland approximation is easy to
implement and economical in computing needs. These issues are especially important with
regard to incorporating internal radiant heat transfer into existing heat transfer codes. However,
this simplicity carries with it a significant source of inaccuracy – the inability of this method
to capture the physics of the thermal boundary layer near the walls. Indeed, this downfall may
lead to significant errors, especially in systems such as this where the thermal boundary layers
are important for driving flow. The P1 model is slightly more difficult to implement, since it
requires solving an additional coupled partial differential equation for an additional field var-
iable, the incident radiation. In addition, the computational effort needed to solve with the P1
model is slighter greater than that needed for the Rosseland diffusion methods.
REFERENCES
[1] Modest, M.F., Radiative Heat Transfer, 3rd edn., Academic Press, 2013.
[2] Seaid, M., Klar, A. & Pinnau, R., Numerical solvers for radiation and conduction in
high temperature gas flows. Journal of Flow Turbulence and Combustion, 75(1), pp.
173–190, 2005.
[Link]
[3] ANSYS CFX-Solver: Theory Guide, 1996-2006 ANSYS Europe, Ltd., ANSYS CFX
Release 11.0, 2006.
[4] Dubroca, B., Seaid, M. & Teleaga, I., A consistent approach for the coupling of radia-
tion and hydrodynamics at low mach number. Journal of Computational Physics,
225(1), pp. 1039–1065, 2007.
[Link]
[5] Deissler, R.G., Diffusion approximation for thermal radiation in gases with jump
boundary condition. Journal of Heat Transfer, 86(2), pp. 240–246, 1964.
[Link]
[6] Goldstein, M.E. & Howell, J.R., Boundary Conditions for the Diffusion Solution of
Coupled Conduction-Radiation Problems, NASA Technical Note: TN D-4618, 1968.
[7] Liu, X.L., Gong, G.C. & Cheng, H.S., Combined natural convection and radiation heat
transfer of various absorbing-emitting-scattering media in a square cavity. Advances in
Mechanical Engineering, 6, p. 403690, 2014.
[Link]
[8] Skerget, L. & Ravnik, J., BEM simulation of compressible fluid flow in an enclosure
induced by thermoacoustic waves. Engineering Analysis with Boundary Elements,
33(4), pp. 561–571, 2009.
[Link]
[9] HriberSsek, M. & Skerget, L., Iterative methods in solving Navier-Stokes equations by
the boundary element method. International Jouenal for Numerical Methods in Engi-
neering, 39(1), pp. 115–139, 1996.
[Link]
NME852>[Link];2-D
358 P. Crnjac, et al., Int. J. Comp. Meth. and Exp. Meas., Vol. 5, No. 3 (2017)
[10] RamSsak, M. & Skerget, L., A subdomain boundary element method for high-Reynolds
laminar flow using stream function- vorticity formulation. International Journal for
Numerical Methods in Fluids, 46(8), pp. 815–847, 2004.
[Link]
[11] Popov, V., Power, H. & Skerget, L., Domain Decomposition Techniques for Bound-
ary Elements, Application to Fluid Flow, Advances in Boundary Element Series, WIT
Press: Southampton and Boston, 2007.
[Link]