7.7 Cavity Radiation

Back

In the present section, we examine what happens if radiation is exchanged among several

surfaces. Radiation is a rather complicated subject meriting a much more extensive treatise.

For more details, the reader is referred to (Incropera and DeWitt 2002).

7.7.1 Governing equations

Consider a differential surface dA1 emitting radiation toward a differential surface dA2 at

a distance r (Figure 7.2). The surface dA2 is perpendicular to the line connecting dA1 with

dA2 and covers a spatial angle dω satisfying

dω = dA2

r2 . (7.95)

x

r

y

z

dA1

dA2

ϕ

ψ

dω

Figure 7.2 Radiation of surface dA1 onto dA2

318 HEAT TRANSFER

x

r

y

z

ψ

dφ

dψ

r sinψ

Figure 7.3 Infinitesimal surface element

The spectral intensity IE in a certain direction (ϕ,ψ) is defined as the radiation power

per unit solid angle dω about this direction, per unit wavelength dλ, per unit emitting area

perpendicular to this direction:

IE = dPE

dω dλ dA1 cosψ

. (7.96)

Accordingly, for the radiation power between two infinitesimal areas dA1 and dA2, the

area dA1 enters in IE in the form of the projected area dA1 cosψ, whereas dA2 enters in

the form of the spatial angle dω. Furthermore, IE depends on the wavelength of emission.

The spectral, hemispherical emissive power Eλ, is defined as the radiation power in all

directions of a hemisphere per unit wavelength dλ per unit emitting area (not projected!).

Hence,

Eλ = _

hemisphere

IE cosψ dω. (7.97)

An infinitesimal solid angle can be written as (Figure 7.3)

dω = sinψ dϕ dψ (7.98)

leading to

Eλ = _ 2π

0

_ π/2

0

IE cosψ sinψ dψ dϕ. (7.99)

HEAT TRANSFER 319

If the emission does not depend on the direction (ϕ,ψ), it is called diffuse emission. Here

and in the section that follows, we assume that we deal with diffuse emitters. In that case,

IE is no function of ϕ and ψ, and Equation (7.99) reads

Eλ = IE

_ 2π

0

_ π/2

0

cosψ sinψ dψ dϕ = πIE. (7.100)

A special kind of diffuse emitter is a blackbody. Its properties are as follows:

1. It emits diffuse, that is, the spectral intensity only depends on the wavelength and

temperature, not on the emission angle.

2. No body can emit more energy than a blackbody for a given wavelength and temperature.

3. All incident radiation is completely absorbed, no reflection takes place.

A blackbody is classically symbolized by a cavity at a uniform temperature with a small

aperture. The spectral intensity of blackbody radiation was first determined by Planck, and

satisfies

IE,b =

2hc2

0

λ5[exp(hc0/λkθ) − 1]

(7.101)

where h = 6.6256 × 10−34 Js is the Planck constant, k = 1.3805 × 10−23 J/K is the Boltzmann

constant, c0 = 2.998 × 108 m/s is the speed of light in vacuum and θ is the temperature

of the blackbody in Kelvin. Since a blackbody is a diffuse emitter, one obtains for

the spectral emissive power

Eλ,b = πIE,b. (7.102)

The total emissive power is the power emitted per unit of emitting area and satisfies

Eb = _

∞

0

Eλ,b dλ. (7.103)

Substituting Equations (7.101) and (7.102) into Equation (7.103) and performing the integration,

one obtains the Stefan–Boltzmann law

Eb = σθ4 (7.104)

where σ = 5.67 × 10−8 W/m2K4 is the Stefan–Boltzmann constant.

The blackbody is an ideal emitter. Real bodies will emit less. The spectral, directional

emissivity is defined as the ratio of the real spectral, directional radiation intensity to the

spectral blackbody intensity at the same temperature:

_λ,ω := IE

IE,b

. (7.105)

Here, we will assume to deal with diffuse emitters and work with averages over all wavelengths.

Therefore, we define the total hemispherical emissivity as the ratio of the total

emissive power to the emissive power of a blackbody at the same temperature:

_(θ) := E

Eb

. (7.106)

320 HEAT TRANSFER

The total emissive power is a function of the radiating surface and the temperature. Using

Equations (7.96), (7.97), (7.103) and (7.106), one can write the radiation power as

dP = _EbdA1 (7.107)

and since radiation power and flux are related by

dP = q dA1 (7.108)

the flux satisfies

