Linearization of Nongray Radiation Exchange: The Internal Fractional Function Reconsidered

Lienhard V, John H.

doi:10.1115/1.4042158

The radiation fractional function is the fraction of black body radiation below a given value of λT. Edwards and others have distinguished between the traditional, or “external,” radiation fractional function and an “internal” radiation fractional function. The latter is used for linearization of net radiation from a nongray surface when the temperature of an effectively black environment is not far from the surface's temperature, without calculating a separate total absorptivity. This paper examines the analytical approximation involved in the internal fractional function, with results given in terms of the incomplete zeta function. A rigorous upper bound on the difference between the external and internal emissivity is obtained. Calculations using the internal emissivity are compared to exact calculations for several models and materials. A new approach to calculating the internal emissivity is developed, yielding vastly improved accuracy over a wide range of temperature differences. The internal fractional function should be used for evaluating radiation thermal resistances, in particular.

Introduction

The use of the radiation fractional function is well known in computing total emissivity from spectral emissivity. Less well known is the use of the internal radiation fractional function for approximating net heat flux from one surface to another at a similar temperature. Edwards was a strong proponent of the latter approach, beginning in the 1960s [1], following work published a few years earlier by Czerny and Walther [2]. Edwards' work on radiative property measurements had well acquainted him with the failure of gray-body exchange approximations at even small temperature differences [3], leading him to seek a better approximation that retained the computational convenience of using only a single radiative property.

Edwards distinguished the traditional fractional function by calling it the external fractional function. The terms external and internal derive from Edwards' work in the aerospace industry, differentiating between radiation exchange on the exterior of a spacecraft, where temperature differences are large, and on the interior of a spacecraft, where surface temperatures are not far apart [4]. The internal radiation fractional function appears in several textbooks by Edwards and coworkers [5–7] spanning a period of decades. Most other textbooks have not taken it up. In part, the lack of adoption may reflect the limited range of temperature difference over which the internal emissivity provides an accurate estimate of the net heat flux.

The present paper reconsiders the approximation leading to the internal radiation fractional function, with the objective of understanding its accuracy and how that may be improved. The results are given in terms of the incomplete zeta function, and a rigorous upper bound on the difference between the external and internal emissivity is obtained. Calculations based upon the internal radiation fractional function are compared to exact calculations, and a procedure is established that sharply increases domain of accuracy of the internal emissivity. The internal emissivity should be used in computing radiation thermal resistance in multimode heat transfer and for other linearizations of nongray radiation exchange.

Background

The spectral hemispherical absorptivity, α(λ, T₁), of a surface at T₁, is related to the total hemispherical² emissivity, ε(T₁), by integration over the Planck distribution of blackbody monochromatic emissive power, e_λ_,_b(T₁)

ε (T_{1}) = \frac{1}{σ T_{1}^{4}} \int_{0}^{\infty} α (λ, T_{1}) e_{λ, b} (T_{1}) d λ

(1)

e_{λ, b} (T_{1}) = \frac{π c_{1}}{λ^{5} [exp (c_{2} / λ T_{1}) - 1]}

(2)

where c₁ and c₂ are the first and second radiation constants. If this surface is irradiated by an effectively black environment at T₂ (for example, if the surroundings are much larger than the surface), the total hemispherical absorptivity is

α (T_{1}, T_{2}) = \frac{1}{σ T_{2}^{4}} \int_{0}^{\infty} α (λ, T_{1}) e_{λ, b} (T_{2}) d λ

(3)

The net radiation leaving this surface is

q_{net} = σ ε (T_{1}) T_{1}^{4} - σ α (T_{1}, T_{2}) T_{2}^{4}

(4)

Nongray Errors

For a gray surface, α(λ, T) is independent of wavelength and so ε(T₁) = α(T₁, T₂). However, some older literature went farther, applying this idea to nongray surfaces for which T₁ is near T₂. Edwards noted that the fact that $\lim_{T_{2} \to T_{1}} α (T_{1}, T_{2}) = ε (T_{1})$ is meaningless in calculating heat exchange; and further, working with Nelson, he documented associated “nongray errors” in radiative property measurements that exceeded 60% [3].

The reason for such errors is not that the limit is incorrect, but that the heat flux is determined by the slope of q_net as a function of temperature difference, which changes if the surface is not truly gray. To see why in a simple way, the last term in Eq. may be linearized

α (T_{1}, T_{2}) T_{2}^{4} \approx α (T_{1}, T_{1}) T_{1}^{4} + {\frac{d}{d T_{2}} (α (T_{1}, T_{2}) T_{2}^{4}) |}_{T_{1}} (T_{2} - T_{1})

(5)

= ε (T_{1}) T_{1}^{4} + 4 T_{1}^{3} [ε (T_{1}) + {\frac{T_{1}}{4} \frac{d α}{d T_{2}} |}_{T_{1}}] (T_{2} - T_{1})

(6)

Thus

q_{net} \approx 4 σ T_{1}^{3} [ε (T_{1}) + {\frac{T_{1}}{4} \frac{d α}{d T_{2}} |}_{T_{1}}] (T_{1} - T_{2})

(7)

For a gray (or black) surface, dα/dT₂ = 0, and we recover a common linearization: $q_{net} \approx 4 σ ε T_{1}^{3} Δ T$ [11]. Thus, the purpose of introducing the internal fractional function is to account for the difference in slope when dα/dT₂ ≠ 0.

External and Internal Radiation Fractions

The fraction of blackbody radiation between wavelengths of 0 and λ is

f (λ T) = \frac{1}{σ T^{4}} \int_{0}^{λ} e_{λ, b} d λ

(8)

= 1 - \frac{15}{π^{4}} \int_{0}^{c_{2} / λ T} \frac{t^{3}}{e^{t} - 1} d t

(9)

= 1 - \frac{90}{π^{4}} ζ (c_{2} / λ T, 4)

(10)

ζ(X,s) is the incomplete zeta function (see derivation in Appendix A). Edwards [1,4,5] refers to f as the external radiation fraction. The distinct internal radiation fractional function, f_i(λT), arises from a linearization that allows ε(T_s) and α(T_s,T_e) to be combined, as described next.

Edwards [1] defines the internal hemispherical emissivity such that

ε^{i} (T_{1}) \equiv \lim_{T_{2} \to T_{1}} \frac{σ ε (T_{1}) T_{1}^{4} - σ α (T_{1}, T_{2}) T_{2}^{4}}{σ T_{1}^{4} - σ T_{2}^{4}}

(11)

Because the numerator and denominator both go to zero in the limit, L'Hôpital's rule (or Taylor expansion) is required

ε^{i} (T_{1}) = \lim_{T_{2} \to T_{1}} \frac{\int_{0}^{\infty} α (λ, T_{1}) \frac{\partial}{\partial T_{2}} e_{λ, b} (T_{2}) d λ}{4 σ T_{2}^{3}}

(12)

= \int_{0}^{\infty} α (λ, T_{1}) [\frac{1}{4 σ T_{1}^{3}} {\frac{\partial}{\partial T_{2}} e_{λ, b} (T_{2}) |}_{T_{1}}] d λ

(13)

Thus, when T₂ is not too much different from T₁

q_{net} \approx ε^{i} (T_{1}) 4 σ T_{1}^{3} (T_{1} - T_{2})

(14)

with

ε^{i} (T) = \frac{1}{4 σ T^{3}} \int_{0}^{\infty} α (λ, T) \frac{\partial e_{λ, b}}{\partial T} d λ = \int_{0}^{1} α (λ, T) d f_{i} (λ T)

(15)

where the internal radiation fractional function is

f_{i} (λ T) \equiv \frac{1}{4 σ T^{3}} \int_{0}^{λ} \frac{\partial e_{λ, b}}{\partial T} d λ

(16)

L'Hôpital's rule simply compares the slopes of the numerator and denominator in the limit, so Eqs. (7) and (14) are identical statements, as a change of variables and integration by parts of Eq. (13) shows.

Equation stands in contrast to the external total hemispherical emissivity, which from Eqs. and maybe written

ε (T) = \frac{1}{σ T^{4}} \int_{0}^{\infty} α (λ, T) e_{λ, b} (T) d λ = \int_{0}^{1} α (λ, T) d f (λ T)

(17)

From these relationships, one can show that

f_{i} (λ T) - f (λ T) = F (X) = \frac{15}{4 π^{4}} \frac{X^{4}}{e^{X} - 1}

(18)

where X ≡ c₂/λT. The functions f, f_i, and F are plotted as a function of λT in Fig. 1.

The numerical implementation is described in Appendix B.

Difference Between External and Internal Emissivities

We now derive an upper bound on the difference of these two emissivities

ε - ε^{i} = \int_{0}^{1} α (λ, T) d f (λ T) - \int_{0}^{1} α (λ, T) d f_{i} (λ T)

(19)

= \int_{0}^{\infty} α (λ, T) \frac{d F}{d X} d X

(20)

= \int_{0}^{X_{z}} α (λ, T) \frac{d F}{d X} d X + \int_{X_{z}}^{\infty} α (λ, T) \frac{d F}{d X} d X

(21)

where X = X_z is the finite zero point of dF/dX, X_z = 3.92069 (see Appendix C). The function dF/dX is positive for X < X_z and negative for X > X_z, and

0 < α (λ, T) ⩽ 1

⁠. Thus, we can form two bounds on the difference

ε - ε^{i} ⩽ \int_{0}^{X_{z}} \frac{d F}{d X} d X = F (X_{z}) if ε - ε^{i} > 0, and

(22)

ε^{i} - ε ⩽ \int_{\infty}^{X_{z}} \frac{d F}{d X} d X = F (X_{z}) if ε^{i} - ε > 0

(23)

A simple calculation shows that

F (X_{z}) = \frac{15}{π^{4}} X_{z}^{3} e^{- X_{z}} = 0.18400

(24)

In other words

\begin{matrix} | ε - ε^{i} | ⩽ 0.18400 \end{matrix}

(25)

The value of X_z corresponds to

{(λ T)}_{z} = \frac{14387.77}{3.92069} = 3699.70 μ m K

(26)

At 300 K, the wavelength λ_z = 12.2323 μm.

The following two extreme examples are designed to illustrate these bounds. First, suppose that we have a 300 K surface that is black for

λ ⩽ 12.2323 μ m

⁠, but perfectly reflective on other wavelengths. Numerical integration of Eqs. and yields

ε = 0.4177, ε^{i} = 0.6017, and ε^{i} - ε = 0.1840

(27)

If instead we have a 300 K surface that is black for

12.2323 μ m ⩽ λ

⁠, but perfectly reflective on other wavelengths

ε = 0.5823, ε^{i} = 0.3983, and ε - ε^{i} = 0.1840

(28)

In the first case, εⁱ is 44% greater than ε; in the second case, εⁱ is 32% less than ε.

For a gray surface, Eq. (20) gives ε − εⁱ = 0, consistent with Eqs. (4) and (7). The variation of the total absorptivity for a 300 K surface irradiated by an effectively black source at T₂ is shown for both example surfaces in Fig. 2. The total absorptivity varies significantly. The net heat exchange is discussed in detail below.

The emissivities obtained from these examples are more general than they may appear: the same values will be obtained for any surface for which α(λ) switches between zero and one at X_z. For example, from the preceding development, for the surface black on short wavelengths

ε = \int_{0}^{λ_{z}} α (λ) d f (λ T_{z}) = \int_{\infty}^{X_{z}} \frac{d f}{d X} d X = f (X_{z}) - 0

(29)

= 1 - \frac{90}{π^{4}} ζ (X_{z}, 4) = 1 - \frac{90}{π^{4}} ζ (3.92069, 4)

(30)

= 0.41771

(31)

Equation (18) gives εⁱ, and the results for the long wavelength example may be obtained by subtracting these values from unity.

Linearization Errors

The linearization itself produces error in calculation of the heat flux, apart from nongray error. To quantify this error, we may compare the linearized and exact values of q_net for a black surface. Their ratio is

\frac{4 T_{1}^{3} Δ T}{T_{1}^{4} - T_{2}^{4}} = \frac{4 T_{1}^{3}}{(T_{1}^{2} + T_{2}^{2}) (T_{1} + T_{2})} = \frac{4}{[1 + {(T_{2} / T_{1})}^{2}] (1 + T_{2} / T_{1})}

(32)

Equation (32) is plotted in Fig. 3. For T₂ = 320 K and T₁ = 300 K, the linearized flux is 90.5% of the exact flux. The ratio is asymptotic to $1 - 3 Δ T / 2 T_{1} + 5 {(Δ T)}^{2} / 4 T_{1}^{2}$ as ΔT → 0.

The linearization around T₁ is considerably less accurate than the traditional approach (i.e., for radiation heat transfer coefficients) of linearizing around T_m = (T₁ + T₂)/2. For that case, the ratio is

\frac{4 T_{m}^{3} Δ T}{T_{1}^{4} - T_{2}^{4}} = \frac{4 T_{m}^{3}}{(2 T_{m}^{2} + Δ T^{2} / 2) (2 T_{m})} = \frac{1}{1 + Δ T^{2} / 4 T_{m}^{2}}

(33)

and has errors on the order of (ΔT/T_m)²/4 [11]. As seen in Fig. 3, this linearization is far more accurate. The reason for the greater accuracy is that linearization about T_m is equivalent to a second-order Runge–Kutta method, with third-order truncation error in heat flux, while linearization about T₁ is simply a forward Euler method (see Appendix D).

From these considerations, we see that εⁱ will also have a limited range of accuracy if evaluated at T₁. Instead, εⁱ should be computed at T_m, following the analysis in Appendix D:

ε^{i} (T_{m}) = \frac{1}{4 σ T_{m}^{3}} \int_{0}^{\infty} α (λ, T_{1}) {\frac{\partial e_{λ, b}}{\partial T} |}_{T_{m}} d λ = \int_{0}^{1} α (λ, T_{1}) d f_{i} (λ T_{m})

(34)

with the heat flux evaluated, to an accuracy of

O (Δ T^{3})

⁠, as

\begin{matrix} q_{net} \approx 4 ε^{i} (T_{m}) \cdot σ T_{m}^{3} Δ T \end{matrix}

(35)

Comparison of Models

We compare the heat flux predicted by the exact result (Eq. (4)), the gray approximation, εⁱ(T₁) (Eq. (14)), and εⁱ(T_m) (Eq. (35)) in Figs. 4 and 5 for the same two example surfaces at T₁ = 300 K. The charts show

(1)
The gray surface approximation, ε = α, has the wrong slope as T₂ → T₁, consistent with Eq. (7). This approximation has poor accuracy for even the smallest temperature differences. The ratio of gray slope to exact slope is ε/εⁱ. Using the values from Eq. (27) or Eq. (28), that ratio is equal to 0.69 and 1.46 for Figs. 4 and 5, respectively. As expected, the gray approximation shows the greatest error for T₂ > T₁, where irradiation dominates net exchange.
(2)
The approximation using εⁱ(T₁) has the correct slope as T₂ → T₁, but is accurate for only a limited range of temperatures. For the surface black on short wavelengths, the magnitude of the linearized flux from Eq. (14) at T₂ = 320 K is 13% less than the exact value, Eq. (4). This case corresponds to the seemly small temperature difference ΔT/T₁ = −0.067.
(3)
The approximation using εⁱ(T_m) remains accurate over a very broad range of temperature variation, to the degree that the curves for this equation and the exact equation are difficult to distinguish. The charts in Figs. 4 and 5 include values as high as ΔT/T₁ = ±0.33.

Although Figs. 4 and 5 use 300 K as a reference point, Eqs. (29)–(31) show that, for a surface at T₁ = T_z that switches from black to reflective behavior at λ_z, both ε(T_z) and εⁱ(T_z) are fixed. The values of $ε^{i} (T_{m}) = ε^{i} (T_{z} - Δ T / 2)$ and α(T_z, T_z − ΔT) depend additionally only upon ΔT/T_z. Hence, the trends seen will be the same at any other temperature level. These comparisons are put into nondimensional form in Appendix F.

Edwards [1,4,5] and, later, Mills [7] were both very specific in recommending to evaluate εⁱ at T₁. Edwards et al. [6] did not specify the temperature, but without comment included an example worked using T_m. Likewise, Czerny and Walther [2] provided an example worked using T_m, again without comment. The present results show clearly that only T_m should be used when evaluating εⁱ.

Application to Representative Materials

The previous example surfaces were contrived to maximize the differences between ε and εⁱ. In this section, we examine the differences for representative nongray materials.

Alumina.

Experimental data for polycrystalline alumina (99.5% Al₂O₃, 6 mm thick, 1 μm roughness) have been reported by Teodorescu and Jones [12] for T = 823 K and $2 ⩽ λ ⩽ 25 μ m$ ⁠. The experimental uncertainty is 3.5%. These data are available in tabular form on 1 μm intervals and as a chart (Fig. 6). These results are generally consistent with other reports on polycrystalline alumina.

In this case, we first integrate the data to obtain the spectral normal emissivity, ε_n. One additional data point is extrapolated, ε(1 μm, 823 K) ≈ 0.128, because f_i at 2 μm still has the finite value 0.059; the extrapolation extends the range to f_i = 0.00012. At the other end of this spectrum, f = 0.987 at 25 μm, but no extrapolation was done. Thus, the value for total normal emissivity ε_n is thus likely to be low by ∼1%. Integration of the discrete data by a simple trapezoidal rule yields the following values, both at 823 K:

ε_{n} = 0.506, ε_{n}^{i} = 0.404, and ε_{n}^{i} - ε_{n} = 0.102

(36)

Our interest, ultimately, lies with estimating heat exchange, so we require the spectral hemispherical emissivity. Teodorescu and Jones [12] also measured the spectral directional emissivity. We have integrated those data, as described in Appendix E, to obtain the spectral hemispherical values from 1 to 25 μm.

We consider exchange with an effectively black environment at temperatures down to 625 K. At 625 K, f = 0.972 at 25 μm so that the absorptivity calculation may be low by 2% or so. The results are shown in Fig. 7. Once again, Eq. (35) using εⁱ(T_m) is in excellent agreement with the exact result, Eq. (4).

Drude/Hagens-Ruben Metal.

The Drude/Hagen–Rubens model for emissivity applies approximately to some metals in the infrared regime [5,10,13]. Baehr and Stephan [13] provide the following equation for the spectral, hemispherical emissivity under this model:

\begin{matrix} ε (λ, T) = 48.70 \sqrt{\frac{r_{e}}{λ}} {1 + [31.62 + 6.849 \ln (\frac{r_{e}}{λ})] \sqrt{\frac{r_{e}}{λ}} - 166.78 \frac{r_{e}}{λ} + \dots} \end{matrix}

(37)

for

r_{e} / λ < 5 \times 10^{- 4} Ω \cdot cm / μ m

⁠, where r_e(T) is the electrical resistivity in Ω cm.

We consider a surface at T₁ = 373 K with r_e = 13.1 × 10⁻⁶ Ω ⋅ cm, a value representative of platinum. At this temperature, $0.01 ⩽ f ⩽ 0.99$ for $3.88 ⩽ λ ⩽ 61.3 μ m$ and f_i = 0.01 at λ = 3.36 μm. Over this range of wavelengths, ε(λ) varies from 0.086 down to 0.022 according to Eq. (37). We may consider heat exchange with a surface at T₂ down to 275 K; at 275 K, f > 0.99 for λ > 83.1 μm. The integrations are, therefore, done from 3 to 100 μm. The results are shown in Fig. 8.

Once again, Eq. using εⁱ(T_m) closely tracks the exact solution over a wide range (down to at least 300 K or ΔT/T₁ = 20%). In this instance, nongray effects are present but less pronounced; the gray approximation, while not the best alternative for even small ΔT, shows less divergence than in previous examples. In addition, at 373 K

ε = 0.0497, ε^{i} = 0.0553, and ε^{i} - ε = 0.0056

(38)

The internal and external emissivities are closer, as would be expected when nongray effects are smaller.

Spectrally Selective Materials.

We may also consider an idealized spectrally selective surface that makes a step transition at λ_c from a short-wavelength absorptivity, α_sw, to a long-wavelength absorptivity, α_lw

α (λ) = {\begin{matrix} α_{sw} & for λ ⩽ λ_{c} \\ α_{lw} & for λ > λ_{c} \end{matrix}

(39)

From Eqs. and

ε (T) = α_{sw} f (λ_{c} T) + α_{lw} [1 - f (λ_{c} T)]

(40)

= α_{sw} + \frac{90}{π^{4}} Δ α ζ (X_{c}, 4)

(41)

where X_c = c₂/λ_cT and Δα = α_lw − α_sw. Similarly, from Eqs. and

ε^{i} (T) = α_{sw} f_{i} (λ_{c} T) + α_{lw} [1 - f_{i} (λ_{c} T)]

(42)

= α_{sw} + Δ α [\frac{90}{π^{4}} ζ (X_{c}, 4) - F (X_{c})]

(43)

Further

ε (T) - ε^{i} (T) = Δ α F (X_{c})

(44)

and also

ε^{i} (T_{m}) = α_{sw} + Δ α [\frac{90}{π^{4}} ζ (X_{c, m}, 4) - F (X_{c, m})]

(45)

where X_c_,_m = c₂/λ_cT_m. Finally, from Eq.

α (T_{1}, T_{2}) = α_{sw} + \frac{90}{π^{4}} Δ α ζ (X_{c, 2}, 4)

(46)

with X_c_,2 = c₂/λ_cT₂. The previous results show that the impact of selectivity will be greatest when X_c and X_z are close.

As an example, we may consider the following crude approximation to the soft-anodized aluminum described in Ref. [5]: α_sw = 0.10, α_lw = 0.85, and λ_c = 7 μm (Fig. 9). If this selective solar reflector is at T₁ = 360 K and exchanges radiation with sky at T₂ = 290 K, we obtain the results in Table 1. At these temperatures, most of the radiant energy is on wavelengths above 7 μm, and nongrayness is not pronounced. Nonetheless, as in previous examples, the agreement between Eq. (4), 371.8 W/m², and Eq. (35), 371.0 W/m², is excellent, whereas Eq. (14), 462.1 W/m², performs quite poorly.

Equations (39)–(46) are used further in Appendix F.

Accuracy Versus Effort

One rationale for introducing εⁱ is simply to avoid computing the total hemispherical emissivity and absorptivity separately when the spectral absorptivity has a significant wavelength dependence. The same labor-saving rationale was used in some older literature to justify a gray body approximation, ε(T₁) ≈ α(T₁, T₂), when the temperatures are close. When the internal fractional function was first developed, personal computers (and pocket calculators, spreadsheets, symbolic integration, etc.) did not exist. Routine analyses required either tedious hand calculations or cumbersome trips to the mainframe computer center (to type up Hollerith cards and then wait for the SysOp to run the job!). So, one attraction of the method was clearly to reduce labor.

While numerical integration is very easy today, εⁱ(T_m) can still save some time in doing basic calculations. The greater challenge may be to obtain complete spectral, hemispherical data spanning the necessary range of wavelengths. If, instead, electromagnetic theory is used, working from the complex refractive index and calculating spectral directional emissivity, much coding has to be undertaken, and the single integration saved by εⁱ will not be at all significant.

However, a second very important consideration is that a linearized expression for radiation heat transfer allows the use of a radiation thermal resistance, perhaps in parallel to a natural convection resistance, in network analyses of heat transfer. For an area A, such a resistance has the form

R_{rad} = \frac{1}{4 ε σ T^{3} A}

(47)

R_rad is routinely evaluated using a gray body approximation. For small temperature differences, using the internal emissivity, εⁱ(T₁), in Eq. will provide a more accurate value than a gray body approximation. The present work has shown εⁱ(T_m) to be far more robust, with excellent accuracy for

| Δ T / T_{1} |

of 20–30%. Hence, the radiation thermal resistance is best calculated as

\begin{matrix} R_{rad} = \frac{1}{4 ε^{i} (T_{m}) σ T_{m}^{3} A} \end{matrix}

(48)

Despite modern computing power, the need for an accurate radiation resistance remains as important today as it was 60 years ago. Thus, we conclude that the internal fractional function has an ongoing value when surfaces do not politely behave as “gray bodies.”

Summary

The internal radiation fractional function has been in the literature since the late 1950s, although not widely adopted in textbooks covering thermal radiation. The most visible recommendations have been from Edwards and coworkers, generally advising the use of εⁱ(T₁) to linearize nongray exchange for enclosures with modest temperature differences. In this study, the following new findings are reported:

(1)
The maximum absolute difference between ε(T₁) and εⁱ(T₁) is 0.18400, Eq. (25).
(2)
The internal emissivity should be evaluated at the mean temperature, T_m, not T₁ as has often been suggested. The differences in accuracy are shown to be significant because the evaluation at T_m using Eq. (35) is effectively a second-order accurate numerical method, with truncation error in the heat flux of $O (Δ T^{3})$ ⁠. In contrast, evaluation at T₁ is a first-order approximation.
(3)
Equation (35) with εⁱ(T_m) has excellent accuracy for $| Δ T / T_{1} |$ up to 20–30%. Beyond this range, the linearization error of Fig. 3 will cause increasing error for any surface.
(4)
Theory and examples for several nongray materials show that the gray-body approximation gives the wrong slope for heat flux as T₂ → T₁.
(5)
Calculations involving both the internal and external fractional functions can be conveniently implemented using the incomplete zeta function.
(6)
εⁱ(T_m) should be used in calculating radiation thermal resistances for nongray surfaces, with Eq. (48).

Acknowledgment

Don Edwards was on the faculty of UCLA's CNTE Department when I studied there, although he was not my instructor. Nevertheless, I felt his positive influence through the curriculum and the other faculty in the heat transfer program.

Nomenclature

c₁ =
first radiation constant, 3.741772 × 10⁻¹⁶ (W m²)

c₂ =
second radiation constant, 14387.77 (μm K)

e_λ_,_b =
blackbody monochromatic emissive power (W m⁻³)

f =
(external) radiation fractional function

F =
f_i − f, Eq. (18)

f_i =
internal radiation fractional function

q_net =
net radiant heat flux from surface 1 to 2 (W m⁻²)

r_e =
electrical resistivity (Ω cm)

R_rad =
radiation thermal resistance, Eq. (47) (K W⁻¹)

T =
temperature (K)

T_m =
mean temperature, (T₂ + T₁)/2 (K)

T₁ =
temperature of surface (K)

T₂ =
temperature of environment (K)

X =
c₂/λT

X_z =
finite value of X for which dF/dX = 0, 3.92069

Y =
function in Appendix C, Eq. (D1)

α_lw, α_sw =
see Eq. (39)

α(T₁, T₂) =
total hemispherical absorptivity

α(λ, T) =
spectral hemispherical absorptivity

δT =
T₂ − T₁ (K)

ΔT =
T₁ − T₂ (K)

Δα =
α_lw − α_sw

ε(T) =
total hemispherical emissivity (external)

ε_n(T) =
total normal emissivity (external)

ε_n(λ, T) =
spectral normal emissivity

ε(λ, T) =
spectral hemispherical emissivity

$ε' (θ, λ, T)$ =
spectral directional emissivity

εⁱ(T) =
internal total hemispherical emissivity

$ε_{n}^{i} (T)$ =
internal total normal emissivity

ζ(s) =
Riemann zeta function

ζ(X, s) =
incomplete zeta function

θ =
polar angle (rad)

λ =
wavelength (m)

σ =
Stefan–Boltzmann constant, 5.670367 × 10⁻⁸ (W m⁻² K⁻⁴)

2

The spectral hemispherical absorptivity and emissivity are equal if the irradiation is diffuse: α(λ, T₁) = ε(λ, T₁) [5,10]. The main points of this work also apply to spectral directional values and are so used in certain parts of this study.

3

Standard theory shows that for multistep integration, the global error of Euler's method is $O (Δ T)$ [30, Sec. 9.3].

Appendix A: Incomplete Zeta Function

The radiation fractional function may be written in terms of the incomplete zeta function for convenience

f (λ T) = \frac{1}{σ T^{4}} \int_{0}^{λ} \frac{2 π h c_{o}^{2}}{λ^{5} [exp (h c_{o} / k_{B} T λ) - 1]} d λ

(A1)

= \frac{1}{σ T^{4}} \frac{2 π k_{B}^{4} T^{4}}{h^{3} c_{o}^{2}} \int_{c_{2} / λ T}^{\infty} \frac{t^{3}}{e^{t} - 1} d t

(A2)

When λT → ∞, f = 1 and the last equation yields the well-known result

σ T^{4} = \frac{2 π k_{B}^{4} T^{4}}{h^{3} c_{o}^{2}} \underset{\equiv ζ (4) Γ (4)}{\underset{︸}{\int_{0}^{\infty} \frac{t^{3}}{e^{t} - 1} d t}}

(A3)

where the Gamma function Γ(4) = 3! and the Riemann zeta function, ζ(4), has the indicated integral representation [14, Sec. 13.12]. A classical result due to Euler [15] gives ζ(4) = π⁴/90 (see also Ref. [16, Sec. 167]), from which we recover the usual definition of the Stefan–Boltzmann constant, σ. Returning to Eq. with this information, we have

f (λ T) = \frac{15}{π^{4}} \int_{0}^{\infty} \frac{t^{3}}{e^{t} - 1} d t - \frac{15}{π^{4}} \int_{0}^{c_{2} / λ T} \frac{t^{3}}{e^{t} - 1} d t

(A4)

= 1 - \frac{15}{π^{4}} \int_{0}^{X} \frac{t^{3}}{e^{t} - 1} d t

(A5)

= 1 - \frac{15}{π^{4}} Γ (4) ζ (X, 4)

(A6)

where X = c₂/λT and we have identified the incomplete zeta function, ζ(X, s) [17, Sec. 8.22]. Hence

f (λ T) = 1 - \frac{90}{π^{4}} ζ (X, 4)

(A7)

Appendix B: Numerical Implementation

The integrals were computed using the GNU Scientific Library [18] with FFI bindings [19] to Lua code [20] under LuaLaTeX [21] using TEXShop over TEX Live [22]. The C code was compiled using XCode under Mac OS X. Integration was performed using GSL's QAG adaptive integration with 61 point Gauss–Kronrod rules. Convergence was checked, and the numerical values given in the text are believed to be accurate to the number of digits shown.

The incomplete zeta function was rendered in terms of the third-order Debye function, D₃(X) = 3 ζ(X, 4) Γ(4)/X³, and computed using GSL's built-in routine. These routines were supplied to PGFPLOTS [23] to generate the charts.

The integrations of the discrete data for alumina were handled outside GSL, as described in the main text and Appendix E.

Readers interested in implementing these numerical methods may refer to the following tutorial publications and posts: [24–27]. All of these packages (other than Mac OS X) are available online at no cost and are actively maintained.

Appendix C: The Constant Xz

X_z is the finite zero point of dF/dX, specifically the nonzero root of

4 (1 - e^{- X_{z}}) = X_{z}

(C1)

The equation may be solved by iteration to find

X_{z} = 3.92069 \dots

⁠. Alternatively, by defining y = 4 − X_z, Eq. becomes

y e^{y} = 4 e^{- 4}

(C2)

This equation may be solved using the Lambert W function

X_{z} = 4 - W (4 e^{- 4}) = 3.92069 \dots

(C3)

We may show by contradiction that X_z is irrational. Assuming X_z is rational, X_z = a/b for nonzero integers a and b. Then, the left-hand side of Eq. (C1) is 4 − 4e⁻^a/b; however, e^a/b is itself irrational (e.g., see Ref. [28, Sec. 4.7]). Hence, 4 − 4e⁻^a/b is an irrational number and cannot equal the assumed right-hand side, a/b.

Diophantine approximations to X_z may be constructed by calculating a continued fraction representation of X_z through the usual process

X_{z} = 3.92069 \dots = 3 + \frac{1}{1 + \frac{1}{11 + \frac{1}{1 + \frac{1}{1 + \frac{1}{⋱}}}}}

(C4)

Successive convergents (finite truncations) of the fraction give an infinite number of rational approximations to X_z

X_{z} \approx {4, \frac{47}{12}, \dots, \frac{149}{38}, \frac{247}{63}, \dots, \frac{1137}{290}, \dots}

(C5)

The second of these is within 0.1% of the exact value; the last agrees to six digits.

Appendix D: Linearizations as Finite Difference Approximations

We may consider q_net to be a function of T that we wish to determine approximately at T = T₂ by evaluating at other temperatures, i.e., T = T₁ and T = T_m. With

Y (T) = σ ε (T_{1}) T_{1}^{4} - σ α (T_{1}, T) T^{4}

(D1)

we seek to estimate

q_{net} = Y (T_{2}) = σ ε (T_{1}) T_{1}^{4} - σ α (T_{1}, T_{2}) T_{2}^{4}

(D2)

A first-order Euler method starting at T = T₁ would approximate Y (T₂) based only on conditions at T₁ as follows:

Y (T_{2}) \approx Y (T_{1}) + Y' (T_{1}) \cdot δ T = Y' (T_{1}) \cdot δ T

(D3)

for δT = T₂ − T₁ (= −ΔT), which gives Eq. or Eq. . The local truncation error of this single-step approximation is

O (Y^{″} Δ T^{2})

[29].³

A second-order Runge–Kutta method works from T_m with expansions toward both T₁ and T₂, subtracting the former from the latter

Y (T_{2}) = Y (T_{m}) + Y' (T_{m}) \frac{δ T}{2} + Y^{″} (T_{m}) \frac{δ T^{2}}{8} + O (δ T^{3})

(D4)

Y (T_{1}) = Y (T_{m}) - Y' (T_{m}) \frac{δ T}{2} + Y^{″} (T_{m}) \frac{δ T^{2}}{8} - O (δ T^{3})

(D5)

Y (T_{2}) = Y (T_{1}) + Y' (T_{m}) \cdot δ T + O (δ T^{3})

(D6)

Y (T_{2}) \approx Y' (T_{m}) \cdot δ T

(D7)

This single-step Runge–Kutta approximation has local truncation error

O (Y^{″'} Δ T^{3})

[29]. In these calculations, α(λ, T₁) is evaluated at the surface temperature, T₁. With Eqs. and

Y' (T_{m}) = - {\frac{d}{d T} (σ T^{4} α (T_{1}, T)) |}_{T_{m}}

(D8)

= - \int_{0}^{\infty} α (λ, T_{1}) {\frac{\partial e_{λ, b}}{\partial T} |}_{T_{m}} d λ

(D9)

= - 4 σ T_{m}^{3} \cdot ε^{i} (T_{m})

(D10)

which results in Eq. . Given lower truncation error of the Runge–Kutta approach, we should expect significantly better accuracy using Eq. than when using Eq. or Eq. ; and indeed the numerical results confirm this conclusion.

Appendix E: Integration of Directional Emissivity for Alumina

The spectral hemispherical emissivity was obtained by integration of the spectral directional emissivity data from Ref. [12] as

ε (λ, T) = \int_{0}^{π / 2} ε' (θ, λ, T) sin (2 θ) d θ

(E1)

The spectral directional emissivity data are in 12 deg increments of polar angle θ from 0 deg to 72 deg. In all cases, the data are essentially constant from 0 to 36 deg, and this range was integrated analytically. From 36 deg to 84 deg, a five-point trapezoidal rule was used, and the integral from 84 deg to 90 deg was approximated as a trapezoid. The value at 90 deg was set to zero, in line with theory. This procedure was found to have a numerical truncation error of 1.0% for a gray surface.

The data for the lowest emissivities showed angular behavior consistent with a dielectric medium, as noted in Ref. [12], decreasing at higher angles. On this basis, the point at 84 deg was interpolated, using a value representative of variations at large angle for a dielectric: ε(84 deg, λ) ≈ 0.75 ε(72 deg, λ). Without large angle measurements or the complex refractive index, we cannot exclude the presence of the kind of peak emissivity above 80 deg predicted by Drude's model for metals [13]; and, in particular, the data at 17 and 23 μm do show a 30–50% increase in emissivity by 72 deg. Even so, a sensitivity analysis letting ε(84 deg, λ) ≈ 2.5 ε(72 deg, λ) increases the hemispherical emissivity by only about 5% of the previous estimate. Absent further data, there is not much basis upon which to adjust the calculation further.

Appendix F: Nondimensional Formulation of Model Surface Heat Exchange

Figures 4 and 5 plot heat exchange for a surface at T₁ = T_z that switches between perfectly black and perfectly reflective behavior at wavelength λ_z, where X_z = c₂/λ_zT_z = 3.92069, for the case T₁ = T_z = 300 K. The four equations used may be nondimensionalized as follows:

exact : \frac{q_{net}}{σ T_{1}^{4}} = ε (T_{z}) - α (T_{z}, T_{2}) {(\frac{T_{2}}{T_{1}})}^{4}

(F1)

gray : \frac{q_{net}}{σ T_{1}^{4}} \approx ε (T_{z}) [1 - {(\frac{T_{2}}{T_{1}})}^{4}]

(F2)

ε^{i} (T_{z}) : \frac{q_{net}}{σ T_{1}^{4}} \approx 4 ε^{i} (T_{z}) \frac{Δ T}{T_{1}}

(F3)

ε^{i} (T_{m}) : \frac{q_{net}}{σ T_{1}^{4}} \approx 4 ε^{i} (T_{m}) {(\frac{T_{m}}{T_{1}})}^{3} \frac{Δ T}{T_{1}}

(F4)

In all cases, $q_{net} / σ T_{1}^{4} = fn (T_{2} / T_{1})$ ⁠, since T_m/T₁ and ΔT/T₁ may both be written in terms of T₂/T₁.

Surface Black Below λz.

From Eqs. , , and with α_sw = 1, Δα = −1, and λ_c = λ_z, we have

ε (T_{z}) = 0.4177 and ε^{i} (T_{z}) = 0.6017

(F5)

ε^{i} (T_{m}) = 1 - [\frac{90}{π^{4}} ζ (X_{c, m}, 4) - F (X_{c, m})]

(F6)

α (T_{z}, T_{2}) = 1 - \frac{90}{π^{4}} ζ (X_{c, 2}, 4)

(F7)

where $X_{c, m} = c_{2} / λ_{z} T_{m} = X_{z} (T_{1} / T_{m}) = 3.92069 (T_{1} / T_{m})$ and $X_{c, 2} = c_{2} / λ_{z} T_{2} = X_{z} (T_{1} / T_{2}) = 3.92069 (T_{1} / T_{2})$ ⁠.

The curves for this case are plotted in Fig. 10. Accuracy of the linearization-based εⁱ(T_m) is excellent over the range T₂/T₁ ≈ 1 ± 0.30 in even this extreme case of nongray behavior. The gray approximation strongly overpredicts q_net when T₂ > T₁ because it underestimates the absorption of irradiation.

Surface Black Above λz.

From Eqs. , , and with α_sw = 0 and Δα = 1, we have

ε (T_{z}) = 0.5823 and ε^{i} (T_{z}) = 0.3983

(F8)

ε^{i} (T_{m}) = [\frac{90}{π^{4}} ζ (X_{c, m}, 4) - F (X_{c, m})]

(F9)

α (T_{z}, T_{2}) = \frac{90}{π^{4}} ζ (X_{c, 2}, 4)

(F10)

The curves for this case are plotted in Fig. 11. As before, the linearization-based εⁱ(T_m) is very accurate over the range T₂/T₁ = 1 ± 0.30. The gray approximation strongly underpredicts q_net when T₂ > T₁ because it overestimates the absorption of irradiation.

The internal emissivity at T_m is shown for both surfaces in Fig. 12.

References

1.

Edwards

,

D. K.

,

1969

, “

Radiative Transfer Characteristics of Materials

,”

ASME J. Heat Transfer

,

91

(

1

), pp.

1

–

15

.

Google Scholar

Crossref

2.

Czerny

,

M.

, and

Walther

,

A.

,

1961

,

Tables of the Fractional Functions of the Planck Radiation Law

,

Springer-Verlag

,

Berlin

.

3.

Edwards

,

D. K.

, and

Nelson

,

K. E.

,

1961

, “

Maximum Error in Total Emissivity Measurements Due to Non-Grayness of Samples

,”

ARS J.

,

31

(

7

), pp.

1021

–

1022

.

4.

Edwards

,

D. K.

,

1970

, “

Thermal Radiation Measurements

,”

Measurement Techniques in Heat Transfer

(AGARDograph, Vol. 130),

E. R. G.

Eckert

and

R. J.

Goldstein

, eds.,

Advisory Group for Aerospace Research and Development of NATO, Technivision Services

,

Slough, UK

, Chap. 9.

5.

Edwards

,

D. K.

,

1981

,

Radiation Heat Transfer Notes

,

Hemisphere Publishing

,

New York

.

6.

Edwards

,

D. K.

,

Denny

,

V. E.

, and

Mills

,

A. F.

,

1979

,

Transfer Processes

,

2nd ed.

,

Hemisphere Publishing

,

Washington, DC

.

7.

Mills

,

A. F.

,

1999

,

Heat Transfer

,

2nd ed.

,

Prentice Hall

,

Upper Saddle River, NJ

.

8.

Rosseland

,

S.

,

1936

,

Theoretical Astrophysics

,

Oxford University Press

,

London

.

9.

Brewster

,

M. Q.

,

1992

,

Thermal Radiative Transfer and Properties

,

Wiley

,

New York

.

10.

Modest

,

M. F.

,

1993

,

Radiative Heat Transfer

,

McGraw-Hill

,

New York

.

11.

Lienhard

,

J. H.

, IV

, and

Lienhard

,

J. H.

, V

,

2018

,

A Heat Transfer Textbook

,

4th ed.

,

Phlogiston Press

,

Cambridge, MA

.

12.

Teodorescu

,

G.

, and

Jones

,

P. D.

,

2008

, “

Spectral and Directional Emittance of Alumina at 823 K

,”

J. Mater. Sci.

,

43

(

22

), pp.

7225

–

7229

.

Google Scholar

Crossref

13.

Baehr

,

H. D.

, and

Stephan

,

K.

,

1998

,

Heat and Mass Transfer

,

Springer

,

Berlin

.

14.

Whittaker

,

E. T.

, and

Watson

,

E. N.

,

1920

,

A Course of Modern Analysis

,

3rd ed.

,

Cambridge University Press

,

Cambridge, UK

.

15.

Euler

,

L.

,

1740

, “

De summis serierum reciprocarum

,” Commentarii Academiae Scientiarum Petropolitanae, Vol. 7, pp. 123–134, First communicated to Daniel Bernoulli in 1734 and read before the St. Petersburg Academy in December 1735.

16.

Euler

,

L.

,

1748

,

Introductio in analysin infinitorum

, Vol. 1,

Apud Marcum-Michaelem Bousquet & Socios

,

Lausannae, Switzerland

.

17.

Olver

,

F. W. J.

,

Olde Daalhuis

,

A. B.

,

Lozier

,

D. W.

,

Schneider

,

B. I.

,

Boisvert

,

R. F.

,

Clark

,

C. W.

,

Miller

,

B. R.

, and

Saunders

,

B. V.

, eds.,

2017

, “NIST Digital Library of Mathematical Functions, Release 1.0.17,” National Institute of Standards and Technology, Gaithersburg, MD, accessed Dec. 22, 2018, http://dlmf.nist.gov/

18.

Galassi

,

M.

,

Davies

,

J.

,

Theiler

,

J.

,

Gough

,

B.

,

Jungman

,

G.

,

Alken

,

P.

,

Booth

,

M.

,

Rossi

,

F.

, and

Ulerich

,

R.

,

2017

, “

GNU Scientific Library Release 2.4

,” Free Software Foundation, Boston, MA, accessed Dec. 12, 2018, https://www.gnu.org/software/gsl/

19.

Pall

,

M.

,

2017

, “FFI Library,” luajit, Muenchen, Germany, accessed Mar. 1,

2018

, http://www.luatex.org/svn/branches/0.79/source/libs/luajit/LuaJIT-2.0.3/doc/ext_ffi.html

20.

Pontifical Catholic University of Rio de Janeiro

,

2017

, “

Lua 5.3 Reference Manual

,” Pontifical Catholic University of Rio de Janeiro, Rio de Janeiro, Brazil, accessed June 18, 2018, https://www.lua.org/manual/5.3/

21.

Pégourié-Gonnard

,

M.

,

2013

, “

A Guide to LuaLaTeX

,” DANTE eV, Heidelberg, Germany, accessed May 5, 2013, http://dante.ctan.org/tex-archive/info/luatex/lualatex-doc/lualatex-doc.pdf

22.

Koch

,

R.

,

2018

, “

TEXShop, Version 4.01

,” University of Oregon, Eugene, OR, accessed Apr. 16, 2018, http://pages.uoregon.edu/koch/texshop/obtaining.html

23.

Feuersänger

,

C.

,

2017

, “

Manual for Package PGFPLOTS, Version 1.15

,” SourceForge Media, La Jolla, CA, accessed June 5, 2017, http://sourceforge.net/projects/pgfplots

24.

Montijano

,

J. I.

,

Pérez

,

M.

,

Rández

,

L.

, and

Varona

,

J. L.

,

2014

, “

Numerical Methods With LuaLaTeX

,”

TUGboat

,

35

(

1

), pp.

51

–

56

.https://tug.org/TUGboat/tb35-1/tb109montijano.pdf

25.

Menke

,

H.

,

2017

, “

Plot Complete Fermi-Dirac Integral in LuaLaTeX

,” Stack Overflow, New York, Nov. 29, 2017, https://tex.stackexchange.com/questions/403766/plot-complete-fermi-dirac-integral-in-lualatex

26.

Menke

,

H.

,

2018

, “

Tutorial: Using External C Libraries With the LuaLaTeX FFI

,”

TUGboat

,

39

(

1

), pp.

37

–

40

.https://www.tug.org/TUGboat/tb39-1/tb121menke-ffi.pdf

27.

Dyer

,

R.

,

2015

, “

Compiling the GSL Library for OSX

,” Virginia Commonwealth University, Richmond, VA, accessed July 29, 2015, https://dyerlab.ces.vcu.edu/2015/07/29/compiling-the-gsl-library-for-osx/

28.

Hardy

,

G. H.

, and

Wright

,

E. M.

,

1979

,

An Introduction to the Theory of Numbers

,

5th ed.

,

Oxford University Press

,

Oxford, UK

.

29.

Press

,

W. H.

,

Teukolsky

,

S. A.

,

Vetterling

,

W. T.

, and

Flannery

,

B. P.

,

1988

,

Numerical Recipes in C

,

Cambridge University Press

,

Cambridge, UK

.

30.

Hornbeck

,

R. W.

,

1975

,

Numerical Methods

,

Quantum Publishers

,

New York

.

Linearization of Nongray Radiation Exchange: The Internal Fractional Function Reconsidered

Introduction

Background

Nongray Errors

External and Internal Radiation Fractions

Difference Between External and Internal Emissivities

Linearization Errors

Comparison of Models

Application to Representative Materials

Alumina.

Drude/Hagens-Ruben Metal.

Spectrally Selective Materials.

Accuracy Versus Effort

Summary

Acknowledgment

Nomenclature

Appendix A: Incomplete Zeta Function

Appendix B: Numerical Implementation

Appendix C: The Constant Xz

Appendix D: Linearizations as Finite Difference Approximations

Appendix E: Integration of Directional Emissivity for Alumina

Appendix F: Nondimensional Formulation of Model Surface Heat Exchange

Surface Black Below λz.

Surface Black Above λz.

References

Contents

Figures & Tables

Supplements

References

Get Email Alerts

Cited By

Journals

Conference Proceedings

eBooks

Resources

Opportunities

X_z	X_c	X_c_,_m	X_c_,2
3.92069	5.70943	6.32429	7.08757
ε(T₁)	εⁱ(T₁)	εⁱ(T_m)	α(T₁, T₂)
0.7258	0.6237	0.6807	0.7964
q_gray	$q_{int, T_{1}}$	$q_{int, T_{m}}$	q_exact
400.2	462.1	371.0	371.8

Linearization of Nongray Radiation Exchange: The Internal Fractional Function Reconsidered

Introduction

Background

Nongray Errors

External and Internal Radiation Fractions

Difference Between External and Internal Emissivities

Linearization Errors

Comparison of Models

Application to Representative Materials

Alumina.

Drude/Hagens-Ruben Metal.

Spectrally Selective Materials.

Accuracy Versus Effort

Summary

Acknowledgment

Nomenclature

Appendix A: Incomplete Zeta Function

Appendix B: Numerical Implementation

Appendix C: The Constant Xz

Appendix D: Linearizations as Finite Difference Approximations

Appendix E: Integration of Directional Emissivity for Alumina

Appendix F: Nondimensional Formulation of Model Surface Heat Exchange

Surface Black Below λz.

Surface Black Above λz.

References

Contents

Figures & Tables

Supplements

References

Related

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only