Predicting the Effective Emissivity of an Array of Aligned Carbon Fibers Using the Reverse Monte Carlo Ray-Tracing Method

Taking the next big step in deep-space travel requires innovative and extremely powerful propulsion systems capable of reaching much higher speeds in space. One of NASA's Grand Challenges is developing an ultra-efficient heat radiator for various advanced propulsion and power generation systems, including nuclear-electric and magnetohydrodynamic propulsion [1]. In a constant effort to minimize the mass in space-based systems, replacing the commonly used solid radiator fins with lightweight, porous materials is of great interest. Bulk carbon fiber (e.g., not embedded in a matrix material) has a much lower density and higher thermal conductivity than metals and composites, thus, it has been proposed as a solution to the space radiator challenge [2]. One of the thermal uncertainties of a fibrous radiator is the effect of porosity on the thermal radiation characteristics. Both experimental and numerical evaluation of the relationship between porosity and radiative properties presents certain challenges.

Fibrous media is often characterized as pseudo-continuous, with complex integro-differential radiative heat transfer equations that have only numerical solutions [3,4]. Due to the scattering effects, the radiation intensity at a given point cannot be determined by the local properties and their derivatives like most transport relations (e.g., Fourier's law). Instead, the intensity of radiation at each surface in the domain depends on the intensities at every other surface, which makes radiative heat transfer problems involving porous materials challenging or impossible to solve through numerical integration methods. In the absence of empirical data or suitable analytical approximations of the radiative properties of porous media, a stochastic numerical method called Monte Carlo ray-tracing (MCRT) is the best available method.

Of all of the radiative heat transfer modeling techniques, MCRT can be applied with relative ease to a wide range of geometrically complex problems [5,6]. The MCRT method directly simulates fundamental radiative processes, including emission, absorption, and reflection. The basis of MCRT assumes geometric optics with no diffraction or near-field effects [7]. These assumptions are valid when the features of the domain are large compared to the wavelengths of the photons carrying the thermal energy.

Studies in the literature on radiative heat transfer through fiber insulation have employed the MCRT methods to determine the surface-to-surface view factors, which are useful for solving full radiosity problems [8,9]. Arambakam et al. [9] used an MCRT method to calculate the surface-to-surface view factors needed to predict the temperature gradient across a slice of fibrous insulation to investigate the effects of fiber diameter and packing density on the transmittance of incident radiation through a fiber array [9]. This work used a forward MCRT method where rays were emitted from the hot side of a slice of insulation and traced through the fibrous array.

Other studies have employed MCRT techniques to predict the effective emissivity of cavities and self-viewing bodies, many of which use “reverse MCRT” methods [7–15]. Reverse methods use the principle of reciprocity to obtain solutions for the emissivity of a set of surfaces by measuring the absorptivity of the surfaces. The basic principle of reciprocity uses the assumption that the effective total absorptivity of an isothermal cavity will be equal to the effective total emissivity due to Kirchhoff's Law (⁠$ελ=αλ$ for all materials) and the constitutive equations for all ideal gray materials (⁠$ελ=εT$⁠).

In the reverse method rays are emitted from a simple, typically planar, target surface, and traced back into the domain. This method is especially useful when the emitting surfaces in the domain have complex geometries, are part of a porous medium, or make up the boundaries of a cavity, in which case the reverse method drastically reduces computational cost as compared with the forward method [3,16,17]. Several applications of this method appear in the literature, where the target surface is much smaller or simpler than the emitting surface(s), such as predicting the effective emissivity of blackbody cavities and scattering packed-beds [11,14,15]. One such study by Ishii et al., evaluates the effective emissivity of black body cavities with complex inner geometries using a reverse MCRT method [15].

It should be noted that very few experiments have been conducted on scattering behavior of bulk porous materials because they, too, pose considerable challenges. Zhu et al., however, showed experimentally that porous, optically thick carbon fiber felt is more emissive than a smoother carbon surface [18], presumably due to the increased total surface area exposed to the object's external boundary. This relationship between porosity and increased effective emissivity is of interest to this work.

The objective of this work was to use a framework similar to the aforementioned study by Arambakam to build a reverse MCRT analysis to solve problems where the fibers themselves are the sources of heat. One can intuit that for very low fiber volume fractions, the total emissivity tends to zero, and that as the fiber volume fraction increases, the effective emissivity increases asymptotically toward that of the fibers as the array becomes densely packed. The reverse MCRT method developed here predicts the effective emissivity of a carbon fiber medium made up of aligned, opaque fibers with evacuated void spaces. The interest of this work is to investigate the relationship between the fiber volume fraction and the effective emissivity in the midrange fiber volume fractions, and to determine whether or not this dependence is monotonic.

To develop a better understanding of the effective emissivity and other thermal radiation characteristics of an array of isothermal carbon fibers, we combined the approaches of Arambakam and Ishii with a novel implementation of domain simplification and diffuse ray scattering described in Sec. 2.1.

Due to the geometric complexity, a type of reverse MCRT was the selected modeling technique. The fibers are assumed to be arranged in a parallel, hexagonal-close-packed pattern where the third dimension (z-direction) has no effect on the depth of penetration of radiation, so the model was simplified to two dimensions. In this arrangement, the centers of three nearest-neighbor fibers are placed at the vertices of an equilateral triangle. The primary variable in this study is the fiber volume fraction, which was varied across the entire physical range of nearly zero to the maximum packing density of hexagonal-close-packed cylinders; given by the equation $ρ=(1/6)π3≅0.9069$⁠.

For computational efficiency, the array domain was simplified in the lateral (y-axis) direction such that there are only two half-fibers across the width of the domain and symmetry boundary conditions along lateral edges of the array as shown in Fig. 1. The symmetry boundaries reflect rays specularly, preserving the true depth of ray penetration. For a given fiber volume fraction, the array may be created by fixing the array thickness or by fixing the number of fibers through the array thickness. For a fixed array thickness, the array width and number of fibers (i.e., volume fraction) vary together.

The opaque fibers are assumed to have gray-diffuse surfaces. The emissivity of graphitic materials is typically in the range of 0.7–0.9 depending on the surface preparation and purity. It has been shown experimentally that similar materials (i.e., graphite and carbon–carbon composites), have relatively flat spectral responses over a broad range of temperatures [18,19], and so act as gray bodies. The fiber surfaces have a dull luster so they were assumed to be diffuse reflectors.

In this reverse ray-tracing model, rays are generated continuously across the “source” boundary, from which rays enter the domain. The locations where rays intersect fibers and boundaries within the model domain are calculated directly, one at a time, as described below. Since this method solves for the ray paths directly and continuously over the model domain, no meshes are needed to discretize the radiation source, fibers, nor void spaces.

Each ray is generated at a random x-coordinate along the top boundary and enters the domain at a random angle. The locations at which the rays enter the domain are selected with a uniform probability distribution function (PDF). The angles at which rays enter the domain follow a PDF derived from Lambert's Cosine law for diffuse radiation that was modified for a two-dimensional model. In three-dimensional space, a diffuse surface emits rays uniformly in the azimuthal angle, and with a zenith angle distribution such that the maximum number of rays are emitted in the differential solid angle normal to the surface and no rays are emitted in the differential solid angle tangent to the surface.

For a three-dimensional array of fibers which are infinitely long in the z-direction, the effect of radiation at any z-coordinate is the same. Therefore, the absorption of the whole array may be calculated by integrating all of the rays from a single z-location along the fibers. This integration has the result of reducing the probability distribution function from three dimensions to two dimensions, thus, vastly reducing computational cost. This two-dimensional ray-direction distribution function can be defined by a single altitude angle instead of both azimuthal and zenith angles, as is required for the three-dimensional distribution function. The following is the development of the altitude angle for diffuse radiation in two dimensions for geometries that are uniform in the third dimension.

This PDF was transformed into two-dimensions by integrating all of the rays that strike a fiber along the z-axis for a given x–y location, thereby reducing the spherical coordinates of the ray-direction to a single altitude angle (⁠$θ2$⁠) in the x–y plane, as shown in Fig. 2. This captures the true directional distribution of diffuse hemispherical emission for domains uniform in the z-dimension. The following equation gives the PDF of diffuse radiation for two-dimensional domains that was used in this simulation to generate the directions of incident rays and reflections from fibers:

Figure 3 demonstrates that the PDF of $θ2$ follows Lambert's Cosine Law, which shows that the most probable emission direction is normal to the surface (⁠$θ2$  = π/2), and a zero probability of tangential emission (⁠$θ2$  = 0,π).

This simulation assumes opaque fibers, so radiation is transmitted through only void spaces. In ray-tracing models, there are two common ways to simulate radiation intensity and attenuation. Rays consist of either a single photon or a bundle of photons. When a single-photon ray strikes an absorbing surface, the entire ray is either absorbed or reflected. In this case, the probability that the ray gets absorbed is given by the total absorptivity of the surface, $Pabs=αT$⁠. Thus, the probability of a reflection is $Pref=1−Pabs=1−αT$⁠. For multiphoton rays, each interaction with an absorbing surface reduces the intensity of the ray by a fraction given by the total absorptivity of the surface. The remaining energy is then reflected at angles randomly selected according to the distribution given by Eq. (3) and traced until the fractional energy is reduced to an assumed threshold, or until the ray exits the domain. In general, simulations using multiphoton rays converge with fewer rays than the single-photon method since each ray path contains more information about the radiative behavior of the participating media. Except for the comparison of ray structures briefly discussed in Sec. 3, multiphoton rays were used in this study.

Ray-tracing assumes geometric optics, which is valid when $(πd/λ)>1$⁠, where $d$ is the particle diameter, and $λ$ is the predominant wavelength [13]. According to Planck's law, if the fiber temperature is 1000 K (i.e., the typical operational temperature for the assumed space-based radiator), the fiber diameter must be at least 4 μm for the size parameter criterion to be met by at least 95% of the thermal radiation. The diameter of thermally conductive carbon fiber typically ranges from 10 to 20 μm, well within the range suitable for geometric optics at 1000 K. For smaller-diameter fibers or lower temperatures, less of the thermal radiation would meet the criterion for geometric optics, however, the error in using the assumption is likely insignificant [13].

The model generates a ray and records the ray's path as it travels into the domain until its journey terminates by being absorbed by a fiber, transmitted through the domain, or reflected out of the domain. After thousands of rays are traced, the effective absorptivity is calculated by the ratio of absorbed to incident rays. The effective absorptivity is equivalent to the effective emissivity due to the laws of reciprocity. Figure 4 shows the ray path results for five rays in a model with a fiber volume fraction of 0.127. This demonstrates how each ray has a unique path that may include diffuse reflections from fibers, specular reflections from symmetry boundaries, and absorption at fiber surfaces. For the purposes of clarity, note that no geometrical meshing was required because the ray paths were generated using strictly Monte Carlo techniques, no numerical integration was involved.

Two analyses were conducted for numerical validation of our reverse MCRT model. The first was to use the two-dimensional ray direction PDF given in Eq. (3) to calculate the view factor for a similar problem with a known solution. The second method was to test limiting cases of the full reverse MCRT model for which the solutions are known.

For the first validation analysis, the diffuse radiation PDF used for generating ray emission and reflection directions in two dimensions was tested by using our MCRT algorithm to solve for the known view factor between two infinitely long cylinders with the same diameter, as shown in Fig. 5. Equation (4) is the view factor relation for these cylinders. The value of the view factor obtained using this MCRT method converged to within 2% of the analytical solution in approximately 50,000 rays, as shown in Fig. 6. This confirms that the two-dimensional model generates rays with the correct distribution of angles, and that the ray-tracing is executed properly

For the second validation analysis, the effective emissivity solution of the full fiber array was confirmed by testing our model at limiting cases of fiber solid volume fraction (SVF). Two limiting cases used for validation included an array with the maximum fiber packing density and one with a very low fiber packing density (i.e., more than ten fiber diameters between adjacent fibers). As high fiber packing density tends toward the maximum density, less radiation penetrates the array. Very high fiber densities diminish the scattering effects, and the effective emissivity tends toward that of a flat surface made of the carbon material (i.e., the fiber emissivity). Oppositely, as the fiber packing density approaches zero, the effective emissivity tends toward zero as most of the radiation passes through the array without intersecting any fibers. The results from these two limit tests were zero and approximately the surface emissivity of the individual fibers as expected.

The dependence of the array's effective emissivity on the fiber solid volume fraction (SVF) and the surface emissivity of the individual fibers were determined for a 500 μm-thick array of fibers. The fiber diameter was assumed to be 14-μm, the measured value typical of carbon fiber manufactured by Mitsubishi used in corresponding experimental studies [21]. The fiber SVF in the models ranged from a minimum of 0.0011, which corresponds to two fibers through the array cross section, to the physical maximum of 0.9069, which corresponds to 38 close-packed fibers through the array cross section.

Figure 7 gives results of the total fraction of incident energy absorbed by, transmitted through, and reflected from the array versus solid volume fraction for a set of simulations that assume a surface emissivity of the individual fibers of 0.8, an array thickness of 0.5 mm, and use 50,000 multiphoton rays. Again, the emissivity is equal to the absorptivity. As expected, transmittance decreases monotonically with increasing fiber volume fraction in a smooth manner similar to standard attenuation curves. Reflectivity increases linearly with fiber volume fraction above 0.01, when the fraction of reflected rays becomes statistically significant. The effective emissivity does not increase monotonically with the fiber volume fraction. Figure 8 shows results from a set of simulations similar to Fig. 7 except the array thickness was varied in addition to the fiber volume fraction. The array thickness was varied by specifying the number of rows of fiber in the domain. Figure 8 shows that the effective emissivity of the array varies nonmonotonically with fiber volume fraction, as was demonstrated in Fig. 7, but increases monotonically with array thickness, as more fibers will always result in more ray–fiber interactions and fewer reflected and transmitted rays.

For very small fiber volume fractions, less than about 0.1, the rate of transmission dominates since the array is optically thin. As the solid volume fraction increases, transmittance decreases quickly and the reflectivity increases steadily as an increasing fraction of incident rays reflect out of the domain from the top rows of fiber. Transmittance drops below 0.10 and 0.01 for solid volume fractions above about 0.05 and 0.15, respectively, thus, volume fractions above about 0.15 have optical thicknesses less than 500 μm.

The effective emissivity rises from zero to a maximum at a volume fraction near 0.16 due to the reduction in transmission with an increasing volume fraction as rays strike more fibers as they travel through the domain. After the maximum, the effective emissivity decays asymptotically toward the value of the solid surface emissivity. For optically thick arrays, the effective emissivity is greater than the surface emissivity of the individual fibers. This is due to scattering events that allow radiation emitted from surfaces beneath the top row of fibers to exit the domain after reflecting from other fiber surfaces. This scattering effect yields a greater view factor between the array and its surroundings as compared with a flat emitting surface. Optical thickness is a function of fiber diameter and fiber volume fraction, so for a specific size of fiber, optical thickness can be achieved by increasing the volume fraction and/or increasing the array thickness to at least the optical thickness. Thus, as the array thickness decreases the fiber volume fraction that yields the maximum effective emissivity increases to optimally balance the transmission and reflection, as the volume fraction 0.16 does for the 500 μm thick array of 14 μm-diameter fibers.

It was observed that the peak in effective emissivity occurs at the same fiber volume fraction at which point there is no significant direct transmission (i.e., transmission without striking a fiber). With diffuse radiation entering the array at all hemispherical angles, the only point at which direct transmission ceases is when all fibers touch nearest neighbors. However, since most diffuse radiation from a flat surface is oriented nearly normal, the rate of transmission through the array drops dramatically when there is no clear path through the fiber array in the direction normal to the plane of incident radiation. To demonstrate that eliminating direct transmission corresponds with peaking effective emissivity, a simulation using collimated radiation was run. With collimated radiation, there is a single fiber volume fraction above which there is no direct transmission. The results for the collimated radiation simulation are given in Fig. 9, demonstrating that the peak absorptivity occurs at the volume fraction (0.2267) at which point there is no direct transmission from the normal direction. This volume fraction balances low transmittance with low reflectivity. As with normal incident radiation, each angle of incidence will have a unique critical spacing value that allows direct transmission through the array. Thus, for diffuse incident radiation, the effective emissivity curve has a broad peak due to the contributions of the different angles of incidence.

Comparing these results with the real-world conditions, the fiber volume fraction was estimated for as-received bundles of carbon fiber. A sample of Mitsubishi carbon fiber tows with 2000, 14 μm-diameter fibers per tow was approximately 3 mm wide and 500 μm thick, which yields a fiber volume fraction of about 20%. Since the true fiber volume fraction of the Mitsubishi carbon fiber tows was very close to the optimum volume fraction, these tows naturally have a near-optimal spacing for maximum effective emissivity.

Next, the effects of surface emissivity of the individual fibers are discussed. The total fraction of energy absorbed increases with the emissivity of the fiber surfaces; each time a ray strikes a fiber, more of its energy is absorbed. This is demonstrated by the results shown in Fig. 10, which gives the effective emissivity versus fiber volume fraction for fiber surface emissivities of 0.6, 0.7, 0.8, and 0.9 for a 500 μm-thick array.

The effective emissivity of the array is greater than the surface emissivity of the individual fibers for fiber volume fractions above about 0.05. The maximum effective emissivity is proportional to the surface emissivity of the individual fibers and is generally greater than the surface emissivity by slightly less than half of the difference between the fiber emissivity and one, as shown in Fig. 11. By the constitutive relation, this also implies that the total reflectivity of an optically thick array is proportional to the reflectivity of an individual fiber.

At fiber volume fractions greater than about 0.16, where the maximum typically occurs, the effective emissivity curves decrease again, approaching the fibers' surface emissivity. These curves diverge slightly as the total reflectivity of the array becomes increasingly significant. Figure 12 gives the reflectivity of the array versus fiber volume fraction for the same set of simulations, and shows that as the slope of the reflectivity curve increases the fiber emissivity decreases. For lower fiber emissivities, more energy is reflected after each fiber interaction, which increases the sensitivity of the total reflectivity to the fiber volume fraction for the large volume fractions. For optically thick arrays of fibers with a lower emissivity (i.e., higher reflectivity), the effective emissivity of the array is more sensitive to how much internal scattering occurs as compared with higher emissivity cases since more internal scattering increases the total amount of energy that is ultimately emitted from the array.

At the upper limit of the surface emissivity of the individual fibers (i.e., $εf=1$⁠, $ρf=0$⁠), the effective emissivity curve also reaches unity at the fiber volume fraction at which transmission is zero and remains at 1.0 for all greater volume fractions. Finally, Fig. 13 shows that the transmittance versus volume fraction is nearly identical for all fiber emissivities. This similarity arises because most of the transmitted energy is due to the direct transmission of rays due to geometrical factors, rather than being transmitted after a reflection from a fiber. Thus, there is a weak correlation between reflectivity and the surface emissivity of the individual fibers for emissivities greater than 0.5.

Next, results from different types of rays are discussed. The type of ray does not impact the simulation result as long as the solution has converged. Simulations using multiphoton rays required at least 50,000 rays to achieve convergence, while simulations using single-photon rays required at least 100,000 rays. The percent error between the absorptivity results from the two ray types was less than 0.05% for all volume fractions.

One benefit of single-photon rays is that they allow for simple monitoring of intensity attenuation. With single-photon rays, the attenuation versus array depth can be determined by recording the location at which each ray is absorbed, versus multiphoton rays that get absorbed through the multiple fiber interactions. Figure 14 shows the percent of incident energy absorbed versus array depth for a wide range of fiber volume fractions. The stepped nature of these curves is due to the fact that the radiation is not absorbed continuously, but only at depths containing fibers. For low fiber volume fractions where there is only a single fiber at a given depth, there is a distinct step in the attenuation curves in Fig. 14 corresponding to the presence of that fiber. For optically thick arrays, the attenuation curve reaches 100% before the bottom boundary of the array.

While both single- and multi-photon rays can be used for reverse MCRT simulations, each has costs and benefits that make it uniquely appropriate for collecting different types of data.

This study used a reverse Monte Carlo ray-tracing method developed for two-dimensional models to prove that the porosity of an array of aligned carbon fibers can dramatically affect the array's effective emissivity. A novel formula was presented for generating two-dimensional rays that fully represents three-dimensional diffuse radiation for a geometry with two-dimensional symmetry. For optically thick arrays, the effective emissivity is greater than the surface emissivity of the individual fibers for volume fractions above about 0.05 and fiber surface emissivities above 0.5. The maximum effective emissivity is generally greater than the fiber emissivity by slightly less than half of the difference between the fiber emissivity and unity. This maximum emissivity occurs at a volume fraction of 0.16 for 500 μm thick arrays. For volume fractions above the point at which the maximum effective emissivity occurs, the effective emissivity tends toward the surface emissivity of the individual fibers. Even at the maximum fiber packing density, the effective emissivity is greater than the surface emissivity of the individual fibers due to scattering. These methods and results are applicable to a wide range of problems involving radiative exchange among arrays of cylinders, including tube-and-shell heat exchangers, fibrous insulation, and other advanced cooling systems.

This work was supported by a Space Technology Research Fellowship granted by The Office of the Chief Technologist, National Aeronautics and Space Administration (Grant No. NNX11AM70H). Their support is gratefully acknowledged.

d =

fiber diameter

d_a =

array depth

$F1−2$ =

view factor from surface 1 to 2

$Pabs$ =

probability of ray absorption

$Pref$ =

probability of a reflection

PDF =

probability density function

r =

fiber radius

R =

random number between 0 and 1

x =

size parameter

$αT$ =

total hemispherical absorptivity

$αλ$ =

spectral absorptivity

$εT$ =

total hemispherical emissivity

$ελ$ =

spectral emissivity

$εf$ =

total hemispherical emissivity of fiber

$θ$ =

Zenith angle

λ =

wavelength

$φ$ =

azimuthal angle