q = _(θ)Eb = _(θ)σθ4. (7.109)

Comparing Equation (7.109) with Equation (7.24) for θe = 0 (no irradiation) reveals that

A(θ) = _(θ)σ. (7.110)

In reality, we not only have radiation leaving the body but also irradiation entering the

body. The spectral, directional irradiation intensity II in a certain direction (ϕ,ψ) is defined

as the irradiation power per unit solid angle dω about this direction, per unit wavelength

dλ, per unit receiving area perpendicular to this direction:

II = dPI

dω dλ dA1 cosψ

. (7.111)

Likewise, the total hemispherical irradiation power G is defined as the irradiation power

per unit receiving area:

G = _

∞

0

_

hemisphere

II cosψ dω dλ. (7.112)

A part of the irradiation power is absorbed (αG), a part of it is reflected (ρG) and a part

of it is transmitted (τG) (Figure 7.4).

ρG

αG

τG

G

Figure 7.4 Absorption, reflection and transmission of irradiation

HEAT TRANSFER 321

Energy conservation requires that α + ρ + τ = 1. We assume that we are dealing with

opaque materials, that is, materials for which there is no transmission. Accordingly, τ = 0

and τ = 0. Accordingly,

α + ρ = 1. (7.113)

α is the total hemispherical absorptivity and ρ is the total hemispherical reflectivity. In

reality, α and ρ are dependent on the irradiation angle and its spectrum. Therefore, α and

ρ are averaged values in the same sense as _ is an averaged value of _λ,ω.

Here and in the section that follows, we assume that we deal with

1. diffuse surfaces, that is, _ and α are independent of the radiation and irradiation

direction;

2. gray surfaces, that is, _ and α are independent of the wavelength for the actual range

of interest.

Under these conditions, the important relationship

α = _ (7.114)

applies (Incropera and DeWitt 2002), that is, the absorptivity equals the emissivity. Looking

at Figure 7.4, the total radiation leaving the surface is the sum of the total emissive power

E and the reflected total irradiation power. This is called the total radiosity J :

J = E + ρG. (7.115)

Now we arrive at the view-factor concept. The view factor Fij is defined as the fraction

of the radiation power leaving surface i that is intercepted by surface j. It is assumed that

the surface Ai is characterized by a uniform radiosity Ji . The total radiation leaving the

surface Ai amounts to

R = JiAi . (7.116)

Since the radiosity is assumed to be uniform, the directional radiosity Jω,i satisfies

Jω,i = Ji

π

(7.117)

and the radiosity leaving dAi and reaching surface dAj yields (Figure 7.5)

dRij = Jω,i dAi cosψiωij (7.118)

where ωij is the view angle covered by dAj seen by dAi :

ωij = dAj cosψj

R2 . (7.119)

Accordingly,

Fij = 1

AiJi

_

Ai

_

Aj

 Jω,i cosψi cosψj

R2

!dAi dAj

= 1

Ai

_

Ai

_

Aj

 cosψi cosψj

πR2

!dAi dAj . (7.120)

322 HEAT TRANSFER

x

R

y

z

ϕ

dAi

dAj

ψi

ψj

Figure 7.5 Geometry for the view-factor calculation

Important relations are the reciprocity relation

AiFij = AjFji (7.121)

and the summation rule for enclosures

N

_

j=1

Fij = 1. (7.122)

Now consider N surfaces Ai interacting with each other. From Figure 7.4, we obtain the

relationships

qi = Ei − αiGi = Ei − _iGi (7.123)

Ji = Ei + ρiGi = Ei + (1 − _i)Gi . (7.124)

Hence, eliminating Gi from Equations (7.123) and (7.124),

qi = Ei − _iJi

1 − _i

= _i(Ebi − Ji )

1 − _i

(7.125)

where Ebi stands for the blackbody radiation of surface i. Eliminating Ei from Equations

(7.123) and (7.124) leads to

qi = Ji − Gi . (7.126)

HEAT TRANSFER 323

Gi is the irradiation from all other bodies. Conservation of energy requires

AiGi =

N

_

j=1

j =i

FjiAj Jj . (7.127)

Using the reciprocity rule, Equation (7.127) can be rewritten as

Gi =

N

_

j=1

j =i

Fij Jj . (7.128)

Accordingly, Equation (7.126) now reads

qi = Ji −

N

_

j=1

j =i

Fij Jj . (7.129)

Equating Equations (7.125) and (7.129) yields

Ji − (1 − _i)

