10 HEAT TRANSFER BY RADIATION
10.1 Radiative Transfer through Transparent Media
In this Section, the information on radiation fundamentals and properties developed in Chapter 9 is used to find the net rate of energy transfer to/from a material by radiation. Because of the complexities possible in a complete treatment of the spectral and directional characteristics of radiation, it is useful to first develop some simple if less accurate methods for approximating the radiative transfer. If radiation provides only 10 percent of the total heat transfer by all modes, then it is perfectly acceptable to determine the radiative transfer to within 20 percent, as any error is only ± 2 percent of the total heat transfer. If radiation is 90 percent of the total heat transfer (often the case in radiant heaters and coal-fired steam generators,) then a much more accurate approach is indicated.
10.1.1 Transfer between Two Areas
Consider first a surface of area A, temperature T, and total emissivity ε. The rate of radiative energy emission from this surface per unit area is, from Eq. (9.4), E = εσ T 4 . Let the surroundings in all directions be at temperature T∞. If the surroundings are far from surface A or are black, then little of the energy emitted from A will be reflected from the surroundings. The rate of energy incident on A per unit area from the surroundings is denoted by the symbol G, and is given by 4 G = σ T∞ (10.1) The net radiative heat flux on surface A is the difference between the emitted and absorbed radiative fluxes, or 4 q′′ = E − α G = εσ T 4 − ασ T∞ (10.2) This can be simplified if some further approximations or assumptions are made. First, if the surface is gray, then Kirchhoff's Law gives α = ε . If T∞ ≈ T, comparing Eqs. (9.17) and (9.27) also indicates that α ≈ ε . In either of these two cases, eq. (10.2) reduces to 4 q′′ = εσ ( T 4 − T∞ ) (10.3)
794 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
This is an extremely useful relation, but it involves some implicit assumptions in addition to those explicitly mentioned above. For example, surface A and the surroundings are assumed to each be isothermal, and the value of G does not vary across A so that the local heat flux at every point on A is the same. Equation (10.3) is often used as part of a boundary condition for use in mixed mode problems, and it is clearly a simple way for handling radiation when the conditions and assumptions are satisfied. A further simplification of eq. (10.3) can be made by invoking the same approximation used in eq. (9.60), (T 4 − T∞4 ) = (T 2 + T∞2 )(T 2 − T∞2 ) = (T 2 + T∞2 )(T + T∞ )(T − T∞ ) (10.4) If the absolute temperatures T and T∞ are not too far apart, then each can be replaced in the first two parentheses by their average, T = (T + T∞ ) / 2 , resulting in
(T
4
= 2T 2 × 2T ( T − T∞ ) = 4T 3 ( T − T∞ )
4 2 − T∞ ) = ( T 2 + T∞ )( T + T∞ )( T − T∞ )
(10.5)
and eq. (10.3) becomes
q′′ = 4εσ T 3 ( T − T∞ ) = hrad ( T − T∞ )
3
(10.6)
where hrad = 4εσ T is the radiation heat transfer coefficient. The linearized form of eq. (10.6) makes possible the simple addition of radiation and convection terms in some combined-mode problems, if the loss in accuracy caused by the many assumptions and simplifications is justified. Note also that both the surface and surrounding temperatures must be known in order to determine T ; otherwise, an iterative solution must be invoked. Example 10.1 Let T∞= 27oC, and vary T from 300 to 1000K. Determine the error in q" by comparing the heat flux predicted by eq. (10.6) with that predicted by eq. (10.3).
Figure 10.1 Comparison of heat fluxes
Chapter 10 Heat Transfer by Radiation 795
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Solution: The resulting comparison is shown in the Fig. 10.1. Note that the result is independent of the value of emissivity that is chosen. The error is less than 6 percent up to surface temperatures of 500 K, which is within the range of many combined mode problems.
10.1.2 Diffuse Surfaces: The Configuration Factor
In more complex cases than treated in Section 10.1.1, a surface is exchanging radiative energy with multiple surfaces, each with a different temperature and different properties. It is instructive to first examine the radiative exchange between a single pair of these surfaces. These are shown in Fig. 10.2. A line between the two elements has length S1-2, and the angles between this connecting line and the normals to the elements are θ1 and θ2. Both surface elements are assumed to be black. In the following relations, subscripts such as A1 or dA1 are replaced by 1 or d1 to avoid two levels of subscripts. Equation (10.7) gives the radiation leaving element 1 that is incident on element 2, assuming there is no attenuation caused by a medium separating the elements, as dA cos θ deλ ,d1− d 2 = I λb ,1dA1 cos θ1dΩ 2 dλ = I λb,1dA1 cos θ1 2 2 2 dλ (10.7)
S1− 2
Because the element is black, the intensity is independent of θ. The solid angle subtended by surface 2, dΩ2, has been replaced by the projected area of element 2 over S1-22. Also, the rate of energy incident on element 2 is also the rate of energy absorbed by element 2 because it is black.
θ2 θ1
S1-2
Surface element dA1 at T1; ε1 = 1
Figure 10.2 Radiative exchange between two differential surface elements.
796 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Similar reasoning gives the rate of energy emitted by element 2 that is absorbed by element 1 as dA cos θ deλ ,d 2 − d1 = I λb ,2 dA2 cos θ 2 1 2 1 d λ (10.8)
S1− 2
The net energy exchange between the two elements is then
d 2 qλ ,d1− d 2 = deλ ,d1− d 2 − deλ , d 2 − d1 dA cos θ dA cos θ = I λb ,1dA1 cos θ1 2 2 2 dλ − I λb ,2 dA2 cos θ 2 1 2 1 d λ S1− 2 S1− 2 cos θ1 cos θ 2 dA2 = ( I λ b,1 − I λb ,2 ) dA1 dλ S12− 2 cos θ1 cos θ 2 dA2 = ( Eλb ,1 − Eλ b,2 ) dA1 π S12− 2 dλ
(10.9)
The final form uses the identity Eλb = πIλb. The net radiative flux on element 1 due to exchange between surfaces 1 and 2 is then d 2 qλ ,d1− d 2 cos θ1 cos θ 2 dA2 d 2 q "λ ,d1− d 2 = = ( E λb,1 − E λb,2 ) d λ (10.10) dA1 π S12− 2 Now, observe that the effects of geometry are all contained in a wavelength- and temperature-independent term defined as cos θ1 cos θ 2 dA2 dFd1− d 2 ≡ (10.11) π S12− 2 so that eq. (10.10) becomes d 2 q "λ ,d1− d 2 = ( Eλb,1 − Eλb,2 ) dFd1− d 2 dλ (10.12) The dFd1− d 2 is variously called the configuration factor, shape factor or angle factor. The configuration factor is interpreted as the fraction of radiation leaving surface 1 that is incident on surface 2. Equation (10.11) is the basic definition of the configuration factor; it will now be extended to show some important properties. It will then be applied to configuration factors for areas of finite size. The derivation above also applies to diffuse nonblack surfaces, so that the use of the configuration factor is valid for radiative exchange between black or diffuse surfaces as well as gray or spectrally dependent surfaces because the factor is independent of wavelength. Reciprocity. The rate of energy transfer from element 1 to element 2 must be the negative of the energy transfer rate from element 2 to element 1, or qλ ,d1− d 2 = ( Eλb ,1 − E λb,2 ) dFd1− d 2 dA1d λ = −qλ , d 2 − d1 (10.13) = − ( E λb,2 − Eλ b,1 ) dFd 2 − d1dA2 dλ which results in the identity
dFd1− d 2 dA1 = dFd 2 − d1dA2
(10.14)
Chapter 10 Heat Transfer by Radiation 797
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
This reciprocity relation is valid without restriction for diffuse surfaces of differential area. Exchange between a differential area and a finite area. If surface 2 is finite in extent, (Figure 10.3), then eq. (10.7) becomes cos θ1 cos θ 2 deλ ,d1− 2 = π I λb ,1dλ dA1 (10.15) dA2 = Eλ b,1d λ dA1 A dFd1− d 2 2 A π S1− 2 Defining cos θ1 cos θ 2 Fd1− 2 = (10.16) dA2 = A dFd1− d 2 2 A π S1− 2 eq. (10.15) becomes deλ ,d1− 2 = E λb ,1dλ dA1 Fd1− 2 (10.17)
2 2 2 2
Observe that in the argument of the integral in eq. (10.16), cosθ1, cosθ2 and S1-2 will all vary as the integration across finite area A2 is carried out. We might infer some difficulty in evaluating Fd1− 2 when geometries are complex. Fortunately, there are resources that allow some simplifications. If the intensity leaving surface 2 is uniform across that surface, then the energy transfer from surface 2 to surface 1 is the emitted energy rate times the fraction leaving surface 2 that reaches surface 1, or deλ , 2 − d1 = E λb ,2 A2 d λ dF2 − d1 (10.18) Reciprocity for element-area factors: If black element dA1 and black area A2 are at the same temperature, then the rate of energy exchange between them must be equal. Equating eqs. (10.17) and (10.18) gives a second reciprocity relation dA1 Fd1− 2 = A2 dF2 − d1 (10.19) Surface element dA2 at T2; ε2 = 1
θ2
θ1 S1-2 Surface area A2 at T2; ε2 = 1
Surface element dA1 at T1; ε1 = 1
Figure 10.3 Exchange between surface element dA1 and finite surface A2.
798 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Exchange between two finite areas. If both surface 1 and 2 are finite in extent, then eq. (10.18) can be extended to give deλ ,2 −1 = Eλ b,2 A2 dλ dF2 − d1 = E λb,2 A2 d λ F2 −1 (10.20) A
1
Equation (10.15) gives
cos θ1 cos θ 2 deλ ,d1− 2 = E λb ,1dλ dA2 dA1 = Eλb ,1dλ A A dFd1− d 2 dA1 2 A1 A2 1 2 π S1− 2
(10.21) Substituting eq. (10.19) into eq. (10.21) results in
deλ ,d1− 2 = E λb ,1dλ dFd1− d 2 dA1 = E λb,1dλ Fd1− 2 dA1 = E λb,1 A1 F1− 2 dλ A1 A2 A1
(10.22) Equation (10.22) also indicates that a general definition of F1-2 is cos θ1 cos θ 2 dA2
F1− 2 = Fd1− 2 dA1 =
A1 A1
A2
dFd1− d 2 dA1 = A1 A 2
S12− 2
dA1
(10.23) If T1 = T2, then equating eqs. (10.20) and (10.22) gives a final reciprocity relation A1 F1− 2 = A2 F2 −1 (10.24) Summation relations: As an aid to further use of configuration factors, consider surface j that is part of an enclosure that completely surrounds surface Aj (Figure 10.4). The enclosure is made up of a total of N surfaces. Based on the concept that Fj-k is the fraction of radiation leaving j that is incident on k, it is clear that
F
k =1
N
j −k
=1
(10.25)
A2
Fj-2
Aj
Fj-1 A1 Fj-N Fj-k
AN
Ak
Figure 10.4 Configuration factors in an enclosure of N surfaces
Chapter 10 Heat Transfer by Radiation 799
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
If surface j is concave, then it is possible that some radiation leaving j may strike another portion of j, so that Fj-j ≠ 0 and must be included in the summation. A corollary to eq. (10.24) is that the fraction of energy leaving j that is incident on two surfaces k and l must equal to the fractions incident separately on the two surfaces, or F j − ( k + l ) = F j − k + F j −l (10.26)
10.1.3 Configuration Factor Algebra
The reciprocity equations [(10.14), (10.19) or (10.24)], the summation relation, eq. (10.25), and the associative relation eq. (10.26) are the basis of configuration factor algebra. Configuration factors for simple geometries can be found directly by applying configuration factor algebra, eliminating the integrations that may otherwise be necessary. The method is shown in the following two examples. Example 10.2 Find the configuration factors F1-2, F2-1, and F2-2 for the case of infinitely long concentric circular cylinders (Fig. 10.5). The inner cylinder has outer diameter D1 and the outer cylinder has inside diameter D2. Solution: By definition, all radiation leaving surface 1 is incident on surface 2, so F1-2 = 1. A π D1 L D1 Using reciprocity, F2 −1 = 1 F1− 2 = = . Finally, using the π D2 L D2 A2 summation relation, concave, so F2-2
F
k =1
N
2−k
= F2 −1 + F2 − 2 = 1 . In this case, surface 2 is
must
be
included
in
the
summation.
Thus, F2 − 2 = 1 − F2 −1 = 1 −
D1 . The full set of configuration factors was D2
found for this geometry without recourse to integration.
D2
D1
Figure 10.5 Infinitely long concentric cylinders
800 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Example 10.3 An infinitely long enclosure has a cross-section made up of a scalene triangle with planar sides L1, L2, and L3 (Fig. 10.6). Find an expression for the configuration factor F1-2 in terms of the side lengths. Solution: eq. (10.25) can be written for each surface. In this case, for an enclosure of all planar surfaces, each term Fj-j = 0, so the summation rule, eq. (10.25), gives
F1− 2 + F1− 3 = 1 F2 −1 + F2 − 3 = 1 F3−1 + F3− 2 = 1
This results in three equations involving six unknowns: not a good portent! However, three of the unknowns may be eliminated by using reciprocity. Multiply each equation by its respective area (assuming an equal length of channel), resulting in
L1 F1− 2 + L1 F1− 3 = L1 L2 F2 −1 + L2 F2 − 3 = L2 L3 F3−1 + L3 F3− 2 = L3
Using reciprocity to eliminate three of the unknowns gives
L1 F1− 2 + L1 F1− 3 = L1 L1 F1− 2 + L2 F2 − 3 = L2 L1 F1− 3 + L2 F2 − 3 = L3
Solving this set results in F1− 2 =
L1 + L2 − L3 . 2 L1
L2 L3
L1
Figure 10.6 Infinitely long enclosure with three planar sides
Once the factor for the triangular enclosure is known, factors for other geometries can be derived. For example, if L1 = L2, then the factor between two sides of a wedge is given by
F1− 2 = 2 L1 − L3 L = 1 − 3 = 1 − sin α 2 L1 2 L1
(10.27)
where α is the half-angle between the sides 1 and 2. Clearly, a number of other factors can also be generated. However, all factors can be generated from
Chapter 10 Heat Transfer by Radiation 801
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
configuration factor algebra only for two-sided enclosures if one of the sides is concave as in Example 10.2, or for a three-sided enclosure if all sides are planar as in Example 10.3. There are exceptions to this general rule if symmetry allows some of the factors to be equal. In all cases, configuration factor algebra allows reduction in the number of factors that must be determined independently. The number of independent factors U that must be found in order to determine all factors in an N-sided enclosure made up of planar or convex surfaces is (Siegel and Howell, 2002): U = N ( N − 3) / 2 (10.28) and if M of the N surfaces are concave, then U = [ N ( N − 3) / 2] + M (10.29) For the case of concentric cylinders in Example 10.2, N = 2 and M = 1, so from eq. (10.29), U = [ 2(2 − 3) / 2] + 1 = −1 + 1 = 0 and, as found in the example, no independent factors must be found; all are determined by configuration factor algebra. For the triangular enclosure of planar surfaces in Example 10.3, N = 3 and M = 0, so U = 3 ( 3 − 3) / 2 = 0 , and again no independent factors must be found. The crossed-string method. We now derive an extremely useful extension of the factor for the triangular enclosure derived in Example 10.3 (Hottel, 1954). The restriction still applies that the result will be for an enclosure that is infinitely long in one dimension. Observe an enclosure cross-section made up of four surfaces shown in Fig. 10.7(a). We seek F1-2. Connect the corners a and b, and corners b and c with lines that follow around any convex areas as shown in Fig. 10.7(b). Do the same with lines connecting points a and d, and b and c. Then connect points a and c, forming triangles (a-b-c) and (a-c-d). These triangles are made up of all convex surfaces composed of the constructed lines. Similarly, connect points b and d with a dotted line, making a second pair of triangles (b-c-d) and (a-b-d) also made up of
A1 A1 a b
A3
A4
A3
A4
A2 A2 a) Enclosure d b) Enclosure divided into triangular portions c
Figure 10.7 Enclosure of four parallel infinitely long surfaces
802 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
convex or planar surfaces. Now, applying the result of Example 10.3 for a triangle of planar or convex surfaces to the two constructed triangles (a-b-d) and (I) in Fig. 10.7(b):
Fcd − bc = Fcd − ad = Lcd + Lbc − Lbd 2 Lcd Lcd + Lad − Lac 2 Lcd
(10.30)
where the L values are the lengths of the constructed lines connecting the end points. Equation (10.25) applied to this geometry gives Fcd − ad + Fcd − ab + Fcd − bc = 1 (10.31) Noting that all radiation reaching the constructed boundary ab will also be incident on L1, eq. (10.31) becomes Fcd − ab = Fcd −1 = 1 − ( Fcd − ad + Fcd − bc ) (10.32) Applying reciprocity gives
F1− cd = cd cd Fcd −1 = 1 − ( Fcd − ad + Fcd − bc ) L1 L1
(10.33)
Again, observe that F1− cd = F1− 2 because all radiation incident on cd will finally be incident on surface 2. Substituting eq. (10.30) into (10.33) gives the final result
F1− 2 = = Lcd L L + Lad − Lac Lcd + Lbc − Lbd + Fcd −1 = cd 1 − cd L1 L1 2 Lcd 2 Lcd ( Lbd + Lac ) − ( Lad + Lbc ) 2 L1
(10.34)
The final form of the equation can be interpreted as follows; consider the values of L as the length of taut strings that connect the corners a, b, c, d. The final form shows that the factor F1-2 is equal to the sum of the lengths of the strings connecting the corners of surface 1 and 2 that cross one-another, minus the sum of the lengths of the uncrossed strings that connect the corners, all divided by the length of surface 1. Example 10.4 Find the configuration factor between two infinitely long directly opposed plates as shown Fig. 10.8(a).
W 2 H 1 (a) Figure 10.8 Two infinitely long directly opposed plates (b) H W
Chapter 10 Heat Transfer by Radiation 803
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Solution: Connect the corners with the solid (uncrossed) and dotted (crossed) lines shown in Fig. 10.8(b). Equation (10.34) then gives
F1− 2 =
( crossed ) − ( uncrossed )
2 L1
1/ 2
=
2( H 2 + W 2 )1/ 2 − 2 H 2L
H 2 = + 1 W
H − W
Two methods have now been examined for finding values of configuration factors in terms of the particular variables that describe a given geometry: Configuration factor algebra, and the crossed-string method (remembering that the latter method is restricted to geometries that are infinitely long in one dimension.) For more general geometries, there are two additional possibilities for determining the configuration factors that are necessary in calculating radiative energy transfer. First, configuration factors for many geometries have been derived, and the results published in the open literature. These are collected in a web-based catalog (Howell, 2003) that includes well over 300 geometries with element-element, element-area, and area-area factors. Some of these are given in the Appendix in Table F.1. Second, methods for carrying out the required integrations in the general definitions of eqs. (10.16) or (10.23) can be implemented. These methods include direct numerical integration using commercial codes, contour integration (which replaces integrations over the receiving area with integrations around the area boundary), and differentiation of known finite-area factors to find element-area factors. Detailed information on these methods is in standard radiative transfer texts (Siegel and Howell, 2002; Modest, 2003).
10.2 The Enclosure; The Net Radiation Method for Diffuse Surfaces
Consider an enclosure as shown in Fig. 10.9. Such an enclosure is useful because it allows a mathematical description that makes use of the configuration factor summation rules, resulting in a set of final equations that can be solved for the required unknowns. The enclosure must be constructed so that any surface k views other enclosure surfaces regardless of the direction from k. This would seem restrictive for many problems; however, a virtual enclosure can always be constructed by replacing openings in a physical enclosure with a virtual surface that has the radiative behavior of the opening.
804 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
A2
A1
T∞
A3
Figure 10.9 Open-ended cylindrical enclosure closed with virtual surface A3.
Consider the open-ended cylindrical enclosure shown in Figure 10.9. To make the system into a radiative enclosure, the opening at the left end is replaced with a virtual black surface A3 at T∞. This virtual surface has the same radiative behavior as the opening, because no radiation from the internal cylinder surfaces will reflect from the opening; the opening acts as a perfect absorber (i.e., it is black), and radiation entering the enclosure from the surroundings will be at T∞. Thus, the opening can be replaced with a virtual black surface at T3 = T∞. This strategy can almost always be used to construct a virtual enclosure from an open radiating system.
10.2.1 Radiosity, Irradiation, and Net Energy Transfer
An individual surface k within the enclosure will have energy incident on it from all other surfaces in the enclosure. Let this incident radiation per unit area Ak in the wavelength interval dλ from all surfaces in the enclosure (possibly including from k itself if it is concave) be given the symbol Gλ,kdλ. The G is called the irradiation on surface k. A portion of the irradiation striking surface k will be reflected from k. Let k be a diffuse surface. In that case, the reflected radiative flux from Ak will be (10.35) dq "λ ,ref ,k = ρλ ,k Gλ , k dλ The spectral radiative energy flux leaving Ak will be made up of this reflected flux plus the radiation emitted by the surface. The emitted plus reflected radiation leaving the surface is called the radiosity of the surface, and is given the symbol J, so that J λ ,k = emitted + reflected flux=ε λ ,k E λb, k + ρλ ,k Gλ , k (10.36)
Chapter 10 Heat Transfer by Radiation 805
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Gλ,jdλ
Jλ,kdλ Tk, ελ,k, Ak
dq"λ,k
Figure 10.10 Radiative balance on surface Ak
Because the surfaces are assumed to be diffuse, both the emitted and reflected flux are independent of direction, and the use of configuration factors is justified in finding energy exchange among surfaces. We are now in a position to perform a radiative energy balance on surface k. The factors we can use are shown in Fig. 10.10. The energy balance is then ′′ dqλ ,k = (outgoing-incoming radiation) = J λ ,k dλ − Gλ , k d λ (10.37)
′′ The dqλ ,k is the heat flux in wavelength interval dλ that must be provided to surface k in order to balance the difference between the outgoing and incoming radiation fluxes. Equations (10.36) and (10.37) provide two equations involving the factors J, G, q" and Eb. It is required that either Eb or q" (effectively temperature or radiative heat flux) be specified as a boundary condition. Nevertheless, we are left with two equations in three unknowns for each surface k, so an additional relation is required. Here, we can invoke the concept of an enclosure to provide the additional relation. The factor Gλ,k is made up of the radiation incident on Ak from all N surfaces making up the enclosure. The radiation leaving each of these surfaces is made up of the emitted plus reflected radiation (the radiosity.) Thus, the irradiation on Ak is Gλ ,k Ak = J λ , j A j F j − k
j =1 N
(10.38)
or, invoking reciprocity from eq. (10.23),
Gλ ,k = J λ , j Fk − j
j =1 N
(10.39)
806 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
The three equations (10.36), (10.37) and (10.39) are the basis of the net radiation method for treating radiative transfer. Note the restrictions that are implicit in the method so far: all surfaces making up the enclosure are diffuse; otherwise the configuration factors could not be used in eqs. (10.38) and (10.39). Further, the values of Gλ,j and Eλb,j must be invariant across each surface j, so that Jλ,j will be uniform and the configuration factors can be used. This condition is difficult to meet exactly in enclosures with large surfaces, where in particular the irradiance G will probably vary across a large surface. In that case, surfaces must be subdivided until the requirement of uniform radiosity from each surface is met. An alternative approach to radiative exchange in enclosures is the radiative network method (Oppenheim, 1956). This approach is often used in undergraduate texts (e.g., Incropera et al., 2007). It is based on the same fundamental assumptions and relations as in the net radiation approach used here, but is somewhat more difficult to extend to more than a few surfaces and to cases with wavelength dependence.
10.2.2 Gray Surfaces
Because wavelength-dependent properties are often not available, the net radiation method is quite often applied under the further assumption that all surfaces in the enclosure are gray. In that case, the three basic equations can be integrated over wavelength to become J k = ε k Eb, k + ρ k Gk = ε kσ Tk4 + (1 − ε k ) Gk (10.40) ′′ qk = J k − Gk (10.41)
Gk = J j Fk − j
j =1 N
(10.42)
The equations can be written for each of the N surfaces in an enclosure, resulting in 3N equations in 3N unknowns (Jk, Gk and either q"k or Tk for each surface). Before proceeding to some examples, the method can be simplified still further. In most radiative transfer problems, the explicit value of G is not needed, so G can be algebraically eliminated from the equation set. The Gk can be found from eq. (10.41) and substituted into eqs. (10.40) and (10.42), giving ′′ (10.43a) J k = ε kσ Tk4 + (1 − ε k )( J k − qk ) or (1 − ε k ) ′′ J k = σ Tk4 − qk (10.43b) εk and
′′ J k = qk + J j Fk − j
j =1 N
(10.44)
Using eqs. (10.43b) and (10.44) for each of the N surfaces reduces the set to 2N equations in 2N unknowns. The radiosity is sometimes needed in an engineering
Chapter 10 Heat Transfer by Radiation 807
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
solution, because it is the quantity that is observed by radiometers or radiation pyrometers viewing a surface in an enclosure. If the radiosity is not required, then it can also be eliminated from the equation set by substituting eq. (10.43b) into eq. (10.44) to give N (1 − ε j ) q′′ F 1 ′′ k− j σ Tk4 − qk = σ T j4 − (10.45) j εk εj j =1 To place the temperatures on one side of the equation, eq. (10.45) is rearranged to the form N N δ (1 − ε j ) kj (10.46) Fk − jσ (Tk4 − T j4 ) = ε − Fk − j ε q′′j j =1 j =1 j j where δ kj = 1 if k = j and δ kj = 0 otherwise. Equation (10.46) is written for each surface k in an enclosure, resulting in N equations in N unknowns. One consequence of the form of eq. (10.46) is that the temperature of at least one surface in the enclosure must be specified; if only q" is specified for each surface, the solution for Tk becomes indeterminate. For pure radiative transfer (no conduction or convection present), eq. (10.46) results in a set of linear equations with unknowns T4 or q" for each surface; any standard method for solving a set of linear equations can be used to determine the unknowns. The matrix of coefficients that results from the formulation may be fully populated if every surface in the enclosure can exchange energy with every other element. Example 10.5 Derive an expression for the radiative heat flux between two infinite parallel plates that are at (ε1, T1) and (ε2, T2) using the net radiation method. Solution: For this case, eq. (10.46) has N=2,
F
j =1
2
k− j
σ ( Tk4 − T j4 ) =
2
δ kj
j =1 ε j
− F1− j
(1 − ε ) q′′
j
εj
j
Expanding this for k = 1 gives F1−1σ ( T14 − T14 ) + F1− 2σ ( T14 − T24 )
δ (1 − ε1 ) ′′ δ12 (1 − ε 2 ) ′′ = 11 − F1−1 − F1− 2 q + q ε1 ε1 1 ε 2 ε2 2
Because F1-1 = 0 and F1-2 = 1 for this geometry, an energy balance requires that q"1 = -q"2, and δ11 = 1 and δ12 = 0, this reduces to σ ( T14 − T24 ) 1 (1 − ε 2 ) σ ( T14 − T24 ) = + q1′′ = q1′′ or ε2 1 1 ε1 + − 1 ε1 ε 2
808 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
By invoking the energy balance that q"1 = -q"2, it wasn't necessary to write eq. (10.46) for surface k = 2. Generally, it's better to apply eq. (10.46) for each surface, and then use the overall energy balance as a check on the final results. The value of q" is the rate of energy per unit area transferred by radiation; the sign convention chosen [Fig. (10.6)] is that q" is the heat flux added to the surface to balance the radiative transfer. Applying eq. (10.46) to the case of infinitely long concentric circular cylinders results in another useful relation σ (T14 − T24 ) q1′′ = (10.47)
D 1 + 1 − 1 ε1 D2 ε 2 1
Example 10.6 A perfectly adiabatic large horizontal surface (1) is exposed to air at a temperature of Tair = 77 oC. Surface 1 has an emissivity of ε1 = 0.2. A second large surface (2) is parallel to the adiabatic surface. The second surface has a black coating, and is at temperature T2 = 450K. The free convection between surface 1 and the air is given by the relation q"conv = h(ΔT)(T1-Tair), where h(ΔT ) = [(T1 − Tair ) / 3] (W / m 2 ⋅ K ) . a) What is the temperature of surface 1? b) What temperature would surface 1 have if surface 2 were at the same temperature as surface 1, with all other conditions the same? c) What temperature would surface 1 have if the air temperature were at the same temperature as surface 1, with all other conditions the same? Solution: The energy balance on surface 1 is (using the result of Example 10.5); σ ( T14 − T24 )
′′ ′′ qrad + qconv =
4 1
1
= ε1σ ( T − T
ε1
+
1
4 2
) + h ( ΔT )(T − T ) = 0
1 air
ε2
−1
+ h ( ΔT )( T1 − Tair )
0.2σ ( T14 − T24 ) + 0.2σ T14 +
(T1 − Tair )
3
2
=0
T2 T12 2 − T1Tair − 0.2σ T24 − air = 0 33 3
Rewrite as
T1 = T 0.3 T 2 0.3 T2 0.3 σ T14 + 1 − σ T24 − air = σ T14 + 1 + (173) Tair 2Tair Tair 2 700 700
Chapter 10 Heat Transfer by Radiation 809
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Trial and error (plug in guess on RHS, get new guess) or the use of a root finding routine gives T1 = 382.1K b) By inspection, T1 = Tair. c) By inspection, T1 = T2 Equation (10.46) can be solved for a larger number of surfaces as follows. Taking the case of the temperatures of all N surfaces known, the equation becomes N δ (1 − ε j ) q′′ = N F σ T 4 − T 4 (10.48) kj − Fk − j ( k j) ε ε j j j =1 k − j j =1 j Let
ε
j =1
N
δ kj
j
− Fk − j
(1 − ε ) q′′ = a
j
εj
j
′′ kj q j ;
F
j =1
N
k− j
σ ( Tk4 − T j4 ) = C k .
(10.49)
Then akj =
δ kj εj
− Fk − j
(1 − ε ) .
j
εj
For an enclosure of N surfaces, the set of
equations becomes
′′ a11q1′′ + a12 q2 + ... + a1 j q′′ + .......a1N q′′ = C1 j N ′′ ′′ ′′ a21q1 + a22 q2 + ... + a2 j q′′ + .......a2 N qN = C 2 j . . . . . . . . . . . . . . . . . . . .
(10.50)
. ′′ aN 1q1′′ + aN 2 q2 + ... + aNj q′′ + .......aNN q′′ = C N j N
In matrix form, this is Aq" = C, where the A matrix is
a11 + a12 + ... + a1 j + .......a1N a21 + a22 + ... + a2 j + .......a2 N .. .. .. aN 1 + aN 2 + ... + aNj . . . . . . .. .. .. + .......aNN
A=
(10.51)
The solution for q"j is found using standard matrix inversion routines to find the inverse of matrix A, denoted by A-1, and the solution for q" is q" = A -1C (10.52) Example 10.7 An infinitely long rectangular enclosure is shown in the Fig. 10.11. Find the heat transfer q at each surface. Solution The configuration factors can be found using the crossed-string method (Example 10.4). This gives
810 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
A1, T1 = 600K, ε1 = 1.0
A2, T2=1000K ε2 = 0.80 A3, T3=750K ε3 = 0.65 W=2m
H=1 m A4, T4=300K ε4 = 0.50
Figure 10.11 Rectangular enclosure
H 2 1/ 2 H F2 − 4 = + 1 − = 0.52 + 1 − 0.5 = 0.618 . W W It follows that F2 −3 = F2 −1 = (1 − F2 − 4 ) / 2 = 0.191 . Similarly,
F1− 3
1/ 2
W 2 = + 1 H
1/ 2
1/ 2 W − = 22 + 1 − 2 = 0.236 and F1-4 = F1-2 = 0.382. H
All remaining factors are easily found using reciprocity. δ (1 − ε j ) , the A matrix becomes Using akj = kj − Fk − j εj εj 1 A = 0 0 0 -0.095 1.25 -0.095 -0.155 -0.127 -0.103 1.538 -0.103
N
-0.382 -0.618 -0.382 2
The C vector has elements [eq. (10.49)] C k = σ Fk − j ( Tk4 − T j4 ) , giving
j =1
-1.872 × 104 5.159 × 104 C= -5.628 × 103 -3.942 × 104 Then, inverting the A matrix to give A-1, q" = A-1C results in q" = -2.312×104 3.211×104 -6.019×103 -1.754×104
Chapter 10 Heat Transfer by Radiation 811
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
To check the results, an overall energy balance on the enclosure gives (remember that energy rate and not heat flux must sum to zero, so we need to multiply each flux by its respective area)
q j = q′′j A j
j =1 j =1
4
4
= −2.312 × 104 × 1 + 3.211 × 104 × 2 − 6.019 × 103 × 1 − 1.754 × 104 × 2 = 0
Radiation shields An important application of radiative transfer between infinite parallel plates and between infinitely long concentric circular cylinders is in the construction of multilayer insulation. Multiple closely-spaced sheets of low-emissivity material are used to make highly insulating light-weight blankets for many uses, including thermal blankets for campers and fire protection blankets as an emergency shelter for forest fire fighters. Consider N shields of very thin material coated on each side with a low emissivity coating of emissivity ε as shown in Fig. 10.12. Let the upper surface 1 be at temperature T1, and the lower surface at T2. The top and bottom plates have the same emissivity as the shields. Treating the radiative transfer between any two layers as a series resistance heat transfer problem, the result of Example 10.5 can be placed in the form σ ( T j4 − T j4+1 ) σ ( T j4 − T j4+1 ) σ ( T j4 − T j4+1 ) ′′ = = = (10.53) q
1 1 ε + ε − 1 2 ε − 1 R
where R = [(2/ε)-1] is the radiative resistance. For N=1, there will be two gaps σ ( T14 − T24 ) , or, in general for N shields between between the plates, so q′′ =1 = N
2R
the two bounding plates,
′′ qN shields =
σ ( T14 − T24 )
( N + 1) R
(10.54) 1 S1
Sj-1 Sj Sj+1 SN 2
Figure 10.12 Multiple radiation shields
812 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
The ratio of heat transfer with N shields in place to that with no shields is
′′ qN shields 1 = ′′ qno shields ( N + 1)
(10.55)
One shield cuts the radiative transfer in half, and 10 shield layers drops the radiative transfer to 9 percent of the unshielded value. Spacecraft use this type of insulation blanket, since in the vacuum of space there is no gas between the shield layers and conduction and convection heat transfer that may be present between the layers on Earth are suppressed. (Actually, thin blankets may use a dimpled surface on each layer to provide spacing between layers, and the small dimples provide a minor conduction path even in space.) Insulation blankets to minimize temperature swings of electronics and other temperature-sensitive packages exposed to sun-shade orbits usually have 5 to 30 layers. The multilayer blankets are also used to shield cryogenic vessels, where environmental temperatures near 300K can otherwise provide a significant radiation flux to cold cryogenic tank surfaces. The NASA James Webb space telescope will view objects mostly in the infrared spectrum (0.6 to 28 μm), and so must be maintained at a very low temperature (below 50K) to minimize IR signal interference from local radiation sources. A large five-layer radiation shield is used to minimize solar-induced temperature gradients and any resulting thermal distortions of the telescope mirrors as well as to maintain the required low temperature. The shield layers are made of Kapton with an aluminum/doped-silicon coating on the front face (for good radiative properties as well as micrometeoroid protection) and a thin aluminum layer on the other surfaces. The orbit of the telescope will allow the shield to block radiation from the Sun, Moon, and the Earth. A photo of the development shield is shown in Fig. 10.13.
Figure 10.13 Multilayer shield development for the James Webb space telescope. The full-scale sunshield is about 22 by 12 m.(Courtesy NASA)
Chapter 10 Heat Transfer by Radiation 813
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
10.2.3 Nongray Surfaces
If the surfaces exchanging energy have wavelength-dependent properties, then the absorptivity and emissivity of the surfaces are not equal as for the gray surfaces considered in Section 10.2.1. If all surfaces in the enclosure are close in temperature, then the total absorptivity and total emissivity of a given surface may not differ significantly; however, in systems such as solar concentrators, radiant furnaces with high-temperature heating elements, and systems with combustion sources, this will not be the case, and the wavelength dependence must be considered. The net radiation equations for a small wavelength interval dλ were given in eqs. (10.36), (10.37) and (10.39): J λ ,k = emitted + reflected flux=ε λ ,k E λb, k + ρλ , k Gλ , k (10.56) ′′ (10.57) dqλ ,k = outgoing-incoming radiation = J λ ,k d λ − Gλ ,k dλ
Gλ ,k = J λ , j Fk − j
j =1 N
(10.58)
As for the simplified case of gray surfaces, these can be combined (retaining the spectral dependence) to give (1 − ε λ ,k ) dq′′ J λ ,k dλ = Eλb dλ − (10.59) λ ,k ε λ ,k and
′′ J λ k dλ = dqλ ,k + J λ , j Fk − j d λ
j =1 N
(10.60)
The configuration factors are not wavelength dependent, but their use does require that all surfaces be treated as diffuse. As for gray surfaces, the radiosity can be eliminated from the equation set, resulting in N N δ (1 − ε λ , j ) kj Fk − j ( Eλb,k −Eλb, j )dλ = ε − Fk − j ε dqλ′′, j (10.61) j =1 j =1 λ , j λ, j If all temperatures Tk are given as boundary conditions, then eq. (10.61) can be solved by dividing the spectrum into M finite wavelength intervals Δλm, and solving for the Δq"j,m in each interval. The total heat flux for each surface is then found from
q′′ = Δq′′,m j j
m =1 M
(10.62)
For the case of all surface temperatures known, the problem is reduced to solving the gray problem in M wavelength intervals, and then summing the results. Example 10.8 Two parallel plates are of width W = 2m and are separated by a distance H = 1m as shown in Fig. 10.14. Plate 1 is at T1 = 800 K and plate 2 is at T2 = 300K. The spectral emissivities of the diffuse plates are
814 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
shown in Fig. 10.15. Find the heat fluxes q"1 and q"2. Assume for simplicity that the radiosity of each plate is uniform. Solution: The net radiation equation (10.61) is written for the wavelength intervals 0-3 μm, 3-9 μm, and 9-∞ μm. These spectral regions are called A, B, and C. For region A, using the blackbody fractions [Eqs. (9.11)-(9.13)], eq. (10.61) becomes, using (h = H/W)
F
j =1
N
k− j
( Eλ b,k −E λb, j )dλ
Spectral region A
= Fk − j
j =1
N
3
λ =0
( E λb,k − Eλb , j )d λ
N N δ (1 − ε A, j ) Δq′′ kj A, j = Fk − jσ (Tk4 F0 − 3Tk −T j4 F0 −3T j ) = − Fk − j ε A, j j =1 j =1 ε A , j
The configuration factor between the plates is
F1− 2 = F2 −1 = h 2 + 1
1/ 2
− h = 0.52 + 1
1/ 2
− 0.5 = 0.6180 .
x2
W
T2, ε2
T∞ = 0 H
T1, ε1
x1
Figure 10.14 Radiation heat transfer between two parallel plates
1.0 0.8 ελ,1 0.5 ελ,2 0.3 3.0 0.8
0
λ, μm
0
9.0
λ, μm
Figure 10.15 Spectral emissivity of the two plates
Chapter 10 Heat Transfer by Radiation 815
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Expanding the net radiation equation for surface k = 1 gives F1− 2σ ( T14 F0 − 3T − T24 F0 −3T ) + 2 F1− 3σ ( T14 F0 − 3T ) = σ (T14 F0 − 3T − F1− 2T24 F0 − 3T )
1 2 1 1 2
1 (1 − ε A,2 ) Δq′′ = Δq′′ ,1 − F1− 2 A A ,2 ε ε A ,2 A ,1
Use was made of the configuration factor algebra relations F1-1 = F2-2 = 0, F1-3 = F1-4, and F1-2 +2F1-3 = 1. Surfaces 3 and 4 are the open ends of the enclosure. For surface k = 2, 1 (1 − ε A,1 ) Δq′′ σ ( T24 F0 −3T − F2 −1T14 F0 −3T ) = Δq′′ ,2 − F2 −1 A A ,1 ε ε A ,1 A ,2
2 1
The two relations for k = 1 and 2 in spectral region A have the two unknowns Δq"A,1 and Δq"A,2, so it is unnecessary to write relations for k = 3,4 unless Δq"A,3 and Δq"A,4 are desired for an energy balance check. Similar relations can be derived for spectral regions B and C to give the final equation set: 1 (1 − ε A,2 ) Δq′′ σ (T14 F0 −3T − F1− 2T24 F0 −3T ) = A: Δq′′ ,1 − F1− 2 A A ,2 ε ε A ,2 A ,1
1 2
σ ( T24 F0 −3T − F2 −1T14 F0 −3T ) =
2 1
1 (1 − ε A,1 ) Δq′′ Δq′′ ,2 − F2 −1 A A ,1 ε A ,1 ε A ,2 1
2 2
B:
σ (T14 F3T −9T − F1− 2T24 F3T −9T ) = Δq′′ − F1− 2 ε B ,1 B ,1
1 1
(1 − ε ) Δq′′
B ,2
ε B ,2
B ,2
σ (T24 F3T −9T − F2 −1T14 F3T −9T ) = Δq′′ − F2 −1 ε B ,2 B ,2
2 1 1
1
(1 − ε ) Δq′′
B ,1
ε B ,1
B ,1
C:
σ (T14 F9T −∞ − F1− 2T24 F9T −∞ ) = ε
1 2 2
1 (1 − ε C ,2 ) Δq′′ ′′ ΔqC ,1 − F1− 2 C ,2 ε C ,2 C ,1
σ (T24 F9T −∞ − F2−1T14 F9T −∞ ) = Δq′′ − F2−1 ε C ,2 C ,2
2 1 1
1
(1 − ε ) Δq′′
C ,1
ε C ,1
C ,1
Solving for the heat flux in each spectral range results in Spectral Δq"1 Δq"2 (W/m2 ) (W/m2 ) range 0-3 μm 2444 -1308 3-9 μm 7547 -3989 9-∞ μm 1737 - 347
ΣΔq"
11,728
-5644
The total heat flux from summing the flux in each spectral range is thus ′′ q1′′ = 11,728 W/m2, and for surface 2, q2 = –5644 W/m2.
816 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Suppose now that, rather than specifying the temperatures of all surfaces, the total heat flux q"k is given as a boundary condition on one or more surfaces. Equations (10.56) through (10.62) still apply: Now, however, the left hand side of eq. (10.62) is a known value for the one or more surfaces with specified radiative heat flux. To solve for the unknown temperatures on these surfaces requires an iterative solution. The temperatures for the flux-specified surfaces must be guessed, and the M values of all heat fluxes are then found. The computed spectral flux values for the flux-specified surfaces are substituted into the right-hand side summation in eq. (10.62), and the result checked against the specified radiative flux for that surface. If it differs, then new temperatures for the flux-specified surfaces are assumed, and the solution is repeated until a converged result is reached.
10.2.4 Surfaces with Varying Temperature, Radiative Flux, or Properties
One of the stipulations for use of the configuration factor is that the radiosity J must not vary across the surface described by the configuration factor. In many cases of finite surfaces, this requirement is not met; in particular, near the intersection of two surfaces, the irradiation G is likely to vary considerably, so that even if ε and T are constant across the surface, the radiosity J = εσ T 4 + ρ G will vary. In such a case, a surface must be divided into ever smaller increments so that each increment will meet the requirement of uniform radiosity as shown in Fig. 10.16. In the limit, the net radiation equations must be rewritten for a differential element dAk on surface k at location (rk ) . Equations (10.40) and (10.41) for that element k are J k ( rk ) = ε k ( rk ) σ Tk4 ( rk ) + (1 − ε k ( rk ) ) Gk ( rk ) (10.63) (10.64) Because the radiosity may vary over each surface, the irradiation arriving at surface k from a given surface j is now in the form of an integral, G j →k (rk ) = J j (r j )dFdk − dj (rk , r j ) (10.65) A
j
′′ qk (rk ) = J k (rk ) − Gk (rk )
and Equation (10.42) must now sum the irradiation coming from each finite surface j in the enclosure to become
Gk (rk ) = J j (r j )dFdk − dj (rk , r j )
j =1 Aj N
(10.66)
The result is that we now must solve 3N simultaneous integral equations rather than 3N simultaneous algebraic equations as for the uniform-radiosity case. The problem can again be reduced to eliminate the irradiation Gk (rk ) to give
1 − ε k ( rk ) q′′ r J k ( rk ) = σ Tk4 ( rk ) − k( k) ε k ( rk )
(10.67)
Chapter 10 Heat Transfer by Radiation 817
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
A2
Aj
rj
dAj, εj(rj) Tj(rj) or q"j(rj)
A1
dFdk-dj
AN
Ak rk
Figure 10.16 Enclosure geometry for surfaces with varying properties and radiation variables.
and
′′ J k ( rk ) = qk ( rk ) + J j ( r j )dFk − j ( rk , r j )
j =1 Aj N
(10.68)
Finally, eliminating the radiosity from eqs. (10.67) and (10.68) gives the set
σ T ( r ) − T ( r )dF ( r , r )
4 k
N
k
4 j
j
j =1
Aj
dk −dj
k
j
=
ε k ( rk )
′′ qk ( rk )
−
j =1
N
1 − ε j ( rj )
Aj
(10.69)
q′′ ( r j ) dFdk −dj ( rk , r j ) j
ε j ( rj )
Example 10.9 Two infinitely long parallel plates of finite width W are separated by a distance H (Fig. 10.17). Plate 1 has uniform emissivity and temperature ε1 and T1, while plate 2 has uniform emissivity and temperature ε2 and T2. The environment around the plates is at T∞ = 0 K. Derive expressions for the heat flux distribution on each plate. Solution: Write eq. (10.69) for surface k = 1:
σ T
j =1 Aj
4
4 1
− T j4 dFd1−dj ( x1 , r j ) −
j =1
4
=
q1′′( x1 )
1− ε j
Aj
ε1
εj
q′′ ( r j ) dFd1−dj ( x1 , r j ) j
818 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
x2 W
T2, ε2 T∞ = 0 H
T1, ε1
x1 Figure 10.17 Radiation heat transfer in a system with two parallel plates
For the open ends, ε = 1 and T = 0 K; also dFdk-dk = 0 for a plate. In addition, the temperatures and properties are constant across all surfaces, so they do not depend on x1 or x2 and can be taken outside the integrals. By using the summation rule and the definition of element-area onfiguration factors, the result for the left hand side is
σ T
j =1 Aj
4
4
4 1
− T j4 dFd1−dj ( x1 , r j )
4
= σ T14 dFd1−dj ( x1 , r j ) − σ T j4
j =1 A j j =1
Aj
dF ( x , r ) = σ T
d 1−dj
1
j
4 1
− σ T24 Fd1− 2
With those simplifications, the equation reduces to q′′( x ) 1 − ε 2 σ T14 − σ T24 Fd1− 2 = 1 1 − q′′ ( x ) dF ( x1 , x2 ) ε1 ε 2 A 2 2 or q1′′( x1 ) 1− ε2 q′′ ( x ) dF = σ T14 − σ T24 Fd1− 2 + ( x1 , x2 ) ε1 ε 2 A 2 2 A similar derivation for surface k = 2 gives ′′ q2 ( x 2 ) 1 − ε1 q1′′( x1 ) dF = σ T24 − σ T14 Fd 2 −1 + ( x2 , x1 ) ε2 ε1 A This gives two integral equations in the two unknown distributions q"1(x1) and q"2(x2). The configuration factors for this geometry (two infinitely long parallel strip elements) can be found using the crossed strings method, or from Configuration C-3 of the configuration factor catalog (Howell, 2003), using h = H/W,
2 d 1−d 2 2 d 1− d 2
1 d 2− d 1
Chapter 10 Heat Transfer by Radiation 819
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
dFd1− d 2 ( x1 , x 2 ) =
h2 2 h 2 + ( x 2 − x1 ) 2
3/ 2
dx 2
The values of q"1(x1) and q"2(x2) can be found by assuming a form for q"2(x2) and substituting it into the first of the two integral equations to generate an approximate distribution q"1(x1). This is then substituted into the second integral equation to generate a new q"2(x2). This process is continued until convergence is achieved. Usually, the surfaces are divided into a finite number of elements and the integrals can then be replaced by summations.. If each surface is divided into M elements of equal size Δx, then the problem becomes quite similar in form to the constant radiosity solutions. For this problem, the two integral equations become q1′′,m 1− ε2 M ′′ = σ T14 − σ T24 Fm − 2 + q2,n dF ε1 ε 2 n =1
m −n
′′ q2, n
ε2 ε1 m =1 or, substituting the configuration factor relation,
′′ q1,m
= σ T24 − σ T14 Fn −1 +
M
1 − ε1
q′′
M
1, m
dFn−m
ε1
= σ T14 − σ T24
n =1
h 2 Δx 2 h 2 + ( m − n ) 2 ( Δx ) 2 h 2 Δx
3/ 2 3/ 2
+
1− ε2
ε2
′′ q2,n
n =1
M
2 h 2 + (m − n ) 2 (Δx ) 2
M
′′ q2,n
ε2
= σ T24 − σ T14 + 1 − ε1
h 2 Δx 2 h 2 + ( m − n ) 2 ( Δx ) 2 h 2 Δx
3/ 2 3/ 2
m =1
ε1
m =1
q′′
M
1, m
2 h 2 + ( m − n ) 2 ( Δx ) 2
These two equations can be written in matrix form (giving 2M equations) and solved directly, or the trial and error method proposed for the integral form can be used. Example 10.10 For the geometry of Example 10.9, let ε1 = ε2 = 1, T1 = 1000K, T2 = 300K, and H = 1 m, W = 1 m, find the heat flux distribution on each surface. Solution: The final two equations of Example 10.9 become
′′ q1,m = σ T14 − σ T24
n =1 M M
h 2 Δx 2 h 2 + (m − n ) 2 (Δx ) 2 h 2 Δx 2 h 2 + (m − n ) 2 (Δx ) 2
3/ 2 3/ 2
′′ q2, n = σ T24 − σ T14
m =1
820 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Figure 10.18 Heat flux distribution
Examining the two equations show that because of the choice of black surface, the heat fluxes on the two surfaces are not coupled, and can be computed independently. For h = H/L = 1 and choosing M = N = 30, the heat flux distributions are found, and shown in the Fig. 10.18 (using x = [(Δx/2) + mΔx]) Because of the sign convention, the q1"(x1) is positive because energy must be added to surface 1 to maintain its higher temperature, and q2"(x2) is negative because energy must be removed from surface 2 to maintain its lower temperature. Examples 10.9 and 10.10 demonstrate why the assumption of uniform radiosity is invoked in many radiation problems. Extension to the required integral equations and numerical solution results in a major increase in computational effort. However, if the loss of accuracy by assuming constant radiosity is not acceptable, the more exact formulation is required.
10.3 Multimode Heat Transfer with Radiation
When radiation is combined with conduction and/or convection, the problems are called conjugate or multi-mode heat transfer. If temperature is one of the unknowns in the problem being studied, then such problems are highly nonlinear, in that the fundamental radiation terms involve temperature to the fourth power, while conduction and convection depend on the first power of temperature or its derivatives. Property variations with temperature may add further nonlinearities.
Chapter 10 Heat Transfer by Radiation 821
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
W
δ
x2 H x1 Surface 1, k, ε Surface 2, k, ε, q2(x2) = 0 Te = 0 K
qe(x1) = 200 W/m2
Figure 10.19 Infinitely long parallel conducting plates with thickness δ and thermal conductivity k
The mathematics of highly nonlinear equations is not well-developed, and almost all combined mode problems are solved numerically; closed-form analytical solutions seldom exist. As an example of a combined-mode radiationconduction problem, consider the geometry and conditions shown in Fig. 10.19. Radiation transfer occurs between the parallel plates and from the plates to the environment. Because the plates are not isothermal, conduction occurs within the plates, and will modify the temperature distributions, which are the unknowns in the problem. For this case, the usual net radiation equations are as derived in Example 10.10: ′′ q1,rad ( x1 ) 1− ε = σ T14 − σ T24 Fd1− 2 + q′′ ( x ) dF (10.70) ( x1 , x2 ) ε ε A 2,rad 2 ′′ q2,rad ( x 2 ) 1− ε = σ T24 ( x 2 ) − σ T14 ( x1 ) Fd 2 −1 + q′′ ( x ) dF ( x2 , x1 ) ε ε A 1,rad 1 (10.71) Note that the q"k,rad is the heat flux required to balance the radiation at the surface. It can be supplied by conduction, convection, internal generation, or a combination of these. In this case, only conduction is present. If conduction is assumed to occur only in the x-direction in the thin plates, then two additional energy equations are written for the conducting plates
2 d 1−d 2 1 d 2− d 1
′′ ′′ q1,rad = qe + kδ
d 2T dx
2 1
1
(10.72) (10.73)
′′ q2,rad = kδ
d 2T dx
2 2
2
Equations (10.72) and (10.73) contain second-order derivatives, so two boundary conditions are required for each. The appropriate ones are
822 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
dT dT = = 0 at x , x = W / 2 dx dx
1 2 1 2 1 2
(10.74)
2
k
dT dT = εσ T ( x = 0); k = εσ T ( x = 0) dx dx
1 4 2 4 1 1 2 1 2
(10.75)
Equation (10.74) is a symmetry condition, and eq. (10.75) equates the conduction into the plate ends with the radiative loss from the plate end. Equations (10.72) and (10.73) can be substituted into eqs. (10.70) and (10.71), resulting in two simultaneous integral-differential equations with the unknowns T1(x1) and T2(x2). The resulting Equations (10.76) are highly nonlinear, involving temperatures to the fourth power as well as second derivatives in temperature. 2 W dT 1 1 − ε d T2 kδ q + kδ = σ T14 ( x ) − σ T24 ( x ) − dFd1− d 2 0 ε dx dx ε (10.76) W d T 1 d 2T 1− ε 4 4 kδ = σ T2 ( x ) − σ T1 ( x ) − q + kδ dFd1− d 2 0 ε dx dx ε To simplify eqs. (10.76), define kδ 1/ 4 Tref = ( qe / σ ) ; N CR = ,θ = T / Tref ; X1 = x1 / W ; X 2 = x 2 / W 3 σ Tref W 2 (10.77) and eqs. (10.76) become 1 dθ d 2θ 2 4 4 1 + N dX = εθ1 ( X ) − 0 εθ 2 ( X ) − (1 − ε ) N CR dX dFd1− d 2 (10.78) 2 1 d θ dθ 4 4 N = εθ 2 ( X ) − 0 εθ1 ( X ) − (1 − ε ) 1 + N dFd1− d 2
2 1 e 2 1 2 2 1 2 2 2 1 2 2 1 e 2 2 1
2
1
CR
2
1
2
2
1
2
2
2
1
CR
dX
2 2
2
1
CR
dX 1
2
with boundary conditions [from eqs. (10.74) and (10.75)]
dθ1 dX1 N CR N CR
X1 =1/ 2
= dθ1 dX1 dθ 2 dX 2
dθ 2 dX 2
X 2 =1/ 2
=0
4 = εθ1 ( X1 = 0 ) X1 = 0 4
δ W δ W
(10.79)
= εθ 2 ( X 2 = 0 )
X2 =0
In terms of the nondimensional variables, the required configuration factor is (using h = H/W),
dFdX − dX =
1 2
h 2 dX 2
2 ( X 2 − X 1 ) + h 2
2 3/ 2
(10.80)
The parameters that appear in the problem are the emissivity ε, the conduction/radiation parameter NCR, and the two geometric parameters δ/W and h. The NCR or modifications of it often appear in radiation-conduction problems.
Chapter 10 Heat Transfer by Radiation 823
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Its meaning can be appreciated by multiplying the numerator and denominator by Tref to give kδ ( Tref / W ) kδ Tref N CR = = (10.81) 4 2 4 σ Tref W σ Tref W The last form in eq. (10.81) is the ratio of a conduction heat transfer to a radiative emission. The conduction is through the plate thickness δ with a temperature gradient of (Tref – 0)/W. The radiation is emission from a black surface of width W and temperature Tref. The NCR is thus a measure of the importance of conduction to radiation in affecting the temperature distributions on the surfaces. For example, if NCR = 0, eqs. (10.78) reduce to the radiation-only solutions. If NCR is very large, then the problem becomes a pure conduction problem. (This is shown be reverting to eqs. (10.72) and (10.73) and placing them in dimensionless form). It is unlikely that eqs. (10.78), subject to boundary conditions eqs. (10.79), will have a closed-form analytical solution, and a numerical approach must be employed.
10.3.1 Numerical Methods
To solve the radiation-conduction problem of eqs. (10.78)-(10.80), the integrals in the equations are replaced by summations over elements of width Δx, i.e.,
N ⋅Δxn
0
f ( xm , xn )dxn ≈ f m , n Δxn
n =1
N
(10.82)
Substituting into eq. (10.78) and using the finite difference form for the second derivatives results in N θ1,m −1 − 2θ1,m + θ1, m +1 4
1 + N CR = εθ1,m − Φ1,m ,n ΔX 2 2 n =1 ( ΔX1 ) (10.83) M θ 2,n −1 − 2θ 2,n + θ 2,n +1 N CR = εθ 2,n 4 − Φ 2,m , n ΔX1 2 m =1 ( ΔX 2 )
θ 2, n −1 − 2θ 2, n + θ 2, n +1
where, including the configuration factor relations, the Φs are
4 Φ1, m , n = εθ 2, n − (1 − ε ) N CR
( ΔX )
2
2
3/ 2 2 ( n ⋅ ΔX − m ⋅ ΔX )2 + h 2 2 1 3/ 2 2 2 ( m ⋅ ΔX1 − n ⋅ ΔX 2 ) + h 2
h2
h2
Φ 2, m , n = εθ1,4m − (1 − ε ) 1 + N CR
θ1, m −1 − 2θ1, m + θ1, m +1
( ΔX )
1
2
(10.84) The boundary conditions become
824 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
θ
1,(
N 2
−1)
=θ
1,(
N 2
+1)
;θ
2,(
M 2
−1)
=θ
2,(
M 2
+1)
; M and N even
N CR N CR
θ1, 2 − θ1,1
( ΔX )
1 2
2
4 = εθ1,1
N δ / W + 1 − 1 − Φ1,1, n ΔX 2 n =1 ΔX1
(10.85)
θ 2, 2 − θ 2,1
( ΔX )
2
4 = εθ 2,1
δ / W M + 1 − Φ 2, m ,1 ΔX1 ΔX 2 m =1
Because the boundary elements at n = 1 and m = 1 are finite with width = ΔX, the radiative exchange to the ΔX surfaces must be included in addition to the radiation-conduction balance at the end surfaces. The ΔX approaches zero in the analytical boundary conditions of eq. (10.79) so these terms do not appear. The result is (M+N) nonlinear algebraic equations with N values of θ1,n and M values of θ2,m. This set can, in principle, be solved by solvers available for treating a set of non-linear algebraic equations. However, experience indicates that such solvers often fail unless some further steps are taken to rearrange the form of the equation. Let ΔX1 = ΔX2 (equal sized increments on each plate, so that eqs. (10.83) (10.85) become
N θ1, m −1 − 2θ1, m + θ1, m +1 4 1 + N CR = εθ1,m − Φ1, m , n ΔX 2 n =1 ( ΔX ) M θ 2, n −1 − 2θ 2, n + θ 2, n +1 4 N CR = εθ 2,n − Φ 2, m , n ΔX 2 m =1 ( ΔX )
4 Φ1, m , n = εθ 2, n − (1 − ε ) N CR
(10.86)
θ 2, n −1 − 2θ 2, n + θ 2, n +1
( ΔX )
3/ 2
2
×
h2
2 ( n − m ) ΔX 2 + h 2
2
Φ 2, m , n
×
4 θ1, m −1 − 2θ1, m + θ1, m +1 = εθ1, m − (1 − ε ) 1 + N CR ( ΔX )2
h2
2 ( m − n ) ΔX 2 + h 2
2 3/ 2
(10.87)
θ
1,(
N 2
−1)
=θ
1,(
N 2
+1)
;θ
2,(
M 2
−1)
=θ
2,(
M 2
+1)
; M and N even
N CR N CR
2
θ1,2 − θ1,1
( ΔX )
2
4 = εθ1,1
N δ / W + 1 − 1 − Φ1,1, n ΔX ΔX n =1
(10.88)
θ 2,2 − θ 2,1
( ΔX )
2
4 = εθ 2,1
δ / W M + 1 − Φ 2, m ,1 ΔX ΔX m =1
The parameter NCR/(ΔX) now appears in eqs. (10.86) - (10.88), and is a guideline for how to approach a numerical solution.
Chapter 10 Heat Transfer by Radiation 825
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
10.3.2 Conduction Dominated Problems
If conduction dominates (i.e., NCR/(ΔX)2>>1), then eq. (10.86) can be rearranged to reflect a conduction problem that is perturbed by the presence of radiation:
θ1,m −1 − 2θ1,m + θ1,m +1 = θ 2,n −1 − 2θ 2, n + θ 2,n +1 =
( ΔX )
N CR
2
N 4 εθ1,m − 1 − Φ1,m ,n ΔX n =1 M 4 εθ 2,n − 1 − Φ 2, m ,n ΔX m =1
( ΔX )
N CR
2
(10.89)
and for this case the right-hand-side of eq. (10.89) will be small and the solution will be close to the conduction-only solution. Equation (10.89) (including the insertion of the boundary conditions) is a set of equations of the matrix form Aθ(x) = C(x) (10.90) where A is a tridiagonal matrix of coefficients that need only be inverted once. The solution is then θ(x) = A -1C(x) (10.91) Equation (10.91) is solved by assuming the distributions of θ1,n and θ2,m, using these to evaluate Φ1,m,n and Φ2,m,n, which in turn are used to evaluate C(x). The matrix multiplication indicated in eq. (10.91) gives a new set of θ values; the process is repeated until convergence. Effectively, this method is solving a set of linear equations by interation in the unknown values of θ, and amounts to a series of matrix multiplications.
10.3.3 Radiation Dominated Problems
If radiation dominates so that NCR/(ΔX)2 1, then eq. (10.86) can be rearranged to
θ1,m 4 = θ 2,n 4
− 2θ1, m + θ1,m +1 θ 1N Φ1, m ,n ΔX − 1 + N CR 1,m −1 2 ε n =1 ( ΔX )
− 2θ 2, n + θ 2, n +1 θ 1M = Φ 2,m ,n ΔX − N CR 2,n −1 2 ε m =1 ( ΔX )
(10.92)
Again, the set of equations can be arranged as a matrix equation of the form A1 ⋅ θ4 = C1 (θ) (10.93) Matrix inversion of A1 gives -1 θ4 = ( A1 ) C1 (θ) (10.94) Again, the A1 matrix need be inverted only once. Equation (10.94) is solved by assuming the distributions of θ1,n and θ2,m, using these to evaluate Φ1,m,n and Φ2,m,n and the second derivative terms, which in turn are used to evaluate C1(x). The matrix multiplication indicated in eq. (10.94) gives a new set of θ values; the
826 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
process is repeated until convergence. Effectively, this method is iteravely solving a set of linear equations in the unknown values of θ4.
10.3.4 Problems with Both Modes Significant
When NCR is not near either the large or small limit, then the problem is truly nonlinear, and the solution methods described for small or large NCR/(ΔX)2 values often fail. In that case, eq. (10.86) can be arranged as N θ1,m −1 − 2θ1,m + θ1, m +1 4
εθ1,m − N CR = 1 + Φ1, m , n ΔX 2 n =1 ( ΔX ) − 2θ 2, n + θ 2, n +1 M θ εθ 2,n 4 − N CR 2,n −1 = Φ 2,m ,n ΔX 2 m =1 ( ΔX )
A 2 θ 4 + Bθ = C 2 ( θ )
3
(10.95)
This equation in matrix form becomes
(10.96)
We could define A 3 ( θ ) = A 2 ⋅ θ + B so that eq. (10.96) becomes (10.97) If this approach is followed, an initial guess for the distributions of θ1,n and θ2,m is used to evaluate A3(θ) and C2(θ), A3(θ) is inverted, and a new set of θ values is found. This method has two problems: First, the A3(θ) matrix must be inverted at each iteration, and A3(θ) is probably a full matrix that can require significant time to invert. Second, this method is often quite unstable, and requires a small relaxation factor to be imposed between iterations to avoid divergence. To avoid this, a modified Newton-Raphson iteration method can be used. In this case, eq. (10.96) is rewritten as
A 3 (θ)× θ = C2 (θ)
gj =
M+N k =1
[A
4( 0 ) jk k
θ
+ B jkθ k(0) ] − C j ≡ residual
(10.98)
The residual is a measure of convergence of the solution, and will approach gj = 0 when solution is complete. Next, the function gjk is found
3 ∂ ( LHS of Eq.(10.96) g jk = 4 A jk θ k(0) + B jk ≈ ∂θ k
(10.99)
and gjk is seen to be the gradient in the residual. Now solve the auxiliary equation −1 g jk [ λk ] + g j = 0; [ λk ] = − g jk g j (10.100) The next iterative value for θk is found from θ k( p ) = θ k( p −1) + λk (10.101) (p) The values of θk are now used to find a new residual, and the process is repeated until convergence. The λk values are seen to be those that would make the residual equal to zero on the next iteration if the problem were linear. This
Chapter 10 Heat Transfer by Radiation 827
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
approach has been found to be quite stable, although sometimes a relaxation factor is needed on θk between iterations. In this section, some of the complexities of multi-mode heat transfer when radiation is present have been outlined. A detailed discussion of more advanced numerical techniques is in Hogan and Gartling (2007). There are commercial software packages that can be used for these problems, and they can be particularly useful when convective transfer must be linked to radiation.
10.4 Inverse Problems
It is usually stated that to solve any mathematical problem described by a differential equation, we must have the number of boundary conditions equal to the order of the PDE that is to be solved; in addition, if the problem is transient in nature, an initial condition must also be prescribed. Further, the number of unknowns being sought must be equal to the number of available equations. This is discussed in general in Section 3.3.6. For integral equations, the boundary conditions are usually implicit in the limits of the integrals. Suppose, however, that these conditions are not met; for example, more than the required number of boundary conditions are prescribed on some surface or surfaces, and it is desired to find the missing boundary conditions on other surfaces. In some cases of this type, the number of equations may not match the number of unknowns. These problems fall under the class of inverse problems, because the requirements for standard solution are not met. Problems of this type may be very sensitive to the imposed values of the given boundary conditions, and solutions may vary wildly. This type of problem is referred to as ill-posed. Inverse solutions in heat transfer may be required when experimental observations of temperature or heat flux are not available at the physical location where they are needed, or radiative property distributions in a participating medium must be obtained from remote measurements. For example on the reentry heat shield of a spacecraft (Fig. 10.20), it may be impossible to measure temperature Ts(t) or heat flux q"s(t) on the ablating surface of the shield, but transient temperatures T(t) and heat fluxes q"(t) can be measured on the interior surface. Can these remote measurements be used to determine the unknown temperature and heat flux on the heat shield surface? Intuitively, it would seem so: It isn't so easy! The governing equations for inverse problems tend to be mathematically illposed, and predicting conditions on the remote boundary can result in multiple solutions, physically unrealistic solutions, or solutions that oscillate greatly in space and time. In this Section, we examine why these difficulties arise, and show one solution method for handling them. Various methods may be applied for overcoming the ill-posed nature of the governing equations. Tikhonov (1963) and Phillips (1962) are often credited with developing the first systematic treatment for these types of inverse problems. For problems dominated by conduction, there are texts and monographs available
828 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
T(t) q"(t) Ts(t) q"s(t)
Figure 10.20 Determining the heat shield hot surface temperature (Ts(t) and heat flux q"s(t) ) from measurement of remote temperature and heat flux (T(t) and q"(t))
that demonstrate many of these methods (Tikhonov, 1963; Alifanov, 1994; Alifanov et al. 1995; Beck et al., 1995; Özişik and Orlande, 2000). For determining radiative properties from remote measurements, various inverse techniques have been employed. This type of problem is required in analyzing temperature and soot distributions in flames, accounting for atmospheric attenuation effects in remote sensing from satellites, and is also closely related to problems in x-ray tomography. Aside from determining conditions at an inaccessible boundary or radiative properties, inverse problems arises in the design and control of thermal systems. In these problems, the designer specifies the desired output of the thermal system being designed; in most cases, this is a desired temperature and heat flux distribution over a product located on a design surface. The designer must then predict the necessary energy inputs to the thermal system that will produce the desired distributions over the design surface. For example, a semiconductor wafer of known properties is to be heated following a specified time-temperature curve, while holding the entire wafer at uniform temperature at every instant to avoid thermal stresses or deformation. Thus, the local temperature and heat flux necessary to heat at the prescribed rate are specified at every time and position on the wafer: Two boundary conditions on a single boundary! The unknown inputs may be the required temperature and energy input distributions to radiant heaters, or heater locations, or oven geometry. For thermal systems dominated by radiative transfer, the problem is complicated because the thermal input at any location on the design surface may be affected by some or all radiant sources in the system, depending on the presence of blocking or shading. The mathematical form of the inverse solution
Chapter 10 Heat Transfer by Radiation 829
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
is the set of algebraic or integral equations derived in Section 10.2.4 that must be solved simultaneously. For the inverse problem, the missing boundary conditions cause some of the integrals to be Fredholm integral equations of the first kind, which are notoriously ill-conditioned [Wing (1991), Hansen (1998)]. In addition, depending on how the design surface and the heater surfaces are subdivided, the number of unknowns in the governing radiative exchange equations may be less than, equal to, or greater than the number of equations that describe the system. These factors imply that design and control of distributed radiative sources may be difficult, especially in problems where both a transient temperature and heat flux distribution are prescribed over the design surface, and where other heat transfer modes play a role. The traditional design strategy has been to use trial-and-error solutions with guidance by experience to attain a viable solution. However, this can lead to sub-optimal solutions and can be quite time consuming. Some problems are simply not amenable to trial-and-error approaches. Reviews of inverse methods for design of radiative transfer systems are in França et al. (2002); Daun et al. (2003; 2006). These discuss various methods for solving inverse radiation problems, and cover the work of many contributors.
10.4.1 An Inverse Design Problem
Analytical techniques used for inverse design and control of radiative systems are similar to those used in inverse data analysis, property determination, and remote measurement, but there are significant differences in how the techniques are implemented. The latter types of problem should have a single solution; i.e., we don’t expect a material to have multiple properties that give the same measurements, nor multiple temperatures or heat fluxes at a remote boundary that give the same measurements at an accessible location. However, design problems may allow significantly wider tolerances in specification of acceptable results. In design, solving the inverse problem may produce multiple solutions that fall within the allowable tolerances but are very different in form. Figure 10.21 shows the usual two-dimensional enclosure. Consider an annealing furnace, where a billet of metal is to be heated to a specified temperature. The heat capacity, volume and density of the billet then impose the required heat flux. As the allowable tolerances on q1′′ (x1) and T1(x1) are relaxed, multiple allowable solutions for the heater power distribution may be found. Multiple acceptable solutions are desirable, as they allow the designer to choose among the solutions based on considerations such as smoothness and ease of implementation. It is possible in design problems to specify conditions on the design surface for which no acceptable physical solution on Surface 2 exists. The designer may specify design surface characteristics of q′′ (x1) and T(x1) that cannot be
830 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
x2
W
T2(x2) = ?
T∞ = 0 H
T1(x1) , q1"(x1)
x1
Figure10.21 Simple radiant furnace with design surface boundary conditions specified.
obtained (at least within acceptable error limits around the desired distributions) by any distribution of heater settings. Just because the designer wants a particular outcome, there is no a priori guarantee that it can be obtained! This is in contrast to data analysis problems, where in most cases a feasible solution to the inverse problem is known to exist because some set of physical variables gave rise to the measured experimental data. Finally, designing radiant systems involves inverse solution of integral equations rather than the differential equations that arise in applications involving other heat transfer modes. Since there are few radiation heat transfer problems where conduction and convection can be completely neglected, the governing energy relations are often highly nonlinear integral-differential equations as seen in Section 10.3. Direct Inverse Solutions: A direct or explicit solution to inverse problems requires use of an inverse formulation. Inverse problems are inherently illposed, and the corresponding discretized set of equations is ill-conditioned (i.e., the matrix of coefficients in the discretized solution is singular or near-singular). Ordinary techniques for solving the discretized set of integral equations (e.g., Gauss-Seidel, Gauss elimination or LU decomposition) are likely to either identify non-physical solutions with high amplitude fluctuations and/or imaginary absolute temperatures, or may completely fail to find a solution. Before considering direct inverse solution techniques, the ill-conditioned behavior of a system of equations of the type arising in radiative systems can be diagnosed by carrying out a singular value decomposition (SVD) (Hansen, 1998). The SVD of an arbitrary M×N matrix A is A = USVT, where U and V are orthogonal matrices and S is the diagonal matrix of singular values, and the elements Si,i of the S matrix are arranged so that S1,1 > S2,2 >….>SN,N ≥ 0. The inverse of A is then given by A−1 = VS′UT, where S′i,i= 1/ Si,i. If the condition number of this matrix (S1,1/SN,N ) is large, small singular values dominate the inverse matrix and the solution becomes unstable. If A is rank-deficient, some
Chapter 10 Heat Transfer by Radiation 831
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
of the singular values equal zero and the conventional inversion process fails completely. Most contemporary mathematical packages (Matlab, Mathcad, etc.) contain commands for extracting U, V, and S and the singular values. Example 10.11 For the geometry in Fig. 10.21 for two black surfaces and h = H/W = 0.5, let the dimensionless heat flux be specified as Q"1(x1) = q"1(x1)/σTref4 = 16x12-16x1-6 and T1(x1) = Tref. Using direct matrix inversion, find the required dimensionless temperature distribution on surface 2, T2(x2). Solution: The governing relation for the emissive power distribution on surface k for black surfaces is (eq. (10.69)
σ T ( r ) − T ( r )dF ( r , r ) = q′′ ( r )
4 k k 4 j j j =1 Aj
dk − dj
N
k
j
k
k
In dimensionless form using θ = T/Tref, X = x/W this is
θ ( X ) − θ ( X )dF ( X
4 k k 4 j j j =1 A j
dk − dj
N
k
′′ , X j ) = Qk ( X k )
Writing the equation for surface k = 2 and using the summation relation for configuration factors results in
1 X2
′′ θ 4 ( X 2 ) dFdX − dX dX 2 = θ 2,i dFi − j ΔX 2 = 1 − Q1, j ΔX1 =0 2
1 2
30
i =1
This can be written in discretized form as
A θ
ij i =1
m
4 2, j
= Cj h 2 ΔX 2
2 3/ 2
where A,ij = dFi-j ΔX2 =
2 h + X1,i − X 2, j
2
2 ′′ C j = (1 − Q1, j ) ΔX1 = ( 7 − 16 x1, j + 16 x1, j ) ΔX1 and
{
}
,
1 1 X 2,i = i − ΔX 2 ; X1, j = j − ΔX1 . 2 2
The resulting linear matrix equation Aθ42 = C can now be solved. Note that in the discretized integral equation, the right hand side 1 − Q1′′, j ΔX1 is composed of the given dimensionless temperature on surface 1, and the dimensionless heat flux on that surface is a known function. This is a Fredholm integral equation of the first kind (Hansen, 1998; Wing, 1991)! Using a direct matrix inversion routine (in this case, using MATHCAD which employs a SIMPLEX solver) to obtain θ42 = A-1b gives the result shown in Fig. 10.22. If the values for θ42,i in the figure are introduced into the rearranged net radiation equation:
832 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
5 10
9
Dimensionless fourth power temperature, θ
4
2,j
0
-5 10
9
0
5
10
15
20
25
30
Index, i
Figure 10.22 Predicted dimensionless fourth-power temperature of heater surface 1 using direct matrix inversion.
′′ Q1, j = 1 −
θ
i =1
30
4 2, j
dFi − j ΔX 2
ΔX1
to calculate values of Q1′′( x1 ) , exact agreement is found with the prescribed Q1′′( x1 ) . Thus, the solution shown in the figure is mathematically correct; it solves the mathematical problem exactly. However, note that some large negative values of dimensionless emissive power are predicted, implying imaginary absolute temperatures on surface 2. Observe also that the predicted profile is quite asymmetric, even though the physical solution is expected to be symmetric in form. Therefore, although the solution in the figure mathematically satisfies the problem formulation, it is not a useful solution. The reason for the unphysical behavior found in the solution in Example 10.11 can be determined by viewing the 30 singular values of the A matrix for this problem. MATHCAD returns the values shown in Fig. 10.23. The singular values range from s1,1 = 0.63 at k = 1 to s30,30 = 1.055x10-17 at k = 30, giving a matrix condition number of s1,1/ s30,30 = 6.0x1017. A large condition number implies difficulty in matrix inversion and/or large errors in the solution.
Chapter 10 Heat Transfer by Radiation 833
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
0.1 0.001 10
-5
Singular values, s
k
10 10 10 10 10 10
-7
-9
-11
-13
-15
-17
0
5
10
15
20
25
30
Index, k
Figure 10.23 Singular values of the A matrix in Example 10.11.
Because elements in the inverse of A contain the reciprocals of the singular values, it is easy to see why direct inversion fails. Alternative solution techniques are necessary. To achieve an accurate and physically useful solution, the explicit system of equations must instead be regularized by modifying the ill-conditioned set of equations, or some other approach in formulation must be used. The modified solution will be subject to some error because some information has been deleted or distorted in the regularization process, and the level of error must be selected so that the accuracy of the solution satisfies the designer’s needs. Unlike most finite difference solutions, if the number of increments on the surfaces is increased, the condition number of the A matrix increases and the dimensionless emissive power becomes even more erratic than shown in Example 10.11. For a general transient system, the discretized energy equation for an element on the design surface (where the desired distributions of T and q" are specified) can be written in terms of the emissive power E = εσT4 of surface elements elsewhere in the enclosure. This is essentially an expansion of the multi-mode net radiation approach of Section 10.3. Assuming the elements are numbered so that the first NDS elements are on the design surface, at any time t,
834 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
ρi c pi δ i Ai
dTi ( t ) dt
′′ = qcond ,i ( t ) + qconv, i ( t ) Ai − ( J i − Gi ) Ai ′′
N DS
′′ = qcond ,i ( t ) + qconv ,i ( t ) Ai + J j ( t ) A j F j →i + ′′
j =1
j = N DS +1
N
J j ( t ) A j F j →i − J i ( t ) Ai
(10.102) Applying reciprocity and rearranging, this becomes
j = N DS +1
N
J j ( t ) Fi → j dTi ( t ) dt ′′ + J i ( t ) − qcond ,i ( t ) + qconv ,i ( t ) − J j ( t ) Fi → j ′′
j =1 N DS
(10.103)
= ρi c pi δ i
All the terms on the right hand side of eq. (10.103) are for the design surface elements and are therefore known, (the Jj terms on the design surface are known when Tj and q"j are prescribed); the terms on the left hand side are for the heaters or other energy sources and are to be determined. Equation (10.103) is thus a discretized version of an ill-posed Fredholm integral equation of the first kind, and a generalization of the relation used in Example 10.11. The various techniques for solving design problems of this type can be grouped under the three broad headings of regularization, optimization, and metaheuristics. Here, we will gain an understanding by studying one regularization method. Regularization techniques consist of modifying the governing relations to reduce their “ill-posedness,” accepting some loss of accuracy to gain a useful solution.
10.4.2 Regularization
A direct or explicit solution to the inverse design problem requires use of an inverse formulation. As illustrated in Example 10.12 and Fig. 10.22, ordinary techniques (e.g., Gauss-Seidel, Gauss elimination or LU decomposition) are likely to either identify non-physical solutions with large amplitude fluctuations and/or complex absolute temperatures, or completely fail to find a solution. To achieve an accurate and physically reasonable solution, the explicit system of equations may be regularized by modifying the ill-conditioned set to a nearly-equivalent set of well-conditioned equations. The solution is then subject to some error and the level of regularization must be selected so that the accuracy of the solution satisfies the designer’s needs. One regularization technique is truncated singular value decomposition (TSVD). TSVD is based on modifying the singular value decomposition of A. The solution uses the pseudo-inverse matrix that is formed by filtering or truncating some (or many) of the smallest singular values, thus reducing the condition number of the matrix A; the solution to Ax = C using the p largest singular values from the singular value decomposition of A, A = USV T , becomes
Chapter 10 Heat Transfer by Radiation 835
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
xn = Vn , k
k =1
p
C
m
N
m
T U m ,k
Sk ,k
, n = 1, 2, , N
(10.104)
where p has a value less than or equal to the rank of A (Hansen, 1998). Retaining different numbers of singular values yields alternative solutions. Those with acceptable accuracy and physically reasonable characteristics present allowable alternatives. For the case analyzed in Example 10.11 with singular values shown in Fig. 10.23, the solution can be regularized by retaining some small subset of the 30 singular values. The predicted heater emissive power distributions found by retaining 1, 3 or 7 singular values are shown in Fig. 10.24 This illustrates that multiple possible physically reasonable solutions may be generated. The emissive power distribution for p = 7 contains some negative θ 42,j values (which cannot be realized in a physical system), and the results for larger p values become even more unrealistic. The accuracy of the three solutions in predicting the required q1"(x1) are compared in Fig. 10.25. The solution shown in Example 10.11, which results from retaining all 30 singular values of the matrix A, will produce a prediction for Q"1 that exactly matches the desired Q"1(X1) curve, but as mentioned in the example, requires heater inputs that cannot be physically realized. The rms error in the solutions is ′′ Q1′′given ( j ) − Q1,calculated ( j ) , and has values for p = 1,3,7 of 0.163, 0.0884, and , 2 0.0397. Thus, the error becomes smaller as more singular values are retained, even though the results tend to become more non-physical.
40
eb1TSVDi , 1 20 θ42,i, p=1 θ42,i p=3 eb1TSVDi , 3 θ42,i p=7 eb1TSVDi , 7 0
20
0
0.5
1
X2(i) Figure 10.24 Regularized TSVD Solutions to Example 10.11.
836 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
4
6 q2TSVD Q"1,i, p = 1 j , 1 Q"1,i p = 3 q2TSVD Q" p = 7 j , 3
1,i
8
q2TSVD j , 7 10
12
0
0.2
0.4
X1( j) 2 (j)
0.6
0.8
Figure 10.25 Predicted Heat Flux on Surface 1 Using Various TSVD Solutions
A problem that incorporates many attributes of a practical design process for a radiant furnace was proposed and solved by a diverse team to compare various solution methods for inverse problems (Daun et al., 2006). A three-dimensional geometry was considered, and the problem was to determine the radiant heater settings that provide a prescribed transient but spatially-uniform heating of the design surface.
10.4.3 Unresolved Problems in Inverse Cases
Although a considerable amount of progress has been made in developing methods for solution of inverse problems involving radiating systems, there remains significant work to be done. In specifying the conditions on the design surface, for example, the design engineer may choose conditions for which there is no acceptable physical solution, i.e. the solution might not be achievable without unacceptable heater conditions (coolers in place of heaters, excessive heater power or temperature requirements, or even imaginary absolute temperatures on the heaters). Such solutions may satisfy the conditions on the design surface mathematically, but they are not useful engineering solutions. An a priori determination of the existence of acceptable physical solutions would save a lot of fruitless calculation, but a means to do this has not yet been developed. The predictions of the inverse solution will not be exact, since assumptions are invariably built into the forward solution which is being inverted. Such
Chapter 10 Heat Transfer by Radiation 837
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
assumptions may include approximated surface properties (diffuse, gray, specular, etc.), and approximated thermophysical properties (conductivity, specific heat, etc.). Thus, to achieve specified conditions on the design surface in a real system such as a radiant furnace, feedback control is probably necessary. For feedback control, at least two additional factors must be addressed: 1) how feedback of measured temperature and/or radiative flux from the design surface can be used to adjust the inverse predictions and 2) how many feedback signals and their locations on the design surface are necessary to adequately provide feedback information. The first of these has been addressed in a preliminary way (Erturk et al., 2002) through the use of neural nets.
10.5 The Effect of Participating Media
The investigation of radiative transfer so far has contained the assumption that the radiative transfer between and among surfaces is unaffected by any medium that exists between the surfaces. Often this is a good assumption; however, in important cases the radiation is attenuated on its path between an emitting and an absorbing surface. In addition, the medium between the surfaces may emit energy, adding to the radiation incident on enclosing surfaces. The effect of a participating medium is especially important in combustion systems, where combustion products and ash and soot particles cause significant attenuation within the medium as well as emission from high-temperature flames and gas mixtures. Other systems where medium effects are important are in shock-layer emission to the surface of reentry heat shields for spacecraft; attenuation by interstellar dust clouds that affect observations; emission from the hot gases at the surface of the sun; the damage to lenses in optical systems due to even very small absorption of energy from a high-intensity laser beam; absorption in the Earth's atmosphere by the "greenhouse gases" that affects global warming; cancer detection by monitoring of skin-temperature anomalies due to submerged tumors; and attenuation of laser energy by the skin and organs that affects everything from treatment of cataracts in the eye to tattoo removal.
10.5.1 Absorption, Emission and Scattering from a Medium
The intensity I was defined and discussed in Section 9.2.1. A property of the intensity that was not discussed is found by looking at the diagram in Fig. 10.26, where we can calculate the energy from the two equal area elements dA1 and dA2 that is incident on surface dAs. The intensity definition, Eq. (9.5), is used to determine the energy striking dAs from the two equal area elements dA1 and dA2, resulting in dAs cos θ (10.105) deλ ,1 = I λ ,1dA1d Ω1dλ = I λ ,1dA1 dλ 2
S1
838 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
d A2
d Ω2 d A1
θ
d Ω1 d As S1
S2
Figure 10.26 Propagation of intensity
deλ , 2 = I λ ,2 dA2 d Ω 2 dλ = I λ ,2 dA2
dAs cos θ dλ 2 S2
(10.106)
Let dA1 = dA2 = dA, and impose the inverse square law for energy attenuation with distance, deλ ,2 =
S12 deλ ,1 . Making these substitutions results in the intensity 2 S2
from the two surfaces that reach dAs being
(10.107) This shows that, in the absence of attenuation or emission by a medium along the path S, the intensity is constant with distance. The intensity thus serves as a metric for determining the effect of attenuation.
I λ ,1 = I λ ,2
10.5.2 Properties of Participating Media
Absorption is the conversion of the radiation into internal energy of the medium; scattering is the change in direction of the radiation away from the direction of propagation without loss of energy (radiation scattering in cases of engineering interest are elastic scattering events); and emission is radiation originating from the medium that adds to the intensity. Consider intensity Iλ that is propagating through a medium in a particular direction (Fig. 10.27). The spectral dependence is maintained in the relations being developed, because the properties of many media, particularly gases, are highly wavelength dependent. As the intensity moves through an element of
Chapter 10 Heat Transfer by Radiation 839
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Iλ(S)
Iλ(S+dS)
dS
Figure 10.27 Attenuation of intensity in element of thickness dS
thickness dS, it is attenuated by absorption into the medium, and by scattering from the medium. The change in intensity is given by I λ ( S + dS ) = I λ ( S ) + dI λ ,attenuation = I λ ( S ) − β λ I λ ( S )dS (10.108) so that dI λ ,attenuation = − β λ I λ ,attenuation ( S )dS (10.109) The loss in intensity dI λ ,attenuation in the element dS is assumed to be proportional to the local intensity Iλ(S), the thickness of the element dS, and a property of the element called the attenuation coefficient, βλ.. The attenuation coefficient can also be shown to be the inverse of the mean free path for radiation through the medium; i.e., a large value of βλ indicates large attenuation and a thus short mean free path. This derivation is based on an intuitive argument, but the results are confirmed by experiment, and can be rigorously derived using the electromagnetic theory relations for attenuating materials. Examination of eq. (10.109) shows that βλ must have units of inverse length, and is usually in (1/m). Because attenuation is caused by both absorption in the element dS and scattering from that element, the attenuation coefficient is composed of an absorption coefficient, κλ, and a scattering coefficient, σs,λ. (The s subscript on σsλ is used to avoid possible confusion with the Stefan-Boltzmann constant.) Both have units of inverse length, and β λ = κ λ + σ sλ (10.110) If the intensity change due to attenuation over a finite distance is needed, then eq. (10.109) can be rearranged and integrated over a path from S = 0 to a final distance S:
The result is
Iλ ( s )
dI λ ,attenuation Iλ
I λ ( S = 0)
= −
S
S =0
β λ ( S )dS
(10.111)
840 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Iλ ( s )
dI λ ,attenuation Iλ
I λ ( S = 0)
= ln I λ ( S ) − ln I λ ( S = 0)
S I (S ) = ln λ = − S = 0 β λ ( S )dS I λ ( S = 0)
(10.112)
Raising both sides to the exponential power gives
I λ ( S ) = I λ ( S = 0) exp[−
S S =0
β λ ( S )dS ]
(10.113)
Transmittance: If the attenuation coefficient is uniform along S (which is assumed to be the case in this chapter; if it's not, then calculation of radiative transfer in other than the simplest geometries becomes quite computationally intense), then eq. (10.113) reduces to I λ ( S ) = I λ ( S = 0) exp[− β λ S ] (10.114) The intensity is seen to decrease exponentially as it travels along the path. We can define the transmittance of the medium τλ(S) as
τ λ ( S ) = exp[−
S S =0
β λ ( S )dS ] = exp ( − β λ S ) [if β λ ≠ F ( S )] (10.115)
The transmittance depends on the thickness of the medium and the value of the attenuation coefficient, and gives the ratio of the transmitted to the incident intensity. It should not be confused with the transmissivity, which is a property of an interface between two media, and is independent of the thickness. Absorptance: If the attenuation is due only to absorption with no scattering, then the fraction of intensity absorbed plus that transmitted must sum to unity. The absorptance of the medium is then defined as (10.116) α λ ( S ) = 1 − τ λ ( S ) = 1 − exp ( −κ λ S ) [if κ λ ≠ F ( S )] Again, this should not be confused with the surface property absorptivity, αλ, (Section 9.3.1) which is independent of the thickness S. Emittance: Using a similar argument to that for Kirchhoff's Law for surface properties, the ability of an absorbing medium to emit radiation can be used to relate the absorptance to the property emittance ελ(S): ε λ ( S ) = α λ ( S ) = 1 − τ λ ( S ) = 1 − exp ( −κ λ S ) [if κ λ ≠ F ( S )] (10.117) The emittance is used to predict the intensity traveling in a given direction due to emission from an isothermal medium along a path of length S: I λ ( S ) = ε λ ( S ) I λ b ( T ) = [1 − exp(−κ λ S ) ] I λb ( T ) (10.118) and for a differential element of medium dS,
dI λ ,emitted = I λb ( T ) d [1 − exp(−κ λ S )] = κ λ I λb (T ) dS dS
(10.119)
Now consider an isothermal volume element dV (Fig. 10.27). The intensity leaving the small element of length dS and face area dAs is given by eq. (10.119). Adding the intensity emitted by all elements parallel to the one shown in Fig. 10.27 gives the average intensity emitted in that direction by dV, or
Chapter 10 Heat Transfer by Radiation 841
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
dV dAp T dS
dIλ,emitted
dAs
Figure 10.28 Isothermal emitting volume element
dI λ ,emitted dAp =
dAp
κ λ I λb ( T )dSdAs = κ λ I λb ( T )
dAp
dSdAs = κ λ I λb ( T ) dV
(10.120) where the final integral over all parallel volume elements is evaluated by noting that adding the volumes of all of the small elements of length dS and face area dAs (as shown in Fig. 10.28) gives the total volume dV. The factor dAp is the projected area of the entire volume element in the direction of the intensity. To get the total energy emitted by the volume element, we now use the definition of the intensity [Eq. (9.5)] to find the energy emitted in one direction as d 3 eλ = dI λ ,emitted ( T ) dAp dΩd λ = κ λ I λb dVd Ωdλ (10.121) where the cosθ term in Eq. (9.5) is omitted because dAp is already the projected area. To get the emission into all directions, integrate over the 4π of solid angles: d 2 eλ = κ λ I λb ( T ) dVdλ dΩ = 4πκ λ I λ b ( T ) dVd λ = 4κ λ E λb dVd λ (10.122) 4π Finally, we can integrate over all wavelengths to get the total emission from the volume element:
de = 4dV
∞
λ =0
κ λ E λb d λ
(10.123)
If the spectral dependence of κλ is known, then the integral in eq. (10.123) can be evaluated, and we can define a mean absorption coefficient (the Planck mean) as (10.124) σT 4 and eq. (10.123) reduces to an equation for total emission into all directions from d V: de = 4κ P σ T 4 dV (10.125)
=0
κP
=λ
∞
κ λ E λb d λ
842 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
dIλ,s(θ,φ)
φ
d Ωi
d Ωs
Iλ(S)
θ
Forward
dS
Figure 10.29 Scattering from a volume element
Scattering: If radiation is scattered from the direction S into another direction (θ, φ) defined relative to the direction of S (Fig. 10.29) information is needed on how much of the scattered intensity goes into a particular direction. A function giving this information is defined as the scattering phase function, Φ, 4π dI λ , s (θ , φ ) Φ λ (θ , φ ) = (10.126) σ sλ I ( S )dS The numerator in eq. (10.126) is 4π times the intensity scattered into a given direction, and the denominator is the scattered energy (into all directions) from the intensity passing through dS. The intensity scattered into direction (θ,φ) is Φ (θ , φ )σ sλ I ( S )dS dI λ , s (θ , φ ) = λ (10.127) 4π Most scattering systems have randomly oriented particles such as soot or ash, and the scattering coefficient σsλ is assumed to be independent of direction through the medium. There are exceptions such as radiation incident on the carbon filaments held in a transparent epoxy matrix in a filament wound structure such as a golf club shaft; in that case both σsλ and Φ depend on the angle of incidence onto the filaments. The shape of the distribution of scattered energy is usually quite complex, consisting of strongly peaked lobes when the particles and wavelength are close in magnitude. For spherical particles, electromagnetic theory was used by Gustav Mie (pronounced "me") (Mie, 1908) to predict the shape of Φ, which is strongly dependent on both the refractive index of the particles and the parameter ς = π D / λ , where D is the spherical particle diameter. Many approximations are used for Φ to avoid the complication of following the path of
Chapter 10 Heat Transfer by Radiation 843
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
multiply scattered radiation in radiative transfer calculations. The simplest assumption is that the scattered intensity is isotropic, in which case Φ = 1. We have now defined many of the important properties that must be considered in treating radiative transfer in the presence of a medium that affects radiative transfer; a participating medium. Now we can examine the calculation of radiative transfer considering the participating medium effects.
10.5.3 The Radiative Transfer Equation
The fundamental equation that describes the change in intensity at a local position S is found by combining the relations found in Section 10.5.2. The change in intensity is given by Change in intensity = - loss by attenuation (absorption + scattering) + gain by emission + gain by inscattering into direction S from other directions Using the mathematical forms for each term,
dI λ ( S ) = − β λ I λ ( S )dS + κ λ I λ b ( S )dS +
σ sλ dS 4π I λ (θ , φ )Φ (θ , φ ) dΩ i 4π Ω = 0
i
(10.128)
For isoptropic scattering (Φ = 1) and using dΩ = sinθ dθ dφ, dI λ ( S ) = − β λ I λ ( S )dS σ dS 4π + κ λ I λb ( S )dS + sλ I λ (θ , φ )dΩ i 4π Ω = 0 (10.129) = − β λ I λ ( S )dS σ dS π 2π + κ λ I λb ( S )dS + sλ I λ (θ , φ ) sin θ dθ dφ 4π θ = 0 φ = 0 Dividing through by β λ dS ≡ dτ λ and defining ωλ = σ s ,λ / β λ , eq. (10.129) becomes dI λ (τ λ ) ω π 2π = − I λ (τ λ ) + (1 − ωλ ) I λb (τ λ ) + λ I λ (θ , φ ) sin θ dθ dφ dκ λ 4π θ = 0 φ = 0 (10.130) The ωλ is called the single scattering albedo, and is a measure of the importance of scattering in the overall attenuation in the medium, and dτλ is the differential optical thickness or opacity, a dimensionless path length that is a measure of the
i
844 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
attenuation along a path. (This is not the same quantity as the transmittance, τ(S), although the same symbol is commonly used for both.) Often, the form of the RTE is further simplified by defining the source function iλ (τ λ ) , which is the sum of the sources of intensity due to emission and in-scattering, ω π 2π iλ (τ λ ) ≡ (1 − ωλ ) I λb (τ λ ) + λ I λ (θ , φ ) sin θ dθ dφ (10.131) 4π θ = 0 φ = 0 so that eq. (10.130) becomes dI λ (τ λ ) = iλ (τ λ ) − I λ (τ λ ) (10.132) dτ λ This is called the differential form of the radiative transfer equation, or RTE. Interpreting eq. (10.132) shows that the change in local intensity with optical thickness is due to the gain from the source function less the loss by attenuation. The RTE is a first-order integro-differential equation (with the integral buried in the source function). It can be formally integrated by multiplying through by the integrating factor exp(τλ), giving dI (τ ) (10.133) I λ (τ λ ) exp(τ λ ) + λ λ exp(τ λ ) = iλ (τ λ ) exp(τ λ ) dτ λ The left-hand side is equal to
d I λ (τ λ ) exp (τ λ ) , and integrating with respect dτ λ
to optical thickness from τλ = 0 to τλ (S) results in τλ I λ (τ λ ) = I λ (0) exp ( −τ λ ) + i (τ *λ ) exp − (τ λ − τ *λ ) dτ *λ τ * =0 λ
λ
(10.134)
where τ λ ( S ) = S *= 0 β λ ( S *) dS * . If the properties are invariant across the medium, τ λ ( S ) = β λ S .
S
τλ*
Iλ(0) d τλ∗
τλ - τλ *
Iλ(τλ)
τλ
Figure 10.30 Propagation of intensity
Chapter 10 Heat Transfer by Radiation 845
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Equation (10.134) is the integrated form of the RTE. It shows (Fig. 10.30) that the intensity at a position τλ along a path is due to intensity leaving the boundary at τλ =0 and exponentially attenuated, plus added (but attenuated) intensity from the source function integrated along the entire path. Equation (10.134) is the fundamental relation for treating radiative heat transfer in a participating medium. If the medium is non-scattering and nonabsorbing, then examination of eq. (10.132) shows that dIλ = 0, and the intensity is constant as expected. Example 10.12: A black surface at T = 1000 K bounds a nonscattering medium with absorption coefficient κλ = 20 cm-1. The medium is also at T = 1000K (Fig. 10.31). Find the intensity in the medium at a distance of S = 20 cm from the boundary at wavelength λ. Solution: For the black boundary, I λ ( S = 0) = I λb (T = 1000 K ) . The source function [eq. (10.131)] is, for no scattering (ω=0) iλ (τ λ ) = I λ b (τ λ ) = constant = I λb (T = 1000 K ) and the RTE [eq. (10.134)] gives I λ (τ λ ) = I λ (0) exp ( −κ λ S )
+
S
= I λb (T = 1000 K ) exp ( −κ λ S ) −1 + κ λ I λb (T = 1000 K ) ( exp ( −κ λ S ) − 1) κλ = I λb (T = 1000 K )
S *= 0
κ λ I λb (T = 1000 K ) exp −κ λ ( S − S *) dS *
The intensity in this particular case is independent of S and the value of
κλ. This occurs because the attenuation of the intensity leaving the
boundary is exactly compensated by the emission along the path.
κλ = 0.2 cm-1
T = 1000K T = 1000 K Iλb(T) S=20 cm Iλ(S)
Figure 10.31 Propagation of intensity (T = 1000 K)
846 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
n
θ
s
dA dΩ Iλ
Figure 10.32 Finding the radiative flux from the intensity.
Because this is a heat transfer book, the next step is to relate the intensity to the radiative heat flux. This is done by integrating the intensity passing through a plane over all incident angles as shown in Fig. 10.32. The flux in the n direction is
′′ qrad , λ =
4π Ω= 0
I λ sndΩ =
4π
Ω= 0
I λ (θ ) cos θ d Ω
(10.135)
The integration over 4π solid angles accounts for radiation crossing dA in both negative and positive directions, and qr′′ad is the net positive radiative heat flux crossing dA in the positive x direction in the wavelength range dλ at λ. The total radiative heat flux is then found by integrating eq. (10.135) over all wavelengths. For use in an energy equation written on a volume element of the participating medium such as Eq. (3.89), the total (all wavelengths) radiative flux divergence ∇q′′ad is needed. This is found by a radiative balance similar to that r used in derivations of the mass, energy, or momentum equation. Considering the ′′ differences in radiative flux qrad ,λ crossing each face of a volume element in a Cartesian coordinate system gives
∇q′′ = rad ′′ ∂qrad , x ∂x + ′′ ∂qrad , y ∂y + ′′ ∂qrad , z ∂z
(10.136)
Reflecting now on what is needed to solve a complete radiative transfer problem, then, at every volume element in the medium, the total radiative flux divergence ∇q′′ is needed. To find this, the total flux in each coordinate rad ′′ ,i is needed, and this requires knowledge of the spectral intensity in direction qrad each direction passing through the local volume element. That in turn requires evaluation of the RTE at the location of the element at every wavelength and every direction. But solution of the RTE depends on knowledge of the entire temperature field throughout the participating medium. It is perhaps clear that a
Chapter 10 Heat Transfer by Radiation 847
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
full solution of radiative transfer is computationally expensive, and often is the most time-consuming part of solving a mixed-mode heat transfer problem in a participating medium. Approaches in general use for solution of the integral RTE include the zonal method (subdividing the medium into isothermal volumes and boundary surfaces, and pre-calculating factors for the exchange among them that are similar to configuration factors, but that account for attenuation and emission by the medium), discrete ordinates (subdividing the 4π of solid angles into a limited number of discrete solid angle increments, and writing the RTE for each of these directions for each volume element), discrete transfer (somewhat similar to the discrete ordinates method, but integrating along a path through each intervening volume element), Monte Carlo (treating radiation as composed of energy "packets" that each follow the laws of emission, absorption, scattering, and surface interactions, and then following a large number of these packets through individual histories to find their ultimate fate), and direct FEM solution methods. For the differential form of the RTE, the Pn method (using spherical harmonics expansion of the intensity) has proven useful. For details on these and other general solution methods, refer to Siegel and Howell (2002) or Modest (2003). Because the complexity of the general case places these more complete solution methods outside the scope of this book, some useful limiting cases will be examined.
10.5.4 Some Limiting Solutions for Radiative Transfer
The Diffusion Solution: If the participating medium is optically thick, then the solution of the RTE reduces considerably in form. Optically thick means that the mean free path for radiation, 1/ (κλ + σs,λ) = 1/βλ, is much smaller than the dimensions of the system being studied. In that case, radiative transfer becomes a diffusion process. This is analogous to the reduction of a molecular energy transport analysis to a conduction (diffusion) analysis at high molecular densities. Consider a 1-dimensional system of infinite parallel plates containing a participating, isotropically scattering medium with wavelength dependent attenuation coefficient βλ (see Fig. 10.33).
T2
x
θ
S
H
T1
Figure 10.33 Geometry for diffusion solution
848 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
If the plates are separated by a distance H, the RTE, eq. (10.132), reduces in this one-dimensional case to cos θ dI λ ( x , θ ) cos θ dI λ ( x , θ ) = = iλ ( x / H ) − I λ ( x / H , θ ) (10.137) dx H βλ d( x / H ) βλ If 1/Hβλ << 1 (optically thick), then the intensity can be expanded in a series,
I λ ( x ,θ ) = I λ
( 0)
1 ( 2) 1 I λ (1) + + I λ + ...... βλ H βλ H
2
(10.138)
Substituting the expansion into eq. (10.137) and using eq. (10.131) for the source function iλ(x/H), (0 (1 cos θ ∂I λ ) ( x , θ ) 1 ∂I λ ) ( x , θ ) + + ... H β λ ∂( x / H ) H βλ ∂( x / H )
ω = (1 − ωλ ) I λb (τ λ ) + λ 4π
2 1 (1) 1 ( 2) I λ ( 0) + Iλ + I λ + ......dΩ i θ =0 φ =0 βλ H βλ H
π
2π
(10.139)
2 0 1 ( 2) 1 1 − Iλ( ) + Iλ( ) + I λ + ...... βλ H βλ H
Collecting all terms of zero order in 1/ H β λ gives ω 4π (0) (0 I λ ) = (1 − ωλ ) I λb (τ λ ) + λ I dΩ i (10.140) 4π Ω = o λ ( Because neither term on the right-hand side depends on Ωi, it follows that I λ0) also cannot depend on Ωi. In that case, the integral in eq. (10.140) can be evaluated, and the result is ( I λ0) = I λb (τ λ ) (10.141) This identifies the first term in the series in eq. (10.138). Now the first order terms are collected in eqs. (10.139), giving ∂I ( x ) ωλ 4π (1) (1) cos θ λb (10.142) = dΩ =0 I λ dΩi − I λ ∂ ( x / H ) 4π
1
i
Multiply through by dΩi = 2π sin θ dθ and integrate over dΩi. The result is
0=
4π Ωi = 0
( I λ1) d Ω i − ωλ
4π
Ωi = 0
( I λ1) dΩ i
(10.143)
Because ωλ ≠0, the integral in eq. (10.143) must be equal to zero. In that case, eq. (10.142) reduces to
∂I ( x ) 1 I λ ( ) = − cos θ λb ∂( x / H )
(10.144)
Now, the first two terms in the series for intensity [eq. (10.138)] are known. If the series is truncated after these two terms, the result is
I λ ( x ,θ ) = I λb ( x ) − 1
βλ
cos θ
∂I λb ( x ) ∂x
(10.145)
Chapter 10 Heat Transfer by Radiation 849
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
For the diffusion case, then, the local intensity depends only on the local blackbody intensity and its gradient. Now the local spectral radiative flux is found by substituting eq. (10.145) into eq. (10.135):
′′ qrad ,λ ( x ) =
∂I ( x ) 1 cos θ λb I λb − cos θ dΩ i βλ ∂x π /2 π /2 2π 2π 1 ∂I ( x ) 2 λb = ϕ =0 I λb cos θ sin θ dθ dφ − θ =−π / 2 φ =0 βλ ∂x cos θ sin θ dθ dφ θ =−π / 2
4π
Ωi = 0
I λ (θ ) cos θ dΩ =
4π
Ωi = 0
= 0−
4π ∂I λb ( x ) 4 ∂E λb ( x ) =− ∂x 3β λ ∂x 3β λ
(10.146 )
For the diffusion solution, the local spectral radiant heat flux is thus dependent on only the local gradient in the spectral emissive power. This is known as the Rosseland diffusion equation, and was first derived from studies of radiative transfer through the Sun's upper atmosphere. If the properties of the participating medium are gray, then integrating over all wavelengths gives 4 dEb ( x ) 4σ dT 4 ( x ) 16σ T 3 dT ( x ) ′′ =− =− (10.147) qrad = − 3β dx 3β dx 3β dx This is the same form as the Fourier equation for heat conduction, and the radiative conductivity can be defined as 16σ T 3 (10.148) krad = 3β If the medium is non-gray, then eq. (10.146) can be integrated over wavelength to find the total heat flux:
′′ qrad = − 4 ∞ 1 ∂Eλb ( x ) 4 dEb λ =0 β λ ∂x dλ ≡ − 3β R dx 3
(10.149)
Here, βR is an averaged attenuation coefficient known as the Rosseland mean attenuation coefficient. If attenuation is due only to absorption (nonscattering medium), the value of κR is found from
1 ∂E λb ( x ) dλ ∂x = dEb κR dx ∞ 1 dE λb 1 ∞ 1 dEλ b = dλ = dλ λ =0 κ σ λ = 0 κ λ dT 4 λ dE b 1
λ
∞
=0
κλ
(10.150)
The derivative dEλb / dT 4 is found by taking the required derivative of the Planck blackbody relation with respect to T4. The use of the Rosseland mean absorption coefficient has some obvious problems; if there are regions of the spectrum that are transparent, then eq. (10.150) will become indeterminate. This is overcome in practice by finding the Rosseland mean for those parts of the spectrum that meet the optically thick criterion, and applying the diffusion solution in those regions.
850 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Where the medium is transparent, the wavelength dependent solutions for no participating medium of Section 10.2.3 are applied. Example 10.13: A slab of plate glass has thickness L = 1 cm. The temperature on one face of the glass is measured to be 800K, and the other face is at 300 K. The thermal conductivity of this glass is k = 1.4 W/m-K, and the absorption coefficient of the this glass over the important wavelength range is κ = 2.3 cm-1. What is the temperature distribution in the glass? (Assume the glass is non-scattering.) Solution: The energy equation within the glass is 16σ k∇ 2T + ∇(krad ∇T ) = k∇ 2T + ∇(T 3∇T ) = 0 3κ or, in one dimension d 2T 16σ d dT (T 3 )=0 k 2+ 3κ dx dx dx Integrating once gives dT 16σ T 3 dT dT + = (k + krad ) = C1 k dx 3κ dx dx or 16σ T 3 (k + )dT = C1dx 3κ Integrating a second time between T1 = 300K at x = 0 and T2 = 800 K at x = L gives 16σ 4 k (T2 − T1 ) + (T2 − T14 ) = C1 L 3κ
Figure 10.34 Temperature distribution
Chapter 10 Heat Transfer by Radiation 851
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
This allows evaluation of C1, and then dividing through by the thermal conductivity T (x) 16σ T 3 16σ 4 4 (1 + )dT = T ( x ) − 300 + 3κ k T ( x ) − 300 T =300 K 3κ k 16σ 4 = (T2 − T1 ) + (T2 − T14 ) ( x / L ) 3κ k Substituting numerical values
16 × 5.67 × 10−8 (W / m 2 ⋅ K 4 ) 16σ = = 9.39 × 10−10 (1/ K 3 ) gives 3κ k 3 × 2.3(cm −1 ) × 1.4(W / m ⋅ K ) × 100(cm / m)
[T ( x ) − 300] + 9.39 × 10−10 (1 / K 3 )[T 4 ( x ) − 3004 ] = 877( x / L )
Solving for the temperature distribution (the easy way is to pick T values between 300 and 800K and find the corresponding value of x/L) gives the graph shown in Fig. 10.34. The temperature is not linear between the boundaries because of the temperature-dependent radiative conductivity. Also, because this is a one-dimensional steady state problem, the total heat flux at any value of x/L must be constant (and examining the energy equation after one integration shows that it's equal to C1); however the relative amount of radiative and conductive heat flux varies across the plate. The Slip Boundary Conditions. In the absence of conduction or convection, radiation problems have unfamiliar behavior near boundaries. The medium temperature does not approach the boundary temperature as a solid boundary with fixed temperature is approached. This can be seen for the case of a gray medium near a gray boundary (Fig. 10.35). An energy equation can be written on a plane very near the boundary. If the boundary is gray-diffuse and at temperature Tw, then the radiative flux 4 ′′ ′′ (radiosity) leaving the boundary in the +x direction is q+ x = εσ Tw + (1 − ε )q− x . q"-x
q"+x
Figure 10.35 Derivation of slip boundary condition
852 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
The flux in the negative x direction comes from the optically thick medium and can be found by integrating eq. (10.145) over the hemisphere of directions with intensity having negative x components. This gives
′′ q− x = =
2π 2π
φ =0
θ π
=
−π / 2
I ( x , θ ) cos θ sin θ dθ dφ
φ =0
θ π
=
−π / 2
∂I b ( x = 0) 1 I b ( x = 0 ) − cos θ cos θ sin θ dθ dφ β ∂x
−π / 2
2π ∂I b ( x = 0) cos3 θ 2π ∂I b ( x = 0) = π Ib ( x = 0) − = Eb ( x = 0 ) + ∂x ∂x β 3β 3 θ =0 = Eb ( x = 0 ) + 2 ∂Eb ( x = 0) q′′( x = 0) = Eb ( x = 0 ) − ∂x 3β 2
(10.151)
Substituting into an energy balance across the area dA adjacent to the boundary, 4 4 ′′ ′′ ′′ ′′ ′′ q′′ ( x = 0 ) = q+ x − q− x = εσ Tw + (1 − ε ) q− x − q− x = ε (σ Tw − q− x )
4 q′′ ( x = 0 ) = ε σ Tw − σ T 4 ( x = 0) − 2
(10.152)
Substituting eq. (10.151),
4 T 4 ( x = 0) = Tw −
q′′ ( x = 0 ) 1 1 − σ ε 2
(10.153)
Equation (10.153) shows that the temperature of the medium at the boundary is less than the boundary temperature when the heat transfer is positive in the +x direction. Example 10.14 Find an expression for the radiative heat flux between two infinite parallel plates at T1, ε1 and T2, ε2 separated by a gray medium of optical thickness τ. Use the diffusion approximation. Solution: The diffusion solution, eq. (10.147) is 4σ dT 4 ( x ) 4σ dT 4 ( x ) ′′ =− qrad = − 3β dx 3 dτ For the one-dimensional case, the heat flux is constant, so the equation can be integrated between the medium boundaries to give τ 4σ T (τ ) 4σ 4 4 4 ′′ ′′ qrad dτ * = − T =T (τ =0) dT = − T (τ ) − T ( 0 ) = qradτ τ =0
3 3
Using the slip boundary condition eq. (10.153) to eliminate the medium temperatures at the boundaries, T4(0) and T4(τ) in terms of the boundary temperatures, 4σ ′′ T 4 ( 0 ) − T 4 (τ ) qrad = 3τ
= 4σ 3τ 4 qrad 1 1 4 qrad 1 1 ′′ ′′ T1 − − − T2 + − σ ε1 2 σ ε 2 2
Chapter 10 Heat Transfer by Radiation 853
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
or, finally,
σ ( T14 − T24 ) 3τ 1 1 + + −1 4 ε1 ε 2 Comparing with the result of Example 10.5 for no participating media, σ ( T14 − T24 ) ′′ = , it is seen that for the radiation diffusion case, the qrad
′′ qrad = 1
ε1 ε 2 presence of the participating medium adds a radiative resistance of 3τ / 4 .
+
1
−1
The Nearly Transparent Gas: If the medium has a small (but non-zero) absorption coeffficient, then it will still emit from a volume element in accordance with eq. (10.125). However, minimal attenuation due to absorption will occur as the emitted intensity travels along a path, and absorption by and within the volume element can be neglected. This observation is often used in analyzing radiation loss from nonscattering small flames such as for laboratory burners, since the radiation energy loss from the element becomes 4 qrad = 4κ P σ TdV dV (10.154) and this is easily included in an energy balance that includes the combustion energy release for determining flame temperature, and often serves as a firstorder radiative loss correction to the adiabatic flame temperature. Isothermal Media. Consider a high-temperature volume of participating medium contained in an enclosure with cold black walls. If the medium is well-mixed, then its temperature and properties are uniform within the enclosure. Given the temperature and composition of the medium, the Planck mean absorption coefficient can be computed from eq. (10.123) using the spectral absorption coefficient. Consider a hemisphere containing a nonscattering uniform isothermal medium at temperature T (Fig. 10.36). At the center of the base of the hemisphere is an element dA. The intensity striking dA from the direction of the ring element on the surface of the hemisphere is given by eq. (10.134): κλ I λ (κ λ R) = I λ (0) exp ( −κ λ R ) + κ i (κ λ r ) exp −κ λ ( R − r ) dr (10.155) κ * =0 λ λ
λ
Here, r is the distance along R from the hemisphere surface, and for this geometry and a nonscattering medium, τλ has been replaced with κλR and τλ * with κ λr. Since the hemisphere boundary is cold and black, Iλ(0) = 0. For the nonscattering isothermal medium, the source function iλ(κλr) reduces to [eq. (10.131)] iλ (κ λ r ) = I λb (κ λ r ) = I λ b (T ) (10.156)
854 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
θ
R dA
R
Figure 10.36 Black cold hemisphere enclosing a participating isothermal nonscattering medium with uniform absorption coefficient a.
and eq. (10.155) becomes
I λ (κ λ R) = I λ b ( T )
aλ R
=I λb ( T ) exp(−κ λ R) [ exp(κ λ R) − 1] = I λb ( T ) [1 − exp(−κ λ R) ]
aλ r = 0
κ λ exp −κ λ ( R − r ) dr
(10.157)
Note that the term containing the exponential is the spectral emittance ελ(R) of the medium (eq. (10.117). The heat flux on dA is now found from eq. (10.158)
′′ qrad ,λ =
4π
ω =0
I λ (θ ) cos θ d Ω =
4π
4π
ω =0 λ
ε ( R) I λb ( T ) cos θ dΩ
= ε λ ( R) I λ b ( T ) = ε λ ( R) I λ b ( T )
ω =0
2π
cos θ dΩ
φ =0 θ =0
π /2
(10.158)
cos θ sin θ dθ dφ = ε λ ( R)π I λ b ( T )
= ε λ ( R) Eλ b Equation (10.158) is clearly a very simple equation for determining the heat flux for this particular geometry. For more complex geometries, the quantity Iλ(θ) inside the integral will vary with Ω, and the integration to find the radiative flux at a location on the boundary will become quite complex. For any given geometry, however, there is some value of a fictitious length, called the mean beam length, Le, that can be substituted for R in the hemisphere relation, eq. (10.158), and will give the correct value for heat flux that is found from the detailed integration. The spectral radiative heat flux is then found from ′′ qrad ,λ = ε λ ( Le ) E λb (10.159) Integrating over all wavelengths and substituting the Planck mean emittance
ε ( L ,T ) = λ
e
∞ =0
ε λ ( Le ) Eλb ( T , λ ) dλ σT 4
gives the total heat flux as
′′ qrad = ε ( Le )σ T 4
(10.160)
Chapter 10 Heat Transfer by Radiation 855
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Table 10.1 Mean Beam Lengths for Common Geometries (Siegel and Howell, 2002) Geometry Characterizing Dimension Mean Beam Length Le for Finite Optical Thickness 0.65D 0.95D
Sphere radiating to surface Circular cylinder of infinite length radiating to interior surface Circular cylinder of height equal to two diameters radiating to: Plane end Concave surface Entire surface Infinite slab of medium radiating to: Element on one face Both bounding planes Rectangular parallelepiped: 1×1×4 volume radiating to: 1×4 face 1×1 face All faces
Diameter D Diameter D Diameter D
0.60D 0.76D 0.73D Slab thickness H 1.8H 1.8H Shortest side X 0.82X 0.71X 0.81X
Values for the mean beam length have been computed for many common geometries. Some of these are given in Table 10.1. General geometries: eq. (10.125) indicates that the total energy emitted by an isothermal volume element is de = 4aP σ T 4 dV (10.161) If the medium is optically thin,, then there is no attenuation of this energy before striking the cold black boundary, so the average radiative flux over the entire boundary is e 4κ σ T 4V ′′ (10.162) qrad = = P
A A
Also, for the optically thin case, the emittance of the gas can be rewritten using a series expansion of the exponential term to give 2 (κ Le ) ε ( Le ) = 1 − exp(−κ Le ) = 1 − 1 − κ Le + − ... ≈ κ Le ,0 (10.163)
2!
Substituting this result into eq. (10.160) gives ′′ qrad = κ Le ,0σ T 4 (10.164) Comparing eqs. (10.164) and (10.162) shows that for an optically thin medium radiating to a cold black boundary, the mean beam length is
Le ,0 = 4V A
(10.165)
By comparison of the result of eq. (10.165) with values of mean beam length computed using complete integration over the volume of a gas that is not optically thin, it is found that to correct the optically thin approximation for
856 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
absorbing gases, reasonable values for the mean beam length for most geometries are
Le = 0.9 Le ,0 = 3.6V A
(10.166)
Example 10.15 A participating medium is contained in a right circular cylinder of height equal to its base diameter, which is D = 4 m. The gas has emittance of ε ( Le ) =0.35+Le/20 (for Le < 10 m) at a temperature of 1200K. Find the average radiative heat flux on the cold black enclosure boundary. Find the mean beam length from Table 10.1, and compare the value with that from eq. (10.166). Solution: From Table 10.1, Le = 0.6D = 2.4 m. From eq. (10.166), 3.6V 3.6 × π D 2 L / 4 1.8 × DL = = = 2.4m Le = 2 A 2 (π D / 4 ) + π DL ( D ) + 2 L giving some confidence in the approximate relation of eq. (10.166). Then ε ( Le ) = 0.35 + 2.4 / 20 = 0.47 and
′′ qrad = ε ( Le )σ T 4 = 0.47 × 5.6704 × 10−8 (W / m 2 K 4 ) × 1200 4 ( K 4 ) = 55,300 W/m 2
Gas Emittance Values: Participating media have extremely varying wavelength dependence, and the assumption of a gray gas is much less applicable than the assumption of gray surfaces. The wavelength dependence occurs because gases absorb radiation through converting radiative energy into internal energy. For gases the conversion processes are by absorption of photons that causes energy transitions between quantized energy states. At lower gas temperatures, these transitions occur between vibrational-rotational states, and at higher temperatures electronic transitions and transitions from so-called bound electronic states to free electron states through dissociation and ionization can occur. Monatomic gases (argon, krypton, neon, etc.) and most homopolar diatomic molecules (N2, O2, etc) do not have transitions that are activated by radiation in the wavelength regions of engineering interest except at very high temperatures such as in spacecraft re-entry. Thus, the atmospheric gases are essentially transparent and do not participate in radiative transfer. Polyatomic molecules such as H2O, CO2, and hydrocarbons have transitions that can be activated (that is, absorption can occur) in the infrared. These gases are of particular engineering interest in the design of utility and chemical process furnaces, and in the study of global warming. The absorption of radiation results in a transition between a pair of quantized energy states Ei and Ej, and this appears as an absorption line in the absorption spectrum. The wavelength of absorption that causes this transition is given by
Chapter 10 Heat Transfer by Radiation 857
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
E j − E i = hν i − j =
hc
λi − j
(10.167)
where νi-j is the frequency of the absorbed photon and λi-j is its wavelength. The line is very narrow; in fact it would have zero spectral width if the lower and upper quantized energy states were separated by an exact energy interval. Absorption of a photon then would occur only for photons that have exactly the frequency or wavelength given by eq. (10.167). However, the energy difference between the upper and lower energy states is not exact due to several effects that perturb the values of the energy states in eq. (10.167) and cause the lines observed for absorption by a group of molecules to have finite spectral width. These include the uncertainty principal, the velocity distribution of the gas molecules relative to the speed of light (Doppler broadening), the interactions with nearby molecules that distort the energy states (collision broadening), the effects of electric fields from free electrons or related effects (Stark broadening) and others. Because there are very many vibrational-rotational states, there are many lines in the absorption spectrum. For CO2 and H2O, the lines are closely spaced in certain regions of the IR spectrum, and indeed the line widths overlap to the extent that in low resolution spectra, the lines coalesce into absorption bands (see Fig. 10.37). Given the computational difficulty of integrating the spectrum over thousands of individual lines within the bands to find total heat transfer (called a line-by-line calculation), much effort has been expended to develop higher-level engineering approximations to the spectral properties that can be used in practical calculations. Some approaches are based on finding averaged properties over
Band designation, wavelength
αλ(T,P,L)
Wavelength, λ (μm)
Figure 10.37 Low resolution absorption spectrum for CO2 gas at 830K, 10 atm., path length of 0.388 m (Adapted from Edwards, 1976)
858 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
narrow spectral bands based on either experimental measurement or line-by-line calculation (narrow band correlations); others are based on correlation of the characteristics of the wide absorption bands (wide-band approximations), and some correlate the properties over the entire spectrum to provide the total gas emittance. Finding useful approaches for correlating the absorbing properties of gases for use in engineering calculations is an active research area (Modest and Zhang, 2002; Denison and Webb, 1995). Because the emittance of a participating gas mixture (say CO2 and H2O as present in combustion products) depends on gas temperature, total pressure, and concentration, correlation is not simple. An additional complication for mixtures of participating gases is that the absorption lines and bands of the individual gases may overlap, and this must be accounted for in addition to predicting the properties of the individual gases. Hottel (1954) presented curves of the emittance ε(Le) for CO2 and H2O which have been widely used in the mean beam length approach. Hottel's curves were based on experimental data with extrapolations to some temperature and Lepartial pressure regions based on theory. These curves have been updated based on more recent data, and the approach remains very useful for the isothermal gas mean beam length approximation used in this section. Rather than using Hottel's original curves, it is useful to have analytical expressions for the emittance. For water vapor in air, Cess and Lian (1976) give the correlation ε H O ( Le , T ) = a0 1 − exp(−a1 X (10.168)
2
where X = pH O Le ( pair + bpH O ) ( 300 / T ) and b = 5.0 ( 300 / T ) + 0.5 . In these
1/ 2
2 2
relations, T is in K, p in atm, and Le in m. The constants in eq. (10.168) are in Table 10.2.
Table 10.2 Constants for use in eq. (10.168)
T(K) 300 600 900 1200 1500
a0 0.683 0.674 0.700 0.673
a1 (m-1/2atm-1) 1.17 1.32 1.27 1.21
0.624
1.15
Leckner (1972) gives empirical expressions for the total emittance derived from expressions for narrow band behavior summed over the spectrum for both water vapor and CO2. In these correlations and equations, p is in bar, and Le is in cm. The most accurate expressions from Leckner agree within five percent to experimental data for T>400K. The correlation equation is
ε = exp a0 + ai [ log( pLe )]
i
M
(10.169)
i =1
Chapter 10 Heat Transfer by Radiation 859
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
where ai = c0 i + c ji (T / 1000 )
j =1
N
j
and the values of cji are in Table 10.3 for water
vapor and CO2. A plot of the emittance predicted by eq. (10.169) is in Fig. 10.38. Observe that the emittance increases with the pressure-path length product as expected. The trend with temperature is that emittance generally decreases with increasing temperature for water vapor. In contrast, CO2 tends to go through a peak in emittance at about 1200 K.
Table 10.3: Coefficients cji for Water Vapor and CO2 Emittance in eq. (10.169) i 0 1 2 0 1 2 3 c0i c1i c2i Water Vapor, T > 400 K, M=2, N=2 c3i c4i
-2.2118 0.85667 -0.10838 -3.9781 1.9326 -0.35366 -0.080181
-1.1987 0.035596 0.93048 -0.14391 -0.17156 0.045915 Carbon Dioxide, T > 400K, M=3, N=4 2.7353 -1.9822 0.31054 -3.5932 3.7247 -1.4535 0.61766 -0.84207 0.39859 0.31466 -0.19973 0.046532
1
0.015719 0.20132 -0.063356 -0.0033086
(pL,T) Emittance, ε
Η2Ο
0.1
0.01
0.001
0.0001 500
pL=0.05 pL=1 pL=2 pL=8 pL=24 pL=64
1000
1500
2000
2500
3000
Temperature, T(K)
Figure 10.38 Computed emittance of water vapor from eq. (10.169) using the coefficients of Table 10.3. Units of pL are bar-cm.
860 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
The individual emittances for H2O and CO2 in air must be modified when both gases are present in a mixture, which is commonly the case. This is because the individual spectral lines and absorption bands for the two gases overlap in some spectral regions, and simple addition of the emittance will overpredict the effect of the mixture. In some cases, a simple addition can predict a gas absorptance and emittance that is greater than unity at certain wavelengths. The overlap correction equation is ε ( pLe ) = ε H O ( pH O Le ) + ε CO ( pCO Le ) − Δε (10.170)
2 2 2 2
Hottel (1954) presents a graph of the approximate band overlap correction Δε. An empirical expression for the band overlap correction that is in good agreement with the Hottel chart (Leckner, 1972) valid for 1000 < T < 2200 K and all pressures is ζ 2.76 Δε = − 0.0089ζ 10.4 [ log10 ( pLe ) ] (10.171) 10.7 + 101ζ where ζ = pH O / ( pH O + pCO ) , p is in bars, and Le in cm.
2 2 2
The correlations and overlap corrections are for properties of the absorbing gases mixed with air at a total pressure of one atmosphere. If the total pressure differs considerably from one atm, then a pressure correction must be applied to the predicted one atm emittance because of increased pressure broadening of the individual lines that make up the bands that are summed to obtain the total emittance. Hottel (1954) has also presented graphs for this correction, and Leckner (1972) has provided an algebraic expression. Because most engineering equipment where radiation is important operates at near to one atm, the pressure correction is not given here. Example 10.16 The cylinder described in Example 10.15 is filled with a mixture of 15 volume percent of CO2, 20 percent H2O vapor, and the remainder air. The total pressure of the gas mixture is 1 atm, and the gas temperature is 1200 K. What is the average heat flux over the boundary? Solution: The mean beam length is the same as that found in Example 10.15, Le = 2.4 m = 240 cm. The partial pressures of the gases are equal to the mole fraction of each times the total pressure. The mole fraction in an ideal gas mixture is equal to the volume fraction, so the partial pressures are pCO2 = 0.15, pH2O = 0.20, and pair = 0.65 atm. For water vapor, the ai values are ai = c0 i + c ji (T / 1000 ) , giving a1 = -3.599, a2
j j =1 N
= 1.766, and a3 = -0.248. Using eq. (10.169) (remembering to convert the pressures to bars)
ε ( H 2 O ) = exp a0 + ai [ log( pLe ) ] = 0.266.
i
M
i =1
Chapter 10 Heat Transfer by Radiation 861
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
A similar calculation for CO2 gives ε(CO2) = 0.146. The overlap correction is calculated using ζ = pH O / ( pH O + pCO ) = 0.20 / ( 0.15 + 0.20 ) = 0.571 and
2 2 2
(
pH 2O + pCO2 Le = ( 0.15 + 0.20 ) (atm) × 1.01325(bar/atm) × 240(cm)
)
= 85.1(bar ⋅ atm). Substituting results in
ζ 2.76 Δε = − 0.0089ζ 10.4 [ log10 ( pLe ) ] = 0.051 10.7 + 101ζ
The total emittance of the gas mixture is then ε ( pLe ) = ε H O ( pH O Le ) + ε CO ( pCO Le ) − Δε
2 2 2 2
Now the average heat flux on the boundary is ′′ qrad = ε ( Le )σ T 4 = 0.360 × 5.6704 ×10−8 (W/m 2 ⋅ K 4 ) × (1200K) 4
= 4.24 × 10 4 (W/m 2 )
= 0.266 + 0.146 − 0.051 = 0.360
The mean beam length approximation for radiation in an enclosure of isothermal gases provides a straightforward tool for radiative transfer estimates. The restrictions to isothermal gases with cold black boundaries can be relaxed; however, these approximations are often met quite well in utility boilers and furnaces where combustion gases are highly turbulent and well mixed, giving near isothermal gas mixtures; the surfaces are covered with ash and soot, making them near-black; and the boundaries are water-walls at relatively low temperatures compared with the combustion gases.
10.6 Applications of Radiative Transfer
10.6.1 Radiation Measurement and Sensing: IR Cameras, Optical Pyrometers and Remote Sensing
Many devices based on radiative transfer are used for remote sensing of surface characteristics. These include IR cameras, optical pyrometers, satellite sensors and others. All view the radiosity leaving the observed bounding surface. If the surface is not nearly black, then the reflected component of radiosity may be a significant fraction of the surface emission, and the computed surface temperature may be in error. The range of spectral sensitivity of each of these devices is generally selected so that absorption and scattering by media between the detector and the measurement source can be neglected. Finding a spectral "window" is not always possible, in particular for satellite sensing systems, where the atmosphere has scattering particles that cause some
862 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
attenuation at almost all wavelengths, and fog, smog, and clouds may be present. Because the absorption/scattering properties along the path are not known, this is an application where inverse analysis to determine the attenuation properties of the atmosphere becomes critical. If the environment is at a low temperature compared with the object being measured by the pyrometer, then the radiosity is only dominated by emission from the object, and J = εσT4. The temperature can then be inferred from the pyrometer measurement if the emissivity of the object is known.
10.6.2 Atmospheric Phenomena Caused by Scattering
In Section 10.5, scattering was discussed for its effect on the intensity passing through a participating medium. Most discussion was for isotropic scattering. Scattering of solar radiation does occur in the Earth's atmosphere. Even on a clear day, solar radiation is scattered from very small aerosol particles and even from molecules in the atmosphere. When the parameter ς = π D / λ is small, the scattering is a limiting case of Mie scattering called Rayleigh scattering, after Lord Rayleigh (right) who derived the scattering relations using dimensional analysis.
Figure 10.39 Atmospheric scattering causing a blue sky.
Chapter 10 Heat Transfer by Radiation 863
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Figure 10.40 Atmospheric scattering resulting in a red sunset.
In Rayleigh scattering, the wavelength dependent scattering coefficient σsλ is found to be proportional to 1/λ4. The complete relation for nonabsorbing particles in air depends on the particle diameter D, refractive index of the scattering particles n, the number density of particles N (particles/m3) and ς through
σ s,λ =
8 π D2 ς 4 n2 −1 3 4 N n2 + 2
2
(10.172)
The wavelength dependence enters through ς . Remembering from Section 9.1 that the visible portion of the spectrum stretches from the blue portion at about λ = 0.4 μm to red at about 0.7 μm, we expect the atmosphere to scatter much more of the blue portion than the red in the visible portion of the spectrum. Thus after multiple scatterings over the long mean free paths present in the clear atmosphere (Fig. 10.39), what an observer on the ground sees when looking at the sky is multiply scattered solar radiation strongly biased to the shorter wavelengths… a blue sky! Similar reasoning is used to examine what happens to solar radiation in the atmosphere near sunset. Figure 10.40 indicates that near sunset, solar radiation takes a much longer path through the atmosphere, and more scattering occurs. The blue portion of the spectrum is lost through Rayleigh scattering along the path. What is left when viewing the sun is the remaining portion of the visible spectrum, strongly skewed toward the longer wavelengths (red). The effect is even stronger because the refractive index change of the atmosphere with altitude curves the sunlight around the Earth's curvature, making the pathlength even longer.
864 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Figure 10.41 Phase function for Rayleigh scattering.
The phase function Φ(θ) for Rayleigh scattering is given by the smooth function
Φ (θ ) = 3 (1 + cos2 θ ) 4
(10.173)
which has a lobed shape with scattering strongly into the forward direction, an equal amount backscattered, and small scattering into side directions relative to the direction of travel of the intensity (Fig. 10.41). If you look at angles near the sun from Earth around noon, the sky near the sun appears nearly white rather than blue. This is because you are observing solar radiation that has been near-forward scattered once or twice into angles near the solar direction, so the blue-scattering effects have not yet become dominant. You are observing forward scattered (white) light from across the visible spectrum.
10.6.3 Pollution, Greenhouse Gases and the Greenhouse Effect, Atmospheric Radiation and the Global Energy Balance
The band absorption properties of the chief absorbing components of the atmosphere (CO2 and H2O) are clustered in the infrared part of the spectrum (Fig. 10.37). Remember that monatomic (argon) and homopolar diatomic (nitrogen and oxygen) gases do not absorb in the parts of the spectrum important for global energy balances. Both CO2 and H2O have strong absorption bands in the IR, and their ability to absorb depends in part on their concentration. Other polyatomic gases (particularly methane) have similar characteristics. The discussion on shields in Section 10.2.2 shows why the greenhouse gases are important in the discussion of global warming. They act as a shield between the radiation emission from the Earth's surface and the heat sink of outer space,
Chapter 10 Heat Transfer by Radiation 865
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Figure 10.42 Greenhouse gases trapping emitted radiation from the Earth's surface
adding additional resistance to heat rejection by the Earth, without significantly affecting the shorter wavelength solar radiation absorbed by the Earth (Fig. 10.42). Increases in the concentration of either CO2 gas or H2O vapor (or both) tend to cause an increase in the overall equilibrium temperature of the Earth. The same effect occurs in greenhouses and solar collectors, where the glass cover transmits solar energy in the visible spectrum; however, glass is nearly opaque at IR wavelengths, so emission from the plant beds or collector plate cannot escape, and the greenhouse or solar collector temperature remains well above the outdoor ambient temperature. This is actually the origin of the descriptive term "greenhouse" in "greenhouse gases." The Earth's equilibrium temperature is affected by many factors. Those we have touched on in Chapters 9 and 10 are the concentrations of the polyatomic gases in the atmosphere (CO2, water vapor, methane and others) and their effect on trapping solar energy, and the changes in the Earth's average and local solar and IR absorptivity as affected by changing patterns of vegetation and ice fields on the Earth's surface. Other important factors include the coupling of surface temperatures with free convection patterns in the atmosphere and water bodies; the effect of the Earth's rotation; the CO2 concentration in the atmosphere as affected by man's activities, plant cycles, and solution of CO2 in the ocean; the effect of particulates in the atmosphere from dust storms, volcanic activity, and emissions from coal, wood and other fossil-fueled sources. All of these form a coupled system that must be considered in modeling global temperature patterns. The material in this text forms a good basis for understanding the complexity of the interactions.
References
Alifanov, O.M., 1994, Inverse Heat Transfer Problems, Springer, Berlin. Alifanov, O.M., Artyukhin, E.A., and Rumyantsev, S.V., 1995, Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Heat Transfer Problems, Begell House, New York, NY.
866 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Beck, J.V., Blackwell, B., and St. Clair, Jr., C.R., 1995, Inverse Heat Conduction: Ill-Posed Problems, Wiley-Interscience, New York, NY. Cess, R.D. and Lian, M.S., 1976, “A Simple Parameterization for the Water Vapor Emissivity,” J. Heat Transfer, Vol. 98, pp. 676-678. Daun, K.J., Erturk, H., Gamba, M. Hosseini Sarvari, M. and Howell, J.R., 2003, “The Use of Inverse Methods for the Design and Control of Radiant Sources,” JSME International Journal, Series B, 46, no. 4, 470-478, 2003. Daun, K.J. França, F., Larsen, M., Leduc, G., and Howell, J.R., 2006, “Comparison of Methods for Inverse Design of Radiant Enclosures,” J. Heat Transfer, Vol. 128, pp. 269-282. Denison, M.K. and Webb, B.W., 1995, “The Spectral-Line Weighted-Sum-ofGray-Gases Model for H2O/CO2 Mixtures,” J. Heat Transfer, Vol. 117, pp. 788798. Edwards, D.K., 1976, “Molecular Gas Band Radiation,” in T.F. Irvine and J.P. Hartnett (eds.), Advances in Heat Transfer, vol. 12, pp. 115-193, Academic Press, New York, NY. Ertürk, H., Ezekoye, O.A., and Howell, J.R., 2002, “The Use of Inverse Formulation in Design and Control of Transient Thermal Systems,” Heat Transfer 2002: Proc. International Heat Transfer Conf., Grenoble, pp. 729-734. França, F.H.R., Howell, J.R., Ezekoye, O.A., and Morales, J.C., 2002, “Inverse Design of Thermal Systems,” Chap. 1, Advances in Heat Transfer, J.P. Hartnett and T.F. Irvine, eds., Elsevier, Vol. 36, 1-110. Hansen, P.C., 1998, Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion, SIAM, Philadelphia, PA. Hogan, R.E. and Gartling, D.K., 2007, “Solution Strategies for Coupled Conduction/Radiation Problems,” Commun. Numer. Meth. Engng,, Vol. 23, pp. 523-542. Hottel, H.C, 1954, “Radiant Heat Transmission,” Chap. 4, in W.H. McAdams (Ed.), Heat Transmission, 3rd ed., McGraw-Hill, New York, NY. Howell, J.R., 2003, A Catalog of Radiation Configuration Factors, 2nd ed., http://www.me.utexas.edu/~howell/ Incropera, F.P, DeWitt, D.P., Bergman, T.L., and Lavine, A.S., 2007, Fundamentals of Heat and Mass Transfer, 6th ed. John Wiley & Sons, New York, NY. Leckner, B., 1972, “Spectral and Total Emissivity of Water Vapor and Carbon Dioxide,” Combustion and Flame, Vol. 19, pp. 33-48. Mie, G., 1908, “Beiträge zur Optik trüber Medien speziell kolloidaller Metallösungen,” Ann. Phys., Vol. 25, pp. 377-445.
Chapter 10 Heat Transfer by Radiation 867
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
Mills, A.F., 1998, Heat Transfer, 2nd ed., Prentice Hall, New York, NY. Modest, M.M., 2003, Radiation Heat Transfer, 2nd ed., Academic Press, New York, NY. Modest, M.M. and Zhang, H., 2002, “The Full-Spectrum Correlated-k Distribution for Thermal Radiation From Molecular Gas-Particulate Mixtures,” J. Heat Transfer, Vol. 124, pp. 30-38. Oppenheim, A.K., 1956, “Radiation Analysis by the Network Method, Trans. ASME, Vol. 27, pp. 725–736. Özişik, M.N. and Orlande, H.R.B., 2000, Inverse Heat Transfer: Fundamentals and Applications, Taylor and Francis, New York, NY. Phillips, D.L., 1962, “A Treatment for the Numerical Solution of Certain Integral Equations of the First Kind,” J. Assoc. Computing Machinery, Vol. 9, pp. 84-97. Siegel, R. and Howell, J.R., 2002, Thermal Radiation Heat Transfer, 4th ed., Taylor and Francis, New York, NY. Tikhonov, A.N., Solution of Incorrectly Formulated Problems and the Regularization Method, Soviet Math. Dokl., 4, 1035-1038, 1963,. [Engl. trans. Dokl. Akad. Nauk. SSSR, 151, 501-504]. Wing, G.W., 1991, A Primer on Integral Equations of the First Kind, SIAM, Philadelphia, PA.
Problems
10.1 Following the derivation in eqs. (10.7) – (10.12), show that the configuration factor definition of eq. (10.12) is valid for non-black diffuse surfaces as well as for black surfaces. Using the relation F3− (1+ 2) = F3−1 + F3− 2 , show whether the relation F(1+ 2) − 3 = F1− 3 + F2 − 3 is also valid. For the cylindrical enclosure of diameter d and length l = 2d shown in Fig. P10.1, find all configuration factors among the surfaces. Using the factor for an "Infinitely long enclosure formed of three plane areas" derived in Example 10.3, find the factor from one side of an infinitely long enclosure with cross-section of an isosceles triangle with apex angle α to its base (Fig. P10.2). An enclosure is formed of three mutually perpendicular isosceles triangles of short side S and hypotenuse H and an equilateral triangle of side H (Fig P10.3). Find the configuration factor F1-2 between two perpendicular isosceles triangles sharing a common short side.
10.2 10.3 10.4
10.5
868 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
A1 A2
A3
Fig. P10.1
α
H A1 S
S A2
H
S H Fig. P10.3 H= 20 cm
Fig. P10.2
20 cm
220 cm
S
Fig. P10.4
10.6 10.7
A plate A1 in Fig. 10.4 is to be moved along positions from S = 0 cm to S = 100 cm. Plot the configuration factor F1-2 vs. S for this configuration. A black electrically heated rod is in a black vacuum jacket (Fig. P10.5). The rod must dissipate 150 W without exceeding 1200 K. Calculate the maximum allowable jacket temperature (neglect end effects).
Chapter 10 Heat Transfer by Radiation 869
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
0.5 cm
vacuum
150 W 1200 K 2 cm
40 cm
Fig. P10.5
10.8 10.9
Derive eq. (10.47) using the net radiation method. A three-surface enclosure has the properties shown in Fig. P10.6. Assuming each surface can be treated as having uniform radiosity, find the unknown temperatures and heat fluxes. L1 L2 L3 L(m) 0.8 1.0 1.2 1.0 0.67 0.5
L2 L1 L3
Fig. P10.6
ε
T(K) 1000 300
q"(W/m2) 0
10.10
A four-surface enclosure has the properties and temperatures shown in Fig. P10.7. Find the heat flux for each surface, and demonstrate that the radiative heat balance is satisfied. Assume that radiosities can be taken as uniform across each surface.
A1, T1=850K ε1= 1.0 A4, T4=400K ε4 = 0.75 A2, T2=1300K ε2 = 0.85 A3, T3=1000K ε3 = 0.71
H=1 m
L=4m
Fig. P10.7
10.11
Derive a relation for the heat transfer q through an N-layer set of concentric circular cylindrical shields placed between an inner surface at
870 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
T1 and an outer layer at T2. The outer surface has diameter D1, and the inner surface has diameter D2. All surfaces have an emissivity of ε. 10.12 A shield system of N large closely-spaced flat-plate thin shields with all shield surfaces of emissivity εs is placed between an outer surface at T1 with emissivity ε1 and another surface at T2 with emissivity ε2. Derive an expression for the heat flux q" between surfaces 1 and 2. Solar radiation with a flux of q"solar = 1353 W/m2 is normally incident on a front surface of a shield configuration. The front surface has solar absorptivity αsolar = 0.10, and with IR emissivity of 0.86. The back of the surface has emissivity 0.015, and faces a 4-layer shield system; each surface of the shields has an emissivity of εs = 0.015. The outer face of the last layer of the shield faces space which can be taken as having a temperature of 4 K. What is the temperature of each layer (i.e., the front surface and each of the shields)? It is proposed to add a cryogenic cooling loop to the face of the last layer of the shield system described in Homework Problem 10.13. The cooler is needed to reduce the last shield face temperature to below 30 K. What cooling capacity q"c (mW/m2) is required to achieve this temperature? Solve the problem in Example 10.10 under the assumption of uniform radiosity on each surface, and compare with the solution given in Example 10.10. An uncovered styrofoam pan filled with water is placed outdoors on a cloudy night. There is almost no wind, so the heat transfer coefficient between the air and water is h = 15 W/(m2•K). The night sky acts as blackbody surroundings at Te = 240 K. Water is opaque for long wavelength radiation. The index of refraction of water is 1.33. What is the minimum air temperature required to prevent freezing?
1100 K ε=1 300 K gas
10.13
10.14
10.15
10.16
0.15 cm diameter thermocouple
Fig. P10.8
Chapter 10 Heat Transfer by Radiation 871
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
10.17
A copper-constantan thermocouple (ε=0.22) is in a transparent gas stream at 300 K adjacent to a large blackbody surface at 1100 K (Fig. P10.8). The heat transfer coefficient from the flowing gas to the Estimate the bare thermocouple thermocouple is 25 W/(m2•K). temperature. A wire between two electrodes is heated electrically with a total of qe watts. The wire resistivity re ohm-cm is constant, and the current is I. One end of the wire is at T1 and the other is at T2. The immediate surroundings are a vacuum, and the surroundings have a radiating temperature of To. The wire is gray with emissivity ε, diameter D, thermal conductivity k, and length L. Set up the differential equation for the steady-state temperature distribution along the wire, neglecting radial temperature variations within the wire. Integrate the equation, and find an expression (in the form of an integral) for the wire temperature as a function of x. The final results should contain only the quantities given, and should not contain dT/dx. Explain how you would evaluate the expression to find T(x). L T1 D Wire
Fig. P10.9
10.18
T2
ε
To Vacuum
10.19
In the inverse problem posed below, the two black surfaces exchange radiant energy. The surroundings are at T∞ = 0. The general problem is: Given T1(x1) = constant = T1 and q"1(x1) is specified, find the required distribution of temperature T2(x2) and net heat flux q"2(x2) on the upper surface. Show that at x1 on the lower surface, the governing equation for this geometry and boundary conditions is
W
x2 = 0
Eb 2 ( x 2 ) K ( x1 , x 2 )dx 2 = Eb1 − q1′′( x1 )
where Eb 2 ( x 2 ) = σ T24 ( x2 ) and that the energy equation can be nondimensionalized using Q ′′ = q′′ / σ T24 , X = x/L, and θ = Eb / σ T14 to the form 1 X =0 θ 2 ( X 2 ) K ( X1 , X 2 )dX 2 = 1 − Q1′′( X1 )
2
872 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
The kernel of the integral is given by K ( X1 , X 2 ) = 10.20
1 h2 2 h 2 + ( X − X ) 2 3/ 2 2 1
where h =H/W. For the geometry in Problem 10.19, set the heat flux distribution on the lower surface to Q1′′( X1 ) = AX12 − BX1 − C , discretize the governing integral equation, and provide the resulting set of algebraic equations in matrix form for determining the unknown values of θ2(X2). Note that there are parameters h, (A, B, C), and the number of discrete elements m and n of size ΔX1 and ΔX2 on each surface.
W
ε2=1, T (x 2) =?
H x1 T 1, q" 1(x 1), ε1 =1 T∞ = 0
x2 T ∞= 0
Fig. P10.10
10.21 For the geometry in Problem 10.19 and the particular case Q1′′( X1 ) = AX12 − BX1 − C , h = 0.6, A = B = 16 and C = 6, and m = n = 30, solve for θ2,i for the elements on the top surface using a standard matrix inversion routine (Mathcad, Matlab, Maple, etc.) 10.22 10.23 10.24 10.25 For the problem stated in Problem 10.21, use SVD on the A matrix of the problem, and show a plot of the singular values. Choosing an appropriate number of singular values, use TSVD to find a reasonable distribution for θ(x2) for the parameters of Homework 10-21. Prove that for isotropic scattering, the phase function Φ = 1. Two large parallel gray plates are 5.65 cm apart. Their temperatures and emissivities are T1 = 775 K, ε1 = 0.70; T2 = 512 K, ε2 = 0.31. Compute the heat transferred by radiation across the space between the plates when the space is a vacuum, and when the space is filled with a gray medium that has absorption and isotropic scattering with extinction coefficient κ = 0.63 cm-1. (Use the diffusion method as an approximation and neglect heat conduction.) For water vapor at 1000K with a mean beam length of 2 m and a partial pressure of 0.2 atm when mixed with air at a total pressure of 1 atm, compare the water vapor total emittance ε(T, pLe) using the predictions of Cess and Lian, eq. (10.168) and Leckner, eq. (10.169).
10.26
Chapter 10 Heat Transfer by Radiation 873
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press
10.27
Use the Leckner correlation to find the emittance of the Earth's atmosphere for a temperature of 300K, a CO2 concentration of 350 ppm by volume, and a mean beam length through the atmosphere of 60 km. Pure carbon dioxide at 1 atm and 2778 K is contained between parallel plates 0.3 m apart. What is the radiative flux received at the plates as a result of radiation by the gas? (Use the Leckner relations for CO2 total emittance.) A furnace at atmospheric pressure with interior in the shape of a a cylinder with height equal to two times its diameter is filled with a 50:50 mixture by volume of CO2 and N2. The gas temperature is uniform at 1825 K and the walls are cooled. The interior surfaces are black. At what rate is energy being supplied to the gas (and removed from the walls) to maintain these conditions? Use the Leckner correlations for gas properties. A rectangular furnace of dimensions 0.3×0.5×2.2 m has soot covered interior walls that can be considered black and cold. The furnace is filled with well-mixed combustion products at a temperature of 2220 K composed of 40% by volume CO2, 30% by volume water vapor, and the remainder N2. The total pressure is 1 atm. Compute the radiative flux and the total energy rate to the walls using the Leckner correlations for the CO2 and H2O mixture.
10.28
10.29
10.30
874 Advanced Heat and Mass Transfer
Amir Faghri, Yuwen Zhang, and John Howell
Copyright © 2010 Global Digital Press