Second-Law Considerations in Monte Carlo Ray-Trace and Discrete Green’s Function Analysis of Coupled Radiation and Conduction Heat Transfer

Ritchey-Chrétien telescopes, such as the one shown in Figs. 1 and 2, are commonly used to monitor the planetary energy budget from low Earth orbit [1–3]. Light from the Earth scene enters through an aperture, passes through a collimator, is reflected from two mirrors, and is finally collected by a thermal radiation detector, as illustrated in Figs. 1 and 2(c). Unfortunately, the detector is also sensitive to thermal emission from the telescope structure, which means that thermal stability of the telescope while in orbit is essential. The common approach to assuring thermal stability involves the use of strategically placed feedback-controlled electric heaters. In this case, a high-fidelity coupled radiation/conduction heat transfer model—a thermal “digital twin” of the instrument—is required to design and optimize the thermal control strategy.

The inherent geometric complexity of space-borne telescopes means that several thousand volume elements are required to achieve meaningful spatial and temporal resolution. Preliminary steps toward building thermal performance models of the instrument depicted in Figs. 1 and 2 have been reported elsewhere [4,5]. Note that only the radiatively active inner surfaces are depicted in Fig. 2. A cylindrical casing, represented by dashed lines in Fig. 2(c), has been removed, and wall thicknesses are not shown. The outer nonradiating surfaces of each telescope component are considered to be insulated; therefore, thermal interaction of the telescope with its surroundings is limited to radiation entering and leaving through the entrance aperture.

is generally violated. Previous studies [9–11] have shown that violation of the second law in formulations involving only radiation heat transfer can translate to uncertainties of the order of a few percent. However, the present contribution establishes that incorporation of the resulting RDFs into a time-dependent finite volume heat conduction analysis results in an instability whose gravity depends on the degree and spatial distribution of second-law violations.

For the instrument depicted in Figs. 1 and 2, when the interior surface temperatures are initialized at a uniform value of 311 K and the telescope is allowed to view a 311-K blackbody positioned in the entrance aperture, all temperatures are expected to remain constant. However, when computed using the coupled-mode numerical model developed in Sec. 2, the element temperatures are observed to vary with time. The curves representing the 2636 temperature histories, plotted in Fig. 3, gather into three groups. This is because, due to the large thermal diffusivity of the aluminum walls, each of the three components of the telescope is quasi-isothermal at a given instant in time. The uppermost group of curves represents the temperature histories of the volume elements constituting the end cap, while the lowest group represents the temperature histories of the elements making up the radiatively floating secondary mirror holder. The horizontal curve represents the fixed temperature of the entrance aperture (the “scene”) as well as that of the isothermal detector, while the gently down-sloping curve regroups the temperatures of all volume elements representing the collimating baffle. The observed temperature variations with time of the end cap, mirror holder, and baffle clearly constitute a violation of the second law of thermodynamics. The model confirms that the mass-averaged temperature of all 2636 volume elements remains constant at 311 K, as required by the first law of thermodynamics. In the absence of early temperature gradients, the solution is initially radiation-dominated. System equilibrium is approached asymptotically as the solution becomes increasingly conduction-dominated and the importance of radiation diminishes. If the distribution factors were exact, Eq. (2) would hold for all combinations of $i$ and $j$⁠, and all temperatures would remain fixed at 311 K. The approximately three-degree temperature spread as the system approaches equilibrium is a measure of the degree of violation of the second law.

The benefits and methods of correcting, or “smoothing,” RDFs obtained using the MCRT method have been previously addressed by Loehrke et al. [9] and by Duan et al. [10]. Loehrke and his coworkers represent the smoothed $AD$ product as a weighted average of $AiDij$ and $AjDji$⁠, where $Dij$ and $Dji$ are estimated using the MCRT method, where weights are based on knowledge of the fractional error bounds in a binomial probability distribution. Duan and his coworkers used a constrained maximum likelihood approach. Taylor et al. [11] have pointed out similar benefits of smoothing view factors (as opposed to RDFs) when only approximate values are available, as might be the case when reading them from a graph. In all three cited contributions, the authors report smoothing strategies that, when implemented for pure radiation problems, improve the accuracy by from 3% to 5%. None of the cited contributions deal with the case considered here of massive time-dependent couple radiation/conduction problems.

In Eq. (3), $Tnp$ and $Tnp+1$ are the temperatures of volume element $n$ at successive time steps $p$ and $p+1$⁠, $Δt$ is the length of the time-step, $Km,n$ is the thermal conductance between each $m$⁠,$n$-pair of neighboring volume elements, $Dm,n$ is the RDF from surface element $m$ to surface element $n$⁠, $Vn$⁠, and $An$ are the volume and surface areas of element $n$⁠, $εm$⁠, and $εn$ are the diffuse emissivities of surface elements $m$ and $n$⁠, $hn$ is the conductance between surface element $n$ and the surroundings (this could be a convection coefficient where applicable or a physical shunt conductance), $ρ$ and $c$ are the mass density and specific heat, $σ$ is the Stefan–Boltzmann constant, $Pj$ is the power emitted by a source in the scene viewed by the telescope, $gnp$ is the heat generation source during time-step $p$ within volume element $n$⁠, and $q″np$ is the heat flux applied (for example, by an electrical heating element in the current application) during time-step $p$ to the surface of volume element $n$⁠. Not all of these terms are present in a given problem.

where $Ni$ is the number of equal-power rays emitted from surface element $i$ during the ray-trace, and $Nij$ is the number of those rays absorbed by surface element $j$ due to direct radiation and to all possible reflections within the enclosure.

Under the assumption of constant thermophysical properties, the elements of the $G$ matrix are invariant with time. The functional similarities between the RDF and the DGF are remarkable, making them very compatible with this application. Both are a direct measure of the sensitivity of the thermal state of one element to the thermal state of the other elements of a system. Also, energy balance expressions based on both are explicit.

where $GT$ is the transpose of the DGF matrix. This approach isolates the nonlinearity difficulties normally associated with coupled conduction/radiation formulations so that they can be systematically explored and controlled.

Early attempts to implement this formulation in the presence of electrical wall heating seemed to work as expected. However, as illustrated in Fig. 3, when the heat sources are removed, the important role played by RDF reciprocity in coupled conduction/radiation problems emerges. Evidently, the power input associated with sufficiently strong heat sources dominates the power associated with second-law violations to the extent that the related errors, while significant, are not apparent.

One way to minimize second-law-related heat transfer errors would be to increase the fineness of the mesh with proportional increases in the number of rays traced. This would tend to isolate and localize the second-law violators, surrounding them with law-abiding neighbors that moderate the violators’ felonious behavior through heat conduction. However, this approach un-necessarily increases the demand for already strained computer resources.

in Eq. (14) lends credibility to the approach. A similar approach has been studied and reported by Loehrke et al. [9] for the case of purely radiative exchange.

When used together, Eqs. (14) and (16) assure that the corrected distribution factor matrix obeys both the first and second laws of thermodynamics. This approach, which is similar to the approach suggested by Taylor et al. [11] for view factors, can be criticized because nothing prevents it from producing negative values of $DiiC$⁠, a situation that clearly has no physical meaning. The authors of Ref. 11 did not encounter this situation, perhaps because their example problem consisted of only six surfaces. While the novel approach presented in Sec. 3, when applied to the massive problem considered here, can also produce negative values of $Dii$⁠, experience shows that they have a negligible impact on the solution. This may be due to the relatively small number of negative RFDs in a large population that otherwise conforms to the first and second laws of thermodynamics, but also to the fact that the magnitudes of on-diagonal distribution factors in the problem at hand are small relative to the magnitudes of the off-diagonal elements. The influence of negative values of on-diagonal RDFs in highly reflective enclosures needs further investigation.

For radiation distribution factors between any two surfaces, the second law of thermodynamics leads to Eq. (2). The usual approach to determining $Dij$ and $Dji$ involves tracing a large number of rays from all $n$ surface elements. Due to the statistical nature of the process, the reciprocity relationship, Eq. (2), is only approximately satisfied, sometimes with a large discrepancy. With this in mind, we propose the following algorithm:

1. Compute both $Dij$ and $Dji$ using the MCRT method. Then accept as more accurate the one that has the smaller value of $εA$⁠. This differs significantly from the approach in Ref. 11 in which the lower off-diagonal view factors are used, along with the reciprocity expression for view factors, to compute the upper off-diagonal view factors.

2. Assign the other RDF based on reciprocity; e.g., $Dij=(εjAj/εiAi)Dji, where εjAj<εiAi$⁠. This is simpler than using Eqs. (11) and (12). The second law is now automatically obeyed.

3. Use the energy conservation principle, $Dii=1−∑j≠inDij$⁠, to determine the corrected diagonal elements. Now the first law is also automatically obeyed. Results reported in Sec. 4 indicate that the occasional negative value of $Dii$ has a negligible impact on the time-dependent temperature field for the absorption-dominated enclosure under consideration.

The justification for this strategy is established from consideration of the following primitive examples.

Example 1. A convex surface surrounded by a larger concave surface

Clearly, the row-1 values, $D11$ and $D12$⁠, are always correct. However, the row-2 values, $D21$ and $D22$⁠, are generally wrong. Equation (18a) is true only in the limit as $A1≪A2$⁠, while Eq. (18b) can only approach the truth as $A1$ approaches $A2$⁠.

Clearly, the row-1 values, $D11$ and $D12$⁠, are again correct in all three outcomes. However, the row-2 values, $D21$ and $D22$⁠, are once again generally wrong. Equation (19a) is true only in the limit as $A1≪A2$⁠, Eq. (19b) is only true in the limit as $A1$ approaches $A2$⁠, and Eq. (19c) is only true when $A1=12A2$⁠.

We can continue this line of thought for three rays, four rays, and so forth, but the conclusion is always the same. By inductive reasoning, the row corresponding to the smaller area is always the more accurate. We conclude that for this geometry the suggested algorithm yields the best results.

Example 2. Infinite parallel plates

The previous example assumed black surfaces and conclusions were made as a result of area differences. How do different emissivity values come into play? To answer this question, consider the infinite parallel plate geometry shown in Fig. 5, where $A1=A2$⁠. We arbitrarily assume that $ε2=1$ and $ε1<1$⁠, so that $ε1<ε2$⁠.

In the concentric-sphere example, $A1A2$ was a surrogate for $A1ε1A2ε2$⁠, and in the parallel-plate example, $ε1ε2$ is a surrogate for the same quantity. We conclude that the logic leading to the observation that the suggested algorithm yields the best results for concentric spheres is equally applicable—word for word—to the case of parallel plates. In fact, we can generalize the results for any number of two-surface geometries for which $A1ε1<A2ε2$⁠.

The distribution factor $Dij$ never appears by itself in the energy equation; rather the combination $εiAiDij$ is always used. The conclusions drawn from both of these examples are identical. Suppose that $A1ε1≪A2ε2$⁠. According to the proposed strategy, we would compute $D11$ and $D12$ statistically using the MCRT method. Unlike the examples shown, these values would be approximate. We then assign rather than statistically estimate $D21$ using the reciprocity relation, Eq. (2), and assign $D22$ using the energy conservation principle, Eq. (1). Since $A1ε1A2ε2<1$⁠, the error in the assigned value of $D21$should be less than the error in the MCRT estimated value of $D21$⁠.

The MCRT engine, which takes advantage of both vectorization and parallelization inherent to the MCRT method, requires several hours of run time on a multicore PC just to produce the RDF matrix; and then the finite volume solver requires another couple of hours to generate the DGF matrix. However, once the DGF matrix has been obtained, the time-marching computation of the transient temperature field requires mere seconds of CPU time. This means that the approach described here is well-suited for optimization studies in which the initial and boundary conditions and electrical heating strategies—indeed, any quantity appearing in Eq. (7)—can be varied without the need to recalculate the DGF matrix.

Figure 3 shows the temperature history of all 2636 volume elements after they were initialized at 311 K and the aperture was then exposed to a 311-K blackbody at time $t$ = 0. The figure was created using uncorrected RDFs; that is, RFDs computed using the MCRT engine without correction for second-law violations. For the case whose results are shown in Fig. 3, if the model was functioning properly, all temperatures would remain at the initial temperature of 311 K, as mandated by the second law of thermodynamics. However, a disconcerting drift from equilibrium is observed as some temperatures spontaneously rise while others fall. The observed temperature drift with time in this case is unacceptable because of what it portends for problems in which the telescope is exposed to temperatures other than its initial temperature, for problems involving electrical heating, or for a combination of these two situations. An unknown nonphysical component of temperature variation would presumably be present in these cases.

In creating Fig. 3 and when computing the other results presented in this contribution, three-hundred thousand rays were traced from the smallest surface elements, and proportionally more from the larger surface elements. This means that around 800 million rays were traced, with many of them suffering multiple reflections before being absorbed. Even though the workload was shared by five processors working in parallel in a PC running a highly vectorized algorithm, this typically involved an overnight run.

Figure 6, like Fig. 3, shows the calculated temperature drift of all 2636 volume elements after they were initialized at 311 K, and the aperture was then exposed to a 311-K blackbody at time $t$ = 0. However, unlike in Fig. 3, the results shown in Fig. 6 were obtained using RDFs corrected to satisfy both the first and second laws of thermodynamics, as described in Sec. 3. Note that the temperature scales on the vertical axes in these two figures differ enormously. Temperature drifts exceeding a full degree are observed in Fig. 3, which is based on calculations using uncorrected RDFs; while the maximum drift observed in Fig. 6, which is based on corrected RDFs, is of the order of 10⁻⁹ K. Correcting the RFD matrix has virtually eliminated temperature drift.

To this point, we have considered the coupled radiation and conduction problem and found unacceptably large errors when the RDF matrix is not corrected to eliminate second-law violations. This led us to investigate the influence of second-law violations in a purely radiative exchange problem with thermal mass but no conduction. In fact, our analysis tool allows us to investigate this question by simply zeroing out the thermal conductances linking the volume elements, a one-line modification of the code. When we do this while viewing the 311-K blackbody, we obtain the results shown in Figs. 7 and 8.

Figure 7 corresponds to a purely radiative analysis with thermal mass using the uncorrected RDFs. It clearly exhibits a physically unrealistic violation of the second law. Figure 8, on the other hand, obtained using the corrected RDFs, exhibits virtually rock-solid constant temperatures over the same time interval. We may conclude that errors associated with the imprecision of the MCRT method can be significant in such problems, but that these errors can be eliminated using our proposed RDF correction scheme. The effects of the small number of negative on-diagonal distribution factors, i.e., when $Dii<0$⁠; inherent in this correction scheme are clearly negligible. We attribute this to their small representation in the overall distribution factor population and their inherently small magnitude for the problem at hand.

Note that the same number of rays were traced to create both Figs. 7 and 8, but all of the resulting distribution factors were used to create Fig. 7 while only half of them were used to create Fig. 8. Therefore, one is tempted to argue that twice the necessary number of rays were traced to obtain the correct results. The fallacy in this argument is that only after computing the full $n×n$ matrix of distribution factors is it possible to identify the larger, and therefore more accurate, member of each $Dij, Dji$ pair for use with Eq. (2) to replace the smaller, less accurate member.

Figure 9 is a study of the thermal response of the telescope when suddenly shuttered to deep space at absolute zero after first soaking at a uniform temperature of 311 K. The uncorrected RDF results are plotted with dashed lines while the corrected RDF results are plotted with solid lines. The temperature histories in Fig. 9 using both the uncorrected and corrected RDFs appear to be essentially the same, but that is an artifact of the 70-K temperature range applied to the vertical axis. Figure 10, showing the difference between the two calculations, reveals a difference of more than 0.4 K. Since precision of the order of 10 mK is typically desired, this 400 mK discrepancy due to the second-law violation is clearly significant.

Figure 11 shows the temperature distribution after one hour of operation when viewing a 0-K blackbody using the corrected distribution factors. This is the temperature distribution corresponding to the corrected curves in Fig. 9 after one hour has passed. Due to the spreading effect of heat conduction, all the elements on the baffle are at approximately the same temperature. Similarly, the elements on the endcap are also approximately uniform in temperature, but significantly warmer than the baffle. Figure 12 shows the temperature distribution on the baffle section only for the same situation as in Fig. 11. The change in scale reveals the temperature gradients, with the coldest areas near the entrance aperture as expected due to their more direct exposure to the 0-K environment. Although the gradients seem small, with only a 0.16 K difference, they are significant in the current telescope design application.

When modeling radiation heat transfer among thousands of surface elements, violations of the second law of thermodynamics due to imprecisions in estimating RDFs using the Monte Carlo ray-trace (MCRT) method can be significant. Moreover, when uncorrected RDFs are used in a time-marching transient conduction analysis, the impact of these imprecisions can accumulate with time, leading to unacceptable errors. Temperatures are observed to evolve with time in initially isothermal adiabatic enclosures, in stark violation of the second law. The completely successful strategy presented here is fully justified through fundamental examples and then demonstrated for a massively discretized ERB instrument concept consisting of 2636 volume elements and requiring nearly 7 million RDFs. In this new procedure, for any pair of surface elements, the MCRT-based RDF estimate from the surface with the smaller value of $εA$ to the surface with the larger value is used along with the second-law-based reciprocity statement to calculate the RFD from the surface having the larger value of $εA$ to the one having the smaller value. The first law statement for RDFs is then used to compute the on-diagonal members of the RDF matrix. This procedure produces a full RDF matrix completely devoid of violations of both the first and second laws of thermodynamics. A small number of the on-diagonal RFDs turn out to be negative, but their impact on the solution is demonstrated to be negligible.

The authors gratefully acknowledge partial support for the effort reported here by NASA’s Langley Research Center through Science Systems and Applications, Inc. Subcontract No. 21606-16-036. A special debt of gratitude is owed Dr. Anum B. Ashraf at NASA Langley for her technical contributions to this effort.

• National Aeronautics and Space Administration (Funder ID: 10.13039/100000104).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.