N

_

j=1

j =i

Fij Jj = _iEbi . (7.130)

If the temperatures of all the participating surfaces are known, Equation (7.130) constitutes

a set of N linear equations in the N unknowns Ji . This set is not necessarily symmetric.

After solving for Ji , the fluxes qi can be obtained through Equations (7.125) or (7.129).

From qi , an equivalent environmental temperature can be derived for each surface Ai using

Equation (7.24):

θei = _

θ4

i

− qi

Ai(θi)

_1/4

(7.131)

where θi is the mean temperature of surface i. Sometimes a cavity is not completely closed

and part of the radiation escapes to the environment. Considering this environment to

behave as a blackbody and attributing it to the surface k, one obtains

_k = 1 ⇒ Jk = Ebk = Eb,environment (7.132)

and Equation (7.130) now yields

Ji − (1 − _i)

N

_

j=1

j =i,k

Fij Jj − (1 − _i)FikEbk = _iEbi (7.133)

or, since

Fik = 1 −

N

_

j=1

j =i,k

Fij (7.134)

324 HEAT TRANSFER

ϕ

ψ

p1i

p2i

p3i

e1i

e2i

ni

nj

ψij

ϕij

R ψj i ij

Triangle i

Triangle j

ci

cj

Figure 7.6 Local coordinate system in triangle i

one obtains

Ji − (1 − _i )

N

_

j=1

j =i,k

Fij Jj = _iEbi + (1 − _i)(1 −

N

_

j=1

j =i,k

Fij)Eb,environment. (7.135)

7.7.2 Numerical aspects

The time-consuming part in generating Equation (7.130) is the calculation of the view

factors. The method proposed here consists of the following steps:

1. Triangulate the free surface of the structure by defining linear triangles within the

element faces without generating any new nodes. For instance, a face of a 20-node

brick element is divided in six triangles, a face of an 8-node brick element in two

triangles and a face of a 10-node tetrahedral element in four triangles. Number the

nodes within each triangle in counterclockwise direction when viewed from outside

the body.

2. For each triangle i, determine the following:

(a) The center of gravity ci .

(b) The normal ni , pointing away from the body (Figure 7.6):

ni = (p2i

− p1i) × (p3i

− p2i )

_(p2i

− p1i) × (p3i

− p2i )_. (7.136)

(c) The area Ai .

HEAT TRANSFER 325

(d) The unit vector e1i satisfying

e1i

= (p2i

− p1i )

_p2i

− p1i

_ (7.137)

(e) The unit vector e2i

= ni × e1i . The basis (e1i , e2i , n) defines a right-handed

rectangular coordinate system.

(f) The scalar

di = −p1i

・ ni . (7.138)

A point p lies in the plane of triangle i if

p ・ ni + di = 0. (7.139)

It is visible from triangle i if and only if (assuming no other triangles block the

view)

p ・ ni + di ≥ 0. (7.140)

3. For each triangle i:

(a) Perform a loop over all triangles j = i with the following actions:

(i) Check whether cj is visible from triangle i. If it is not, that is, if

cj ・ ni + di < 0 (7.141)

cycle

(ii) Check whether ci is visible from triangle j. If it is not, that is, if

ci ・ nj + dj < 0 (7.142)

cycle. Only those triangles j remain from which triangle i can be seen and

which are visible from triangle i (assuming no other triangles block the

view). In the remainder of the text, they will be called visible triangles.

(iii) Calculate the distance

Rij = _cj − ci_ (7.143)

and the unit vector

ξ ij

= (cj − ci)/Rij . (7.144)

326 HEAT TRANSFER

ϕ

ψ

ϕmin ϕmid ϕmax

_ψ

_ϕ

Figure 7.7 φ − ψ grid

(b) Generate a rectangular grid with ϕ on the x-axis and ψ on the y-axis. A

(ϕ,ψ) pair uniquely defines a direction in the local (e1i , e2i , n) system, where

0 < ϕ < 2π, 0 < ψ < π/2 (cf Figure 7.6). The (ϕ,ψ)-range is meshed with

an N ×M rectangular grid (Figure 7.7). Let k and l be functions such that

k(ϕ) and l(ψ) denote the discrete grid element to which (ϕ,ψ) belongs. If

_ϕ = 2π/N and _ψ = π/(2M), then the functions satisfy

k(ϕ) = int(ϕ/_ϕ) + 1 (7.145)

l(ψ) = int(ψ/_ψ) + 1 (7.146)

where int(x) is the largest integer smaller than or equal to x. Initialize by

assuming that all grid elements are uncovered.

(c) For all visible triangles, order Rij in ascending order.

(d) Perform a loop over all visible triangles j = i in ascending Rij -order with the

following actions:

(i) Calculate the coordinates of ξ ij in the local (e1i , e2i , n) system and determine

the angles ϕij and ψij by inverting the relations

ξ ij

・ e1i

= sinψij cos ϕij (7.147)

ξ ij

・ e2i

= sinψij sin ϕij (7.148)

ξ ij

・ ni = cosψij . (7.149)

HEAT TRANSFER 327

(ii) Determine the grid element k(ϕij ), l(ψij ) and check whether it was already

covered. If so, cycle.

(iii) Calculate the view factor

Fij = cosψij cosψjiAj

πR2

ij

(7.150)

where

cosψji = −ξ ij

・ nj . (7.151)

(iv) Determine which grid elements are covered by triangle j. To that end,

calculate the unit vectors qkij connecting ci with pkj , k = 1, 2, 3:

qkij

=

(pkj

− ci )

_pkj

− ci_. (7.152)

If

n31ij

= q3ij

× q1ij (7.153)

n12ij

= q1ij

× q2ij (7.154)

n23ij

= q2ij

× q3ij (7.155)

then the equations of the planes connecting the edges of triangle j with ci

satisfy

p ・ n31ij

= 0 (7.156)

p ・ n12ij

= 0 (7.157)

p ・ n23ij

= 0, (7.158)

(see Figure 7.8). A unit vector p with coordinates (sinψ cos ϕ, sinψ sin ϕ,

cos ψ) lies in the plane defined by p3j , p1j and ci if

ψ = −tan−1 _ (n31ij

・ ni )

(n31ij

・ e1i ) cos ϕ + (n31ij

・ e2i ) sin ϕ

_ =: f3/1(ϕ)

(7.159)

and likewise for the other planes. The spatial angle covered by triangle j

corresponds to a triangle with curved sides in the ϕ,ψ plane (Figure 7.7).

Its vertices are made up of "ϕ(p1j ),ψi(p1j )#, "ϕ(p2j ),ψi(p2j )# and

"ϕ(p3j ),ψi(p3j )#. Now determine ϕimin and ϕimax defined by

ϕmin = min{ϕ(p1j ), ϕ(p2j ), ϕ(p3j )} (7.160)

ϕmax = max{ϕ(p1j ), ϕ(p2j ), ϕ(p3j )}. (7.161)

328 HEAT TRANSFER

Triangle j

ci

p1j

p2j

p3j

q1ij

q2ij

q3ij

Figure 7.8 Spatial angle covered by triangle j

The one that is left is called ϕmid. Let ψmin be the ψ-value corresponding

to ϕmin and ψmax the ψ-value corresponding to ϕmax. First, the grid

elements (k(ϕmin), l(ψmin)) and (k(ϕmax), l(ψmax)) are marked as covered.

Then, for k(ϕmin) ≤ m ≤ k(ϕmid) and ϕmin ≤ (m − 12

)_ϕ ≤ ϕmax, those

grid elements are marked as covered for which

min $l{fmin/max[(m − 12

)_ϕ]}, l{fmin/mid[(m − 12

)_ϕ]}% ≤ n ≤

max $l{fmin/max[(m − 12

)_ϕ]}, l{fmin/mid[(m − 12

)_ϕ]}% (7.162)

and for k(ϕmid) ≤ m ≤ k(ϕmax) and ϕmin ≤ (m − 12

)_ϕ ≤ ϕmax the elements

for which

min $l{fmin/max[(m − 12

)_ϕ]}, l{fmid/max[(m − 12

)_ϕ]}% ≤ n ≤

max $l{fmin/max[(m − 12

)_ϕ]}, l{fmid/max[(m − 12

)_ϕ]}% . (7.163)

This corresponds to the crossed elements in Figure 7.8.

(v) If less than _% of the grid elements are marked as uncovered, exit the loop.

4. Solve Equations (7.135) to obtain Ji , or, equivalently, qi . The corresponding temperatures

θei (Equation (7.131)), can be used as the thermal boundary condition in

the next iteration. Since the structure deforms, the view factor has to be recalculated

in each iteration. However, the changes are usually small and the algorithm can be

accelerated by making diligent use of the visible triangles from the last iteration.

Навигация