Pairing Directional Solar Inputs From Ray Tracing to Solar Receiver/Reactor Heat Transfer Models on Unstructured Meshes: Development and Case Studies

Detailed computational modeling is used to evaluate and optimize the designs of solar receivers/reactors. Simultaneous capture of heat transfer, fluid dynamics, radiative exchange, and chemical reaction phenomena is well documented [1,2] and supported by commercial computational fluid dynamics (CFD) software [3]. Adapting CFD for concentrated solar irradiation is not straightforward due to intense, directional external radiation, which introduces unique and less-documented modeling challenges. Currently, modeling solar irradiation requires manually approximating complex radiative boundary conditions or including concentrating solar infrastructure in the model domain. Either route reduces model accuracy and/or increases model complexity. An alternative, algorithmically simple method is presented for direct mapping from Monte Carlo ray tracing (MCRT) to finite volume approximations to the energy and radiative transport equations (FV-RTE).

Monte Carlo ray tracing and FV-RTE are two methods for incorporating concentrated solar radiation into receiver/reactor models. MCRT is a statistical method to predict radiative absorption by surfaces and participating media with specified radiative properties [4]. Challenges in implementing MCRT for solar receivers/reactors occur when modeling optically thick media, wavelength/temperature-dependent radiative properties, and complex geometries. These conditions are more tractable with FV-RTE methods, including the discrete ordinates (DO) method which is supported in CFD software [5]. FV-RTE methods are computationally expensive for domains requiring high mesh resolutions, like combined models of solar receivers/reactors and collecting/generating infrastructure. Fortunately, MCRT and FV-RTE have complementary strengths. MCRT can model the external radiation to the receiver/reactor, and FV-RTE can capture internal radiative transport within the receiver/reactor. The models are a powerful tool for modeling radiative heat transfer to and within solar receivers/reactors when paired together.

Previous pairings of MCRT and FV-RTE have been implemented by employing (1) non-overlapping or (2) overlapping domains. In non-overlapping schemes [6–9], models of the external irradiation are mapped to the CFD model at a common boundary surface through an intensive process of translating the directional ray-trace to intensity boundary conditions. Non-overlapping schemes are often necessary in designs with falling/entrained solid particles that participate in the radiative exchange [9] and are often functionally energy conservative. However, fine meshes, directionally aligned with the external radiation, are required to mitigate discretization error, which increases computation time and inhibits convergence. Therefore, non-overlapping schemes typically must be restricted to two-dimensional (2D) domains [7].

In overlapping schemes [10–14], portions of the MCRT and CFD modeling domains are spatially coincident. The external radiation in these schemes is modeled from the collector/generator until it is (1) absorbed by internal receiver/reactor surfaces or media or (2) rejected from the receiver/reactor by reflection, transmission, or scattering. Absorbed external radiation is introduced to the CFD model as a boundary or volumetric source. Re-emission and internal irradiation within the receiver/reactor are captured separately, by the CFD model. For diffuse surfaces, overlapping schemes permit coarser CFD meshes that are not required to be directionally aligned with the external intensities. Errors result during CFD discretization from representing curved geometries in the MCRT domain as planar approximations in the CFD domain [14] and are mitigated by mesh refinement.

Absorbed rays in overlapping schemes have been implemented as constant surface [15] or sub-surface [16,17] averaged fluxes. These approximations maintained energy conservation but reduced spatial accuracy between the MCRT and CFD domains. Alternatively, absorbed rays have been binned within a gridded MCRT modeling domain to produce a spatial external irradiation profile and interpolated to the CFD mesh [10]. This approximation maintained spatial accuracy for sufficiently fine grids but did not guarantee energy conservation.

The proposed direct mapping method for MCRT–CFD introduces no spatial or energy conservation errors beyond the MCRT precision and CFD mesh resolution. The method is compatible with 2D and three-dimensional (3D) domains, structured and unstructured meshes, commercial and in-house ray tracing and CFD models, and surfaces of any shape and orientation. Meshing complex geometries with concentrated energy sources often introduces spatial error and/or instability. The method is applied to a full CFD mesh, without the need for specialized meshing procedures for external radiation, and mitigates error in elements with high skew and size variation. The mapping is performed once per mesh and external input, reducing the additional computational expense.

The direct mapping method is demonstrated in ansys fluent via user-defined functions (UDFs). A hybrid method is also proposed which uses the same mapping algorithm paired with a nearest-neighbor filter to limit consideration to nearby elements. The hybrid option improves the computational efficiency of the mapping. Validation is performed for the cavity-type solar thermochemical inclined granular-flow reactor (STInGR) [18,19], as shown in Fig. 1, with direct radiative input from a seven-lamp high-flux solar simulator (HFSS) [20]. Implementation directions for ansys fluent v19.0 and example UDF source codes written in C are available in Appendices A and B of Ref. [21].

A point mapping algorithm is applied to ray intersections to translate absorbed external irradiation from an MCRT model to boundary sources in a CFD model. The mapping for a group of $r⇀k$⁠, depicted as overlaid points, is given in Fig. 2 for a solar reactor aperture and cavity. The points are absorbed on the surfaces of the reactor, assuming non-participating gases in the cavity. A detail view shows the mapping for two $r⇀k$ on a single mesh face F_i of the CFD model. The $r⇀k$ are transformed to a local, barycentric coordinate system (⁠$t^1$⁠,$t^2$⁠) on F_i to determine whether they fall within the boundaries of, and are to be mapped to, F_i. The process is repeated until all the $r⇀k$ are sorted into a F_i.

S_i is defined using triangular geometry. Any F_i with more complex shapes are subdivided into the minimum number required M triangular sub-elements, based on the number of F_i vertices, which are each tested. The direct mapping method is depicted in Fig. 2 with two example intersections at the given F_i, where one intersection falls within and one falls outside F_i.

The first term on the left-hand side is the conductive heat flux in the direction normal to the surface n. The first and second terms on the right-hand side are the convective and radiative heat fluxes, respectively. The convective heat flux uses an overall heat transfer coefficient U to account for contact resistance and/or thin-wall conduction. The net radiative heat flux $qR′′$ resulting from the internal radiative exchange is computed via the FV-RTE.

Monte Carlo ray tracing has been used to model radiative heat transfer in participating media in solar receivers/reactors with quartz windows [22–25], reticulated porous ceramics/metal foams [26–29], and particulates [30–33]. Participating media attenuates directional radiative intensities by absorption and/or scattering within the modeled medium, augmenting the resultant radiation densities and surface heat fluxes. The D and $s^k$ defining a given $r⇀k$ are additionally dependent upon absorption, scattering, and interface refraction.

Volumetrically absorbed rays are mapped into the discretized volumetric domain by a similar algorithm to surface mapping: globally defined points of absorption $r⇀k$ are transformed into a local, barycentric coordinate system (t₁, t₂, t₃) based on the tetrahedral element C_i. Any C_i with more complex shapes are first subdivided into the minimum number required M tetrahedral sub-elements, based on the number of C_i vertices, which are each tested. Subdivision of a given element C_i and volumetric mapping of two $r⇀k$⁠, depicted as overlaid points, are shown in Fig. 3, where the transparent surface bounds a sub-element.

The direct mapping method is depicted in Fig. 3 for two example intersections at the given C_i, where one intersection falls within C_i and one falls outside.

The hybrid computational time t compared with nearest-neighbor time t_n or original barycentric time t_b is a function of the neighborhood size relative to the size of the domain. The t/t_n relationship and energy losses P are shown in Fig. 4 for the hybrid method as a function of d_n/N scaled by the domain characteristic length L_c. White markers below the barycentric t_b/t_n line are values of d_n/N for which hybrid mapping improves t over solely barycentric mapping. As d_n/N → 0, t approaches t_n and the markers approach the nearest-neighbor line at unity. However, P grows quickly due to a small d_n/N excluding the correct F_i or C_i from S_n,i. As d_n/N → L_c, t surpasses t_b and the markers exceed the barycentric line due to the computational cost of the S_n,i calculation and the fact that a large d_n/N filters few, if any, elements. The energy loss decreases to P = 0, though, as the correct F_i or C_i is always in S_n,i. For an optimally sized d_n/N, hybrid mapping greatly reduces t but still maintains P = 0.

While solely nearest-neighbor direct mapping is more computationally efficient than the solely barycentric and hybrid nearest-neighbor/barycentric methods, amplification of discretization error and the introduction of instabilities by erroneously mapping rays to neighboring cells are possible. This scenario is demonstrated in Fig. 3 for the $r⇀k$ within C_i: since d_i+1 < d_i, nearest-neighbor direct mapping maps the $r⇀k$ to the prism-shaped cell despite the $r⇀k$ falling within C_i. Such errors are more likely to occur for meshes with high skew and large, abrupt element size changes, biasing $r⇀k$ toward nearby smaller cells.

Flowcharts of the direct mapping algorithm are shown in Figs. 5(a) and 5(b) for surfaces and volumes.

Validation of energy conservation and spatial preservation by the method was performed for STInGR. STInGR was designed for 5 kW_th scale reduction of redox-active materials directly heated by an HFSS [11]. The combined HFSS/STInGR system, shown in Fig. 6 with the relevant components of each system labeled, was modeled using an overlapping scheme. The 142 mm diameter, 5 mm thick quartz window was modeled as a specularly reflecting, non-scattering, participating medium [25]. The empty cavity and conical frustum were modeled as diffusely reflecting alumina surfaces.

The radiative input was modeled using an MCRT of an HFSS comprised seven Xe arc lamps mounted in truncated ellipsoidal reflectors and aligned to a common focal point [20]. The MCRT was extended to include the semitransparent quartz window and the surfaces comprising the STInGR aperture and internal cavity. The aperture of STInGR was aligned to the HFSS focal plane in the MCRT domain. Emission by the HFSS was assumed to be within the solar spectrum. The MCRT produced 10⁶ rays for each lamp and predicted 8.77 kW_th of radiation absorption by the STInGR surfaces and window. The ray intersections from the MCRT were mapped to a 3D mesh produced in ansys Mesh for CFD models in ansys fluent. The mesh consisted of 59,209 unstructured triangular and quadrilateral face elements and 617 unstructured tetrahedral, hexahedral, or prismatic volumetric elements. The mesh of the STInGR cavity inclined slope was controlled to a uniform structured grid.

The mapped cavity surface-absorbed irradiation from the HFSS was implemented as boundary face sources at the solid-fluid interfaces (Eq. (6)). The mapped window volume-absorbed radiative intensity was modeled as a thermal energy source (Eq. (10)).

Three case studies were performed for 2D surface mapping. The studies investigated the preservation of spatial variation and energy conservation for direct mapping for three STInGR surfaces: the inclined slope, ceiling, and conical frustum, labeled as surfaces 1, 2, and 3 in Fig. 7, respectively. Structured/unstructured meshes and flat/curved surfaces were represented among the studies. The direct mapping method was compared with the profile interpolation process in ansys fluent. The CFD mesh was structured on surface 1 so that profile interpolation and direct mapping produced identical results for the same MCRT grid. The CFD mesh was unstructured for surface 2 so that direct mapping provided improved energy conservation and spatial accuracy compared with profile interpolation. The CFD mesh was unstructured for surface 3, the surface geometry was complex (conical), and the external irradiation gradients were high due to proximity to the HFSS focal point. For this case, grid and mesh resolution deficiencies potentially produced significant energy losses with profile interpolation.

For the profile interpolation process, ansys fluent requires tabular irradiation profiles (i.e., lists of fluxes at specified surface coordinates) which are interpolated to the CFD mesh in a nearest-neighbor manner. In other words, fluent requires area-averaged values rather than the clouds of points produced by MCRT. To obtain the irradiation profiles, the MCRT $r⇀k$ were binned to structured grids in matlab [10]. Wireframes of the grids used to bin the $r⇀k$ are shown in Figs. 8(a)–8(c) for the incline slope, the ceiling, and the conical frustum. The grid resolutions were controlled to produce similar numbers of elements to the CFD surface meshes. No binning was required for the direct mapping method, and the $r⇀k$ were directly input to the CFD mesh using the mapping algorithm.

For complete energy conservation and preservation of spatial variation, P → 0 and summed square of error (SSE) → 0, respectively.

An additional case study was performed for 3D volume mapping. The study investigated absorption within the STInGR quartz window to demonstrate the application of the direct volumetric mapping algorithm for participating media. Profile interpolation for volumetric sources was not supported in ansys fluent v19.0, requiring a different comparison metric for 3D mapping compared with 2D mapping. The case study was instead compared with a nearest-neighbor sorting algorithm as described in Sec. 2.2, as this was the most similar and intuitive corollary to the fluent method, and, therefore, the most likely alternative method.

An isometric view of (1) the solar reactor inclined slope, (2) walls including ceiling, and (3) conical aperture are given in Fig. 7. Each CFD mesh element is shaded according to the absorbed external irradiation, which was applied in ansys fluent as a boundary source via the direct mapping method.

Localized regions of highly concentrated absorbed irradiation from individual lamps were clear along the internal cavity and external front face surfaces (Fig. 7). The localization was particularly evident along the inclined slope and was not capturable with uniform or spatially averaged heat flux profiles. The total power mapped to the reactor surfaces was $∑i(qsun,i′′Ai)$ = 8.29 kW_th. The total power mapped to the quartz window was $∑i(qsun,i′′′Vi)$ = 0.48 kW_th. The total power mapped to the entire CFD mesh in ansys fluent was the same as the 8.77 kW_th originally predicted by the MCRT, confirming energy conservation between the MCRT and CFD.

The first case study is a scenario in which the direct mapping method performs identically to previous methods. Normal views of the STInGR cavity inclined slope are given in Fig. 9, with the mesh elements shaded by the absorbed external irradiation. The slope was discretized as a structured, uniform mesh in the CFD domain with 4250 quadrilateral elements aligned with the MCRT grid. Four local regions of absorbed irradiation from individual HFSS lamps were captured in Figs. 9(a) and 9(b) for the profile interpolation method and the direct mapping method. Peak fluxes up to 300 kW/m² and a total absorbed power of 1.91 kW_th were predicted.

Nearly identical distributions of absorbed external irradiation for the two methods were observed in Fig. 9 due to effectively exact alignment between MCRT grid and CFD mesh. Both methods achieved complete energy conservation (P = 0) and high spatial preservation, with SSE_interp = 0.003 and SSE_map = 0.711, respectively. The spatial errors resulted solely due to differences in numerical precision between C and matlab, as no interpolation between modeling domains was required. Therefore, for structured meshes that align exactly to the binned MCRT grid, the direct mapping method is identical to the interpolated profile method.

The second case study is a scenario in which the direct mapping method performs equivalently to or better than previous methods. Normal views of the meshed STInGR ceiling are given in Figs. 10(a) and 10(b) for the profile interpolation method and the direct mapping method. The ceiling was discretized as an unstructured, non-uniform mesh in the CFD domain with 2864 elements. Three local regions of absorbed irradiation from individual HFSS lamps were captured in Fig. 10. Peak fluxes of 150 kW/m² and a total absorbed power of 0.97 kW_th were predicted.

Slight degradations in spatial accuracy using the profile interpolation method are evident in Fig. 10(a), particularly for F_i with absorbed external irradiation of 50–100 kW/m². The elliptical profile was better preserved by direct mapping, as observed in Fig. 10(b). Spatial accuracy and energy conservation were achieved, respectively, to SSE_map = 5.8 × 10³ < SSE_interp = 7.5 × 10³ and P_map = 0 < P_interp = 0.004. While both methods approximated the spatial profile shape without significant energy losses, the direct mapping method achieved better spatial accuracy (lower SSE) and complete energy conservation.

The third case study is a scenario in which the direct mapping method preserves spatial accuracy better than previous methods and is critical to prevent significant energy losses. A normal view of the STInGR aperture, a conical frustum shape, is given in Figs. 11(a) and 11(b) for the profile interpolation method and the direct mapping method. An inset in the bottom right of Fig. 11(a) shows the spatial variations. The frustum was discretized as an unstructured, non-uniform mesh in the CFD modeling domain with 1182 elements. Radially uniform profiles are present in Fig. 11 from a 2.3 kW_th spillage of concentrated radiation from the HFSS around the aperture.

The direct mapping method qualitatively preserved the spatial profile shape better than profile interpolation. However, the approximation of the conical surface as planar faces in the CFD mesh produced a mismatch between the locations of the MCRT grid and CFD mesh elements and between the total surface areas in the MCRT grid (0.022 m²) and CFD mesh (0.020 m²). This mismatch prevented a meaningful quantitative SSE comparison before and after coupling.

A significantly smaller magnitude of absorbed external irradiation was captured in Fig. 11(a) compared with Fig. 11(b) due to the highly concentrated irradiation in the focal plane. Large flux gradients and limited grid resolution near the aperture led to underestimation during interpolation, producing a peak $qsun′′$ = 15 ≪ 800 kW/m², as shown in the inset of Fig. 11(a). Energy conservation analysis resulted in P_map = 0 < P_interp = 0.98, indicating that profile interpolation introduced large errors in energy conservation—in this case, as high as 98%—depending on (1) the irradiation gradient and (2) MCRT grid/CFD mesh resolutions. Direct mapping, however, was robust even for sharp irradiation profiles and/or coarse meshes.

A comparison of mapping methods for the three studies pictured in Figs. 9–11 revealed the inherent energy conservative nature of the direct mapping method, with spatial accuracy dependent upon the discretization accuracy of the modeled geometry. The method was shown to be independent of the non-trivial process of matching gridded MCRT and meshed CFD modeling domains. Direct mapping achieved equivalent accuracy to the interpolated profile method for aligned MCRT grids/CFD meshes and improved accuracy for misaligned MCRT grids/CFD meshes.

A final case study was performed to demonstrate the barycentric direct mapping method for participating media and to show improved performance over the nearest-neighbor method. A view of the STInGR window depicting the unstructured CFD mesh cell centroids is given in Figs. 12(a) and 12(b) for nearest-neighbor direct mapping and barycentric direct mapping. The volumetrically absorbed radiative intensity was applied as volumetric heat sources in ansys fluent.

Using both algorithms, the profile of the seven-lamp HFSS, roughly symmetric about (0,0), was visible. In Fig. 12(a), however, localized C_i of high or low $qsun′′′$ not present in Fig. 12(b) were evident. The differences were the result of ray misappropriation by the nearest-neighbor algorithm for neighboring cells with significant differences in volume. The most prominent example of ray misappropriation in Fig. 12(a) occurred at the cell centroid near (0.01,0), producing a local hotspot of $qsun′′′≈$ 1.9 × 10⁴ kW/m². The meshed element and centroid are pictured in greater detail in Fig. 3 as element C_i+1. A smoother $qsun′′′$ profile resulted for barycentric mapping, indicative of improved mapping between the MCRT and CFD modeling domains.

Based upon the methodology and observations from the case studies, the advantages and disadvantages of the method are summarized in Secs. 3.5 and 3.6. Note that all disadvantages are also true for other overlapping modeling domain schemes. The disadvantages may also be reasonably mitigated in most scenarios.

1. The method is energy conservative between the MCRT and CFD modeling domains.

2. The method is spatially accurate to within the MCRT and CFD discretization accuracies.

3. The method is compatible with structured and unstructured meshes of arbitrary polygonal or polyhedral construction, for two and three dimensions.

4. The method uses an algorithm that is programmatically simple and is applied using an external code or directly within ansys fluent via UDFs.

5. The method requires just one run of the MCRT and mapping per ray-trace and mesh. It is, therefore, unlikely to be a significant computational expense relative to the combined momentum, energy, and radiation models for most scenarios.

1. Transient mapping is cumbersome using the method. Such cases occur for overlapping schemes with participating media in the band(s) of the solar radiation. Non-overlapping domains or a single computational domain is more appropriate for media with highly temperature-dependent absorption, transmission, reflection, or scattering.

2. Systematic errors in the absorbed heat flux distribution are introduced by approximating curved geometries from the MCRT model with polygonal elements in the CFD model.

3. The method utilizes boundary sources for surface elements. Boundary sources in ansys fluent require defining a boundary region of nonzero thickness, which introduces additional conductive resistance at the interface [3]. The additional resistance becomes negligible for sufficiently thin boundaries with sufficiently high thermal conductivities.

A method for mapping absorbed irradiation and radiative intensities from Monte Carlo ray tracing simulations to computational heat transfer models was presented. The method is preferable to previous work, due to rapid implementation while maintaining energy conservation and spatial patterns between the two modeling domains within mesh precision.

The direct mapping method was implemented for a model of a windowed, solar thermochemical reactor with irradiation from a high-flux solar simulator. The method captured local hotspots from individual lamps without losses in net energy absorption due to direct mapping from the Monte Carlo ray tracing to the computational fluid dynamics mesh. The method was demonstrated to preserve spatial variation and maintain equal or better energy conservation compared with previous models, for various complex geometries and mesh types. The direct mapping method is a useful tool for developing solar infrastructure, as it enables the accurate modeling of highly directional/spatial concentrated external intensity and irradiation.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1650044. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. This material is also based upon work supported by the U.S. Department of Energy SunShot Initiative under Award No. DE-FOA-0000805-1541. This publication derives from work while all authors were affiliated with the Georgia Tech Solar Fuels and Technology Lab (Peter G. Loutzenhiser). At publication time, H. Evan Bush is a Postdoctoral Researcher at Sandia National Labs in Concentrating Solar Technologies (hebush@sandia.gov). Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. At publication time, Andrew Schrader is an Assistant Professor at the University of Dayton in Mechanical and Aerospace Engineering (aschrader1@udayton.edu).

There are no conflicts of interest.

k =

thermal conductivity

n =

integer number

q =

heat rate

r =

position vector

t =

thickness, barycentric coordinate, computational time

A =

area

E =

thermal energy

H =

irradiation

M =

integer number

P =

energy loss/error, MCRT grid to CFD mesh

Q =

thermal energy

S =

set

T =

temperature

V =

volume, integer number of face/cell vertices

U =

heat transfer coefficient

X =

mesh face/cell vertex location vector

V =

velocity vector

L_c =

characteristic length

N_rays =

total number of rays in ray tracing model

a, b, c =

barycentric coordinate locations

α =

total, hemispherical surface absorptance

κ =

linear absorption coefficient

ρ =

density

∞ =

environment conditions

1, 2, 3 =

triangle/tetrahedron edges and coordinates

i =

face index

k =

ray index

local =

local barycentric coordinate system

n =

normal direction, nearest neighbor

N =

nearest neighbor

nb =

neighboring fluid cell

R =

radiation

s =

surface

sun =

real or simulated concentrated sunlight

″ =

value per unit area, e.g., boundary flux

″′ =

value per unit volume, e.g., volumetric source

^ =

unit vector

. =

time rate of change