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 [57] 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, α(λ, T1), of a surface at T1, is related to the total hemispherical2 emissivity, ε(T1), by integration over the Planck distribution of blackbody monochromatic emissive power, eλ,b(T1) 
ε(T1)=1σT140α(λ,T1)eλ,b(T1)dλ
(1)
 
eλ,b(T1)=πc1λ5[exp(c2/λT1)1]
(2)
where c1 and c2 are the first and second radiation constants. If this surface is irradiated by an effectively black environment at T2 (for example, if the surroundings are much larger than the surface), the total hemispherical absorptivity is 
α(T1,T2)=1σT240α(λ,T1)eλ,b(T2)dλ
(3)
The net radiation leaving this surface is 
qnet=σε(T1)T14σα(T1,T2)T24
(4)

Nongray Errors

For a gray surface, α(λ, T) is independent of wavelength and so ε(T1) = α(T1, T2). However, some older literature went farther, applying this idea to nongray surfaces for which T1 is near T2. Edwards noted that the fact that limT2T1α(T1,T2)=ε(T1) 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 qnet 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. (4) may be linearized 
α(T1,T2)T24α(T1,T1)T14+ddT2(α(T1,T2)T24)|T1(T2T1)
(5)
 
=ε(T1)T14+4T13[ε(T1)+T14dαdT2|T1](T2T1)
(6)
Thus 
qnet4σT13[ε(T1)+T14dαdT2|T1](T1T2)
(7)

For a gray (or black) surface, /dT2 = 0, and we recover a common linearization: qnet4σεT13ΔT [11]. Thus, the purpose of introducing the internal fractional function is to account for the difference in slope when /dT2 ≠ 0.

External and Internal Radiation Fractions

The fraction of blackbody radiation between wavelengths of 0 and λ is 
f(λT)=1σT40λeλ,bdλ
(8)
 
=115π40c2/λTt3et1dt
(9)
 
=190π4ζ(c2/λ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, fi(λT), arises from a linearization that allows ε(Ts) and α(Ts,Te) to be combined, as described next.

Edwards [1] defines the internal hemispherical emissivity such that 
εi(T1)limT2T1σε(T1)T14σα(T1,T2)T24σT14σT24
(11)
Because the numerator and denominator both go to zero in the limit, L'Hôpital's rule (or Taylor expansion) is required 
εi(T1)=limT2T10α(λ,T1)T2eλ,b(T2)dλ4σT23
(12)
 
=0α(λ,T1)[14σT13T2eλ,b(T2)|T1]dλ
(13)
Thus, when T2 is not too much different from T1 
qnetεi(T1)4σT13(T1T2)
(14)
with 
εi(T)=14σT30α(λ,T)eλ,bTdλ=01α(λ,T)dfi(λT)
(15)
where the internal radiation fractional function is 
fi(λT)14σT30λeλ,bTdλ
(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 (15) stands in contrast to the external total hemispherical emissivity, which from Eqs. (1) and (8) maybe written 
ε(T)=1σT40α(λ,T)eλ,b(T)dλ=01α(λ,T)df(λT)
(17)
From these relationships, one can show that 
fi(λT)f(λT)=F(X)=154π4X4eX1
(18)

where Xc2/λT. The functions f, fi, 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=01α(λ,T)df(λT)01α(λ,T)dfi(λT)
(19)
 
=0α(λ,T)dFdXdX
(20)
 
=0Xzα(λ,T)dFdXdX+Xzα(λ,T)dFdXdX
(21)
where X = Xz is the finite zero point of dF/dX, Xz = 3.92069 (see Appendix  C). The function dF/dX is positive for X < Xz and negative for X > Xz, and 0<α(λ,T)1. Thus, we can form two bounds on the difference 
εεi0XzdFdXdX=F(Xz)ifεεi>0,and
(22)
 
εiεXzdFdXdX=F(Xz)ifεiε>0
(23)
A simple calculation shows that 
F(Xz)=15π4Xz3eXz=0.18400
(24)
In other words 
|εεi|0.18400
(25)
The value of Xz corresponds to 
(λT)z=14387.773.92069=3699.70μmK
(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. (15) and (17) 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, εi is 44% greater than ε; in the second case, εi is 32% less than ε.

For a gray surface, Eq. (20) gives ε − εi = 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 T2 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 Xz. For example, from the preceding development, for the surface black on short wavelengths 
ε=0λzα(λ)df(λTz)=XzdfdXdX=f(Xz)0
(29)
 
=190π4ζ(Xz,4)=190π4ζ(3.92069,4)
(30)
 
=0.41771
(31)

Equation (18) gives εi, 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 qnet for a black surface. Their ratio is 
4T13ΔTT14T24=4T13(T12+T22)(T1+T2)=4[1+(T2/T1)2](1+T2/T1)
(32)

Equation (32) is plotted in Fig. 3. For T2 = 320 K and T1 = 300 K, the linearized flux is 90.5% of the exact flux. The ratio is asymptotic to 13ΔT/2T1+5(ΔT)2/4T12 as ΔT → 0.

The linearization around T1 is considerably less accurate than the traditional approach (i.e., for radiation heat transfer coefficients) of linearizing around Tm = (T1 + T2)/2. For that case, the ratio is 
4Tm3ΔTT14T24=4Tm3(2Tm2+ΔT2/2)(2Tm)=11+ΔT2/4Tm2
(33)

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

From these considerations, we see that εi will also have a limited range of accuracy if evaluated at T1. Instead, εi should be computed at Tm, following the analysis in Appendix  D: 
εi(Tm)=14σTm30α(λ,T1)eλ,bT|Tmdλ=01α(λ,T1)dfi(λTm)
(34)
with the heat flux evaluated, to an accuracy of O(ΔT3), as 
qnet4εi(Tm)·σTm3ΔT
(35)

Comparison of Models

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

  1. (1)

    The gray surface approximation, ε = α, has the wrong slope as T2T1, consistent with Eq. (7). This approximation has poor accuracy for even the smallest temperature differences. The ratio of gray slope to exact slope is ε/εi. 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 T2 > T1, where irradiation dominates net exchange.

  2. (2)

    The approximation using εi(T1) has the correct slope as T2T1, 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 T2 = 320 K is 13% less than the exact value, Eq. (4). This case corresponds to the seemly small temperature difference ΔT/T1 = −0.067.

  3. (3)

    The approximation using εi(Tm) 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/T1 = ±0.33.

Although Figs. 4 and 5 use 300 K as a reference point, Eqs. (29)(31) show that, for a surface at T1 = Tz that switches from black to reflective behavior at λz, both ε(Tz) and εi(Tz) are fixed. The values of εi(Tm)=εi(TzΔT/2) and α(Tz, Tz − ΔT) depend additionally only upon ΔT/Tz. 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 εi at T1. Edwards et al. [6] did not specify the temperature, but without comment included an example worked using Tm. Likewise, Czerny and Walther [2] provided an example worked using Tm, again without comment. The present results show clearly that only Tm should be used when evaluating εi.

Application to Representative Materials

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

Alumina.

Experimental data for polycrystalline alumina (99.5% Al2O3, 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 fi at 2 μm still has the finite value 0.059; the extrapolation extends the range to fi = 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,εni=0.404,andεniε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 εi(Tm) 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: 
ε(λ,T)=48.70reλ{1+[31.62+6.849ln(reλ)]reλ166.78reλ+}
(37)
for re/λ<5×104Ω·cm/μm, where re(T) is the electrical resistivity in Ω cm.

We consider a surface at T1 = 373 K with re = 13.1 × 10−6 Ω ⋅ cm, a value representative of platinum. At this temperature, 0.01f0.99 for 3.88λ61.3μm and fi = 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 T2 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. (35) using εi(Tm) closely tracks the exact solution over a wide range (down to at least 300 K or ΔT/T1 = 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 
α(λ)={αswforλλcαlwforλ>λc
(39)
From Eqs. (10) and (17) 
ε(T)=αswf(λcT)+αlw[1f(λcT)]
(40)
 
=αsw+90π4Δαζ(Xc,4)
(41)
where Xc = c2/λcT and Δα = αlw − αsw. Similarly, from Eqs. (15) and (18) 
εi(T)=αswfi(λcT)+αlw[1fi(λcT)]
(42)
 
=αsw+Δα[90π4ζ(Xc,4)F(Xc)]
(43)
Further 
ε(T)εi(T)=ΔαF(Xc)
(44)
and also 
εi(Tm)=αsw+Δα[90π4ζ(Xc,m,4)F(Xc,m)]
(45)
where Xc,m = c2/λcTm. Finally, from Eq. (3) 
α(T1,T2)=αsw+90π4Δαζ(Xc,2,4)
(46)
with Xc,2 = c2/λcT2. The previous results show that the impact of selectivity will be greatest when Xc and Xz 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 T1 = 360 K and exchanges radiation with sky at T2 = 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/m2, and Eq. (35), 371.0 W/m2, is excellent, whereas Eq. (14), 462.1 W/m2, performs quite poorly.

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

Accuracy Versus Effort

One rationale for introducing εi 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, ε(T1) ≈ α(T1, T2), 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, εi(Tm) 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 εi 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 
Rrad=14εσT3A
(47)
Rrad is routinely evaluated using a gray body approximation. For small temperature differences, using the internal emissivity, εi(T1), in Eq. (47) will provide a more accurate value than a gray body approximation. The present work has shown εi(Tm) to be far more robust, with excellent accuracy for |ΔT/T1| of 20–30%. Hence, the radiation thermal resistance is best calculated as 
Rrad=14εi(Tm)σTm3A
(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 εi(T1) to linearize nongray exchange for enclosures with modest temperature differences. In this study, the following new findings are reported:

  1. (1)

    The maximum absolute difference between ε(T1) and εi(T1) is 0.18400, Eq. (25).

  2. (2)

    The internal emissivity should be evaluated at the mean temperature, Tm, not T1 as has often been suggested. The differences in accuracy are shown to be significant because the evaluation at Tm using Eq. (35) is effectively a second-order accurate numerical method, with truncation error in the heat flux of O(ΔT3). In contrast, evaluation at T1 is a first-order approximation.

  3. (3)

    Equation (35) with εi(Tm) has excellent accuracy for |ΔT/T1| up to 20–30%. Beyond this range, the linearization error of Fig. 3 will cause increasing error for any surface.

  4. (4)

    Theory and examples for several nongray materials show that the gray-body approximation gives the wrong slope for heat flux as T2T1.

  5. (5)

    Calculations involving both the internal and external fractional functions can be conveniently implemented using the incomplete zeta function.

  6. (6)

    εi(Tm) 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

     
  • c1 =

    first radiation constant, 3.741772 × 10−16 (W m2)

  •  
  • c2 =

    second radiation constant, 14387.77 (μm K)

  •  
  • eλ,b =

    blackbody monochromatic emissive power (W m−3)

  •  
  • f =

    (external) radiation fractional function

  •  
  • F =

    fi − f, Eq. (18)

  •  
  • fi =

    internal radiation fractional function

  •  
  • qnet =

    net radiant heat flux from surface 1 to 2 (W m−2)

  •  
  • re =

    electrical resistivity (Ω cm)

  •  
  • Rrad =

    radiation thermal resistance, Eq. (47) (K W−1)

  •  
  • T =

    temperature (K)

  •  
  • Tm =

    mean temperature, (T2 + T1)/2 (K)

  •  
  • T1 =

    temperature of surface (K)

  •  
  • T2 =

    temperature of environment (K)

  •  
  • X =

    c2/λT

  •  
  • Xz =

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

  •  
  • Y =

    function in Appendix  C, Eq. (D1)

  •  
  • αlwαsw =

    see Eq. (39)

  •  
  • α(T1, T2) =

    total hemispherical absorptivity

  •  
  • α(λ, T) =

    spectral hemispherical absorptivity

  •  
  • δT =

    T2 − T1 (K)

  •  
  • ΔT =

    T1 − T2 (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

  •  
  • εi(T) =

    internal total hemispherical emissivity

  •  
  • εni(T) =

    internal total normal emissivity

  •  
  • ζ(s) =

    Riemann zeta function

  •  
  • ζ(Xs) =

    incomplete zeta function

  •  
  • θ =

    polar angle (rad)

  •  
  • λ =

    wavelength (m)

  •  
  • σ =

    Stefan–Boltzmann constant, 5.670367 × 10−8 (W m−2 K−4)

2

The spectral hemispherical absorptivity and emissivity are equal if the irradiation is diffuse: α(λ, T1) = ε(λ, T1) [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)=1σT40λ2πhco2λ5[exp(hco/kBTλ)1]dλ
(A1)
 
=1σT42πkB4T4h3co2c2/λTt3et1dt
(A2)
When λT, f =1 and the last equation yields the well-known result 
σT4=2πkB4T4h3co20t3et1dtζ(4)Γ(4)
(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) = π4/90 (see also Ref. [16, Sec. 167]), from which we recover the usual definition of the Stefan–Boltzmann constant, σ. Returning to Eq. (A2) with this information, we have 
f(λT)=15π40t3et1dt15π40c2/λTt3et1dt
(A4)
 
=115π40Xt3et1dt
(A5)
 
=115π4Γ(4)ζ(X,4)
(A6)
where X = c2/λT and we have identified the incomplete zeta function, ζ(X, s) [17, Sec. 8.22]. Hence 
f(λT)=190π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, D3(X) = 3 ζ(X, 4) Γ(4)/X3, 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: [2427]. All of these packages (other than Mac OS X) are available online at no cost and are actively maintained.

Appendix C: The Constant Xz

Xz is the finite zero point of dF/dX, specifically the nonzero root of 
4(1eXz)=Xz
(C1)
The equation may be solved by iteration to find Xz=3.92069. Alternatively, by defining y =4 − Xz, Eq. (C1) becomes 
yey=4e4
(C2)
This equation may be solved using the Lambert W function 
Xz=4W(4e4)=3.92069
(C3)

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

Diophantine approximations to Xz may be constructed by calculating a continued fraction representation of Xz through the usual process 
Xz=3.92069=3+11+111+11+11+1
(C4)
Successive convergents (finite truncations) of the fraction give an infinite number of rational approximations to Xz 
Xz{4,4712,,14938,24763,,1137290,}
(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 qnet to be a function of T that we wish to determine approximately at T = T2 by evaluating at other temperatures, i.e., T = T1 and T = Tm. With 
Y(T)=σε(T1)T14σα(T1,T)T4
(D1)
we seek to estimate 
qnet=Y(T2)=σε(T1)T14σα(T1,T2)T24
(D2)
A first-order Euler method starting at T = T1 would approximate Y (T2) based only on conditions at T1 as follows: 
Y(T2)Y(T1)+Y(T1)·δT=Y(T1)·δT
(D3)
for δT = T2 − T1 (= −ΔT), which gives Eq. (7) or Eq. (14). The local truncation error of this single-step approximation is O(YΔT2) [29].3
A second-order Runge–Kutta method works from Tm with expansions toward both T1 and T2, subtracting the former from the latter 
Y(T2)=Y(Tm)+Y(Tm)δT2+Y(Tm)δT28+O(δT3)
(D4)
 
Y(T1)=Y(Tm)Y(Tm)δT2+Y(Tm)δT28O(δT3)
(D5)
 
Y(T2)=Y(T1)+Y(Tm)·δT+O(δT3)
(D6)
 
Y(T2)Y(Tm)·δT
(D7)
This single-step Runge–Kutta approximation has local truncation error O(YΔT3) [29]. In these calculations, α(λ, T1) is evaluated at the surface temperature, T1. With Eqs. (3) and (13) 
Y(Tm)=ddT(σT4α(T1,T))|Tm
(D8)
 
=0α(λ,T1)eλ,bT|Tmdλ
(D9)
 
=4σTm3·εi(Tm)
(D10)
which results in Eq. (35). Given lower truncation error of the Runge–Kutta approach, we should expect significantly better accuracy using Eq. (35) than when using Eq. (7) or Eq. (14); 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)=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 T1 = Tz that switches between perfectly black and perfectly reflective behavior at wavelength λz, where Xz = c2/λzTz = 3.92069, for the case T1 = Tz = 300 K. The four equations used may be nondimensionalized as follows: 
exact:qnetσT14=ε(Tz)α(Tz,T2)(T2T1)4
(F1)
 
gray:qnetσT14ε(Tz)[1(T2T1)4]
(F2)
 
εi(Tz):qnetσT144εi(Tz)ΔTT1
(F3)
 
εi(Tm):qnetσT144εi(Tm)(TmT1)3ΔTT1
(F4)

In all cases, qnet/σT14=fn(T2/T1), since Tm/T1 and ΔT/T1 may both be written in terms of T2/T1.

Surface Black Below λz.
From Eqs. (27), (45), and (46) with αsw = 1, Δα = −1, and λc = λz, we have 
ε(Tz)=0.4177andεi(Tz)=0.6017
(F5)
 
εi(Tm)=1[90π4ζ(Xc,m,4)F(Xc,m)]
(F6)
 
α(Tz,T2)=190π4ζ(Xc,2,4)
(F7)

where Xc,m=c2/λzTm=Xz(T1/Tm)=3.92069(T1/Tm) and Xc,2=c2/λzT2=Xz(T1/T2)=3.92069(T1/T2).

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

Surface Black Above λz.
From Eqs. (28), (45), and (46) with αsw = 0 and Δα = 1, we have 
ε(Tz)=0.5823andεi(Tz)=0.3983
(F8)
 
εi(Tm)=[90π4ζ(Xc,m,4)F(Xc,m)]
(F9)
 
α(Tz,T2)=90π4ζ(Xc,2,4)
(F10)

The curves for this case are plotted in Fig. 11. As before, the linearization-based εi(Tm) is very accurate over the range T2/T1 = 1 ± 0.30. The gray approximation strongly underpredicts qnet when T2 > T1 because it overestimates the absorption of irradiation.

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

References

References
1.
Edwards
,
D. K.
,
1969
, “
Radiative Transfer Characteristics of Materials
,”
ASME J. Heat Transfer
,
91
(
1
), pp.
1
15
.
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
.
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
.