Abstract

A novel method for pairing surface irradiation and volumetric absorption from Monte Carlo ray tracing to computational heat transfer models is presented. The method is well-suited to directionally and spatially complex concentrated radiative inputs (e.g., solar receivers and reactors). The method employs a generalized algorithm for directly mapping absorbed rays from a Monte Carlo ray tracing model to boundary or volumetric source terms in the computational mesh. The algorithm is compatible with unstructured, two and three-dimensional meshes with varying element shapes. Four case studies were performed on a directly irradiated, windowed solar thermochemical reactor model to validate the method. The method was shown to conserve energy and preserve spatial variation when mapping rays from a Monte Carlo ray tracing model to a computational heat transfer model in ansys fluent.

1 Introduction

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).

Previous pairings of MCRT and FV-RTE have been implemented by employing (1) non-overlapping or (2) overlapping domains. In non-overlapping schemes [69], 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 [1014], 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].

Fig. 1
Fig. 1

2 Methods

Monte Carlo ray tracing allows the directional and spectral modeling of radiative heat transfer by partitioning a radiation source such as concentrated solar collectors (e.g., heliostats and parabolic troughs/dishes) or generators (e.g., HFSS), into rays or individual packets of energy, represented as
$E˙k=Q˙sunNrays$
(1)
where $Q˙sun$ is the radiative power, here evenly partitioned among N rays. Nrays is typically between 105 and 107, with convergence expected for Nrays →∞. The path of a given ray k = 1, 2, …, Nrays is defined via Ref. [4], given as
$r⇀k=r⇀0,k+Ds^k$
(2)
where $r⇀k$ is the ray intersection, located a distance D; direction $s^k$ from the ray origin $r⇀0,k$ (i.e., the location of ray generation or introduction to the model) or previous ray interaction (e.g., reflection and scattering). Ray origin and interaction characteristics are defined stochastically from the domain radiative properties.

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 Fi of the CFD model. The $r⇀k$ are transformed to a local, barycentric coordinate system ($t^1$,$t^2$) on Fi to determine whether they fall within the boundaries of, and are to be mapped to, Fi. The process is repeated until all the $r⇀k$ are sorted into a Fi.

Fig. 2
Fig. 2
The local coordinate system origin for a triangular Fi is defined at any of the three Fi vertices, where the basis vectors t1 and t2 are in the directions from the origin to each of the two remaining Fi vertices, respectively. The two-dimensional, local intersection position vectors are defined as
$r⇀local=at1+bt2$
(3)
where (a,b) are the local coordinate dimensions. The subset of rays that fall within Fi is defined as
$Si={k=1,2,…,Nrays|ai,k,bi,k>0,ai,k+bi,k≤1}$
(4)

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

Computational fluid dynamics model boundary sources for each Fi are computed from mapped, summed ray powers as
$qsun,i′′=(αH0)sun,i=1As,i∑k∈SiE˙k$
(5)
where (αH0)sun is the absorbed component of the external irradiation profile H0,sun and As is the face surface area. The summation is limited only to the rays within the subset S. The boundary sources are incorporated into a general surface-fluid mixed boundary condition to the heat diffusion equation as
$−kn,i∂T∂n|s,i+qsun,i′′=Ui(Ts,i−T∞,i)+qR,i′′$
(6)

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.

2.2 Participating Volumetric Cell Mapping.

Monte Carlo ray tracing has been used to model radiative heat transfer in participating media in solar receivers/reactors with quartz windows [2225], reticulated porous ceramics/metal foams [2629], and particulates [3033]. 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 (t1, t2, t3) based on the tetrahedral element Ci. Any Ci with more complex shapes are first subdivided into the minimum number required M tetrahedral sub-elements, based on the number of Ci vertices, which are each tested. Subdivision of a given element Ci and volumetric mapping of two $r⇀k$, depicted as overlaid points, are shown in Fig. 3, where the transparent surface bounds a sub-element.

Fig. 3
Fig. 3
The three-dimensional, local absorption point position vectors are defined as
$r⇀local=at1+bt2+ct3$
(7)
where (a,b,c) are the local coordinates. The subset of rays that fall within Ci is defined as
$Si={k=1,2,…,N|ai,k,bi,k,ci,k>0,ai,k+bi,k+ci,k≤1}$
(8)

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

Volumetric sources for each Ci are computed from mapped, summed ray powers as
$q˙sun,i′′′=(κ∫4πI(s′)dΩ′)sun,i=1Vi∑k∈SiE˙k$
(9)
where $κ∫4πI(s⇀′)dΩ′$ is the absorbed component of the external intensity from all directions $s⇀′$ and over all solid angles Ω′ for a linear absorption coefficient κ and V is the cell volume. The summation is limited only to the rays within the subset S. The volumetric heat sources are applied at each Ci as a term in the thermal energy transport equation, commonly represented as
$∂Ei∂t+∇⋅(ρVh)nb,i=∇⋅(k∇Ti)+qsun,i′′′+qR,i′′′$
(10)
where the left-hand side includes the transient and advected energy from neighboring fluid cells nb and the right-hand side captures heat diffusion and the energy source terms: $qsun′′′$ and the internal radiative balance computed via the FV-RTE $qR′′′$.

2.3 Hybrid Nearest-Neighbor/Barycentric Mapping.

A hybrid nearest-neighbor/barycentric direct mapping method significantly reduces computation time without altering the direct mapping results. In this method, a neighborhood of the n nearest surface/volumetric elements to each $r⇀k$ is identified and the original barycentric direct mapping algorithm is applied exclusively to the neighborhood. The nearest-neighbor algorithm used to identify the neighborhood, or subset of elements eligible for mapping, is defined as
$Sn,i={k=1,2,…,Nrays|di,k≤dn/N,k}$
(11)
where dn/N,k is the size of the neighborhood of n/Nth percentile nearest elements, defined from distances $di,k=||r⇀k−r⇀Ci||21/2$ calculated between the $r⇀k$ and cell centroids $r⇀Ci$. The barycentric algorithm is applied only to Sn,i, reducing the number of elements over which the algorithm iterates. The results are identical to the solely barycentric mapping for a Sn,i of the entire domain, and identical to solely nearest-neighbor mapping for a Sn,i of a single cell.

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

Fig. 4
Fig. 4

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 Ci: since di+1 < di, nearest-neighbor direct mapping maps the $r⇀k$ to the prism-shaped cell despite the $r⇀k$ falling within Ci. 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.

Fig. 5
Fig. 5

2.4 Computational Models.

Validation of energy conservation and spatial preservation by the method was performed for STInGR. STInGR was designed for 5 kWth 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.

Fig. 6
Fig. 6

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 106 rays for each lamp and predicted 8.77 kWth 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.

Fig. 7
Fig. 7

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.

Fig. 8
Fig. 8
To quantitatively compare the energy conservation of both methods, the total energy loss from the MCRT grid to the CFD mesh sources was calculated as
$P=|(∑iAiqsun,i′′)CFD−(∑jAjqsun,j′′)MCRT|(∑jAjqsun,j′′)MCRT$
(12)
where i and j are the indices of CFD mesh and MCRT grid elements, respectively.
To quantitatively compare the spatial accuracy of the methods, the CFD mesh boundary source fluxes from (1) the profile interpolation process in ansys fluent and (2) the direct mapping method were both linearly interpolated back to the structured surface grids used in the profile binning. These were quantitatively compared with the fluxes from the originally binned $r⇀k$ via the sum of square errors, calculated as
$SSE=∑j(qsun,j′′|CFD−qsun,j′′|MCRT)2$
(13)

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.

3 Results and Discussion

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 kWth. The total power mapped to the quartz window was $∑i(qsun,i′′′Vi)$ = 0.48 kWth. The total power mapped to the entire CFD mesh in ansys fluent was the same as the 8.77 kWth originally predicted by the MCRT, confirming energy conservation between the MCRT and CFD.

3.1 Surface Case Study I: Inclined Slope.

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/m2 and a total absorbed power of 1.91 kWth were predicted.

Fig. 9
Fig. 9

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 SSEinterp = 0.003 and SSEmap = 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.

3.2 Surface Case Study II: Reactor Ceiling.

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/m2 and a total absorbed power of 0.97 kWth were predicted.

Fig. 10
Fig. 10

Slight degradations in spatial accuracy using the profile interpolation method are evident in Fig. 10(a), particularly for Fi with absorbed external irradiation of 50–100 kW/m2. The elliptical profile was better preserved by direct mapping, as observed in Fig. 10(b). Spatial accuracy and energy conservation were achieved, respectively, to SSEmap = 5.8 × 103 < SSEinterp = 7.5 × 103 and Pmap = 0 < Pinterp = 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.

3.3 Surface Case Study III: Reactor Aperture.

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 kWth spillage of concentrated radiation from the HFSS around the aperture.

Fig. 11
Fig. 11

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 m2) and CFD mesh (0.020 m2). 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/m2, as shown in the inset of Fig. 11(a). Energy conservation analysis resulted in Pmap = 0 < Pinterp = 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. 911 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.

3.4 Volumetric Case Study: Quartz Window.

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.

Fig. 12
Fig. 12

Using both algorithms, the profile of the seven-lamp HFSS, roughly symmetric about (0,0), was visible. In Fig. 12(a), however, localized Ci 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 × 104 kW/m2. The meshed element and centroid are pictured in greater detail in Fig. 3 as element Ci+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.

4 Conclusion

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.

Acknowledgment

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).

Conflict of Interest

There are no conflicts of interest.

Nomenclature

• k =

thermal conductivity

•
• n =

integer number

•
• q =

heat rate

•
• r =

position vector

•
• t =

thickness, barycentric coordinate, computational time

•
• A =

area

•
• E =

thermal energy

•
• H =

•
• 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

•
• Lc =

characteristic length

•
• Nrays =

total number of rays in ray tracing model

•
• a, b, c =

barycentric coordinate locations

Greek Symbols

• α =

total, hemispherical surface absorptance

•
• κ =

linear absorption coefficient

•
• ρ =

density

Subscripts

• ∞ =

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 =

•
• s =

surface

•
• sun =

real or simulated concentrated sunlight

Superscripts/Accents

• ″ =

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

•
• ″′ =

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

•
• ^ =

unit vector

•
• . =

time rate of change

References

1.
Patankar
,
S.
,
1980
,
Numerical Heat Transfer and Fluid Flow
,
CRC Press
,
Washington, DC
.
2.
Anderson
,
D.
,
Tannehill
,
J.
, and
Pletcher
,
R.
,
Computational Fluid Mechanics and Heat Transfer, Series in Computational Methods in Mechanics and Thermal Sciences
, 2nd ed.,
Taylor and Francis
,
London
.
3.
2013
, “ANSYS Fluent,” ANSYS, Inc., Canonsburg, PA.
4.
Modest
,
M. F.
,
2013
,
,
,
Oxford, UK
.
5.
ANSYS
,
2017
, “ANSYS Fluent,” ANSYS, Inc., Canonsburg, PA.
6.
Mecit
,
A. M.
, and
Miller
,
F.
,
2014
A Comparison Between the Monte Carlo Ray Trace and the FLUENT Discrete Ordinates Methods for Treating Solar Input to a Particle Receiver
,”
ASME 2014 8th International Conference on Energy Sustainability Collocated With the ASME 2014 12th International Conference on Fuel Cell Science, Engineering and Technology, Volume 1: Combined Energy Cycles, CHP, CCHP, and Smart Grids; Concentrating Solar Power, Solar Thermochemistry and Thermal Energy Storage; Geothermal, Ocean, and Emerging Energy Technologies; Hydrogen Energy Technologies; Low/Zero Emission Power Plants and Carbon Sequestration; Photovoltaics; Wind Energy Systems and Technologies
,
Boston, MA
,
June 30–July 2
, p.
V001T002A013
.
7.
Craig
,
K. J.
,
Moghimi
,
M. A.
,
Rungasamy
,
J.
,
Marsberg
,
J.
, and
Meyer
,
P.
,
2016
, “
Finite-Volume Ray Tracing Using Computational Fluid Dynamics in Linear Focus CSP Applications
,”
Appl. Energy
,
183
, pp.
241
256
. 10.1016/j.apenergy.2016.08.154
8.
Wang
,
F.
,
Tan
,
J.
,
Yong
,
S.
,
Tan
,
H.
, and
Chu
,
S.
,
2014
, “
Thermal Performance Analyses of Porous Media Solar Receiver With Different Irradiative Transfer Models
,”
Int. J. Heat Mass Transfer
,
78
, pp.
7
16
. 10.1016/j.ijheatmasstransfer.2014.06.035
9.
Khalsa
,
S. S. S.
, and
Ho
,
C. K.
,
2011
, “
Radiation Boundary Conditions for Computational Fluid Dynamics Models of High-Temperature Cavity Receivers
,”
ASME J. Sol. Energy Eng.
,
133
(
3
), p.
031020
. 10.1016/10.1115/1.4004274
10.
Bush
,
H. E.
,
Schlichting
,
K.-P.
,
Gill
,
R. J.
,
Jeter
,
S. M.
, and
Loutzenhiser
,
P. G.
,
2017
, “
Design and Characterization of a Novel Upward Flow Reactor for the Study of High-Temperature Thermal Reduction for Solar-Driven Processes
,”
ASME J. Sol. Energy Eng.
,
139
(
5
), p.
051004
. 10.1115/1.4037191
11.
,
A. J.
,
De Dominicis
,
G.
,
Schieber
,
G. L.
, and
Loutzenhiser
,
P. G.
,
2017
, “
Solar Electricity via an Air Brayton Cycle With an Integrated Two-Step Thermochemical Cycle for Heat Storage Based on Co3O4/CoO Redox Reactions III: Solar Thermochemical Reactor Design and Modeling
,”
Sol. Energy
,
150
, pp.
584
595
. 10.1016/j.solener.2017.05.003
12.
Qiu
,
Y.
,
He
,
Y.-L.
,
Wu
,
M.
, and
Zheng
,
Z.-J.
,
2016
, “
A Comprehensive Model for Optical and Thermal Characterization of a Linear Fresnel Solar Reflector With a Trapezoidal Cavity Receiver
,”
Renew. Energy
,
97
, pp.
129
144
. 10.1016/j.renene.2016.05.065
13.
Qiu
,
Y.
,
He
,
Y.-L.
,
Cheng
,
Z.-D.
, and
Wang
,
K.
,
2015
, “
Study on Optical and Thermal Performance of a Linear Fresnel Solar Reflector Using Molten Salt as HTF With MCRT and FVM Methods
,”
Appl. Energy
,
146
, pp.
162
173
. 10.1016/j.apenergy.2015.01.135
14.
He
,
Y.-L.
,
Xiao
,
J.
,
Cheng
,
Z.-D.
, and
Tao
,
Y.-B.
,
2010
, “
A MCRT and FVM Coupled Simulation Method for Energy Conversion Process in Parabolic Trough Solar Collector
,”
Renew. Energy
,
36
(
3
), pp.
976
985
. 10.1016/j.renene.2010.07.017
15.
Lapp
,
J.
,
Davidson
,
J. H.
, and
Lipiński
,
W.
,
2013
, “
Heat Transfer Analysis of a Solid-Solid Heat Recuperation System for Solar-Driven Nonstoichiometric Redox Cycles
,”
ASME J. Sol. Energy Eng.
,
135
(
3
), p.
031004
. 10.1115/1.4023357
16.
Koepf
,
E.
,
,
S. G.
,
Steinfeld
,
A.
, and
,
A. K.
,
2012
, “
A Novel Beam-Down, Gravity-Fed, Solar Thermochemical Receiver/Reactor for Direct Solid Particle Decomposition: Design, Modeling, and Experimentation
,”
Int. J. Hydrogen Energy
,
37
(
22
), pp.
16871
16887
. 10.1016/j.ijhydene.2012.08.086
17.
Alonso
,
E.
, and
Romero
,
M.
,
2015
, “
A Directly Irradiated Solar Reactor for Kinetic Analysis of Non-Volatile Metal Oxides Reductions
,”
Int. J. Hydrogen Energy
,
39
(
9
), pp.
1217
1228
. 10.1002/er.3320
18.
,
A. J.
,
Bush
,
H. E.
,
Ranjan
,
D.
, and
Loutzenhiser
,
P. G.
,
2020
, “
Aluminum-Doped Calcium Manganite Particles for Solar Thermochemical Energy Storage: Reactor Design, Particle Characterization, and Heat and Mass Transfer Modeling
,”
Int. J. Heat Mass Transfer
,
152
, p.
119461
. 10.1016/j.ijheatmasstransfer.2020.119461
19.
,
A. J.
,
Schieber
,
G. L.
,
Ambrosini
,
A.
, and
Loutzenhiser
,
P. G.
,
2020
, “
Experimental Demonstration of a 5 KWth Granular-Flow Reactor for Solar Thermochemical Energy Storage With Aluminum-Doped Calcium Manganite Particles
,”
Appl. Therm. Eng.
,
173
, p.
115257
. 10.1016/j.applthermaleng.2020.115257
20.
Gill
,
R.
,
Bush
,
H. E.
,
Haueter
,
P.
, and
Loutzenhiser
,
P.
,
2015
, “
Characterization of a 6kW High-Flux Solar Simulator With an Array of Xenon Arc Lamps Capable of Concentrations of Nearly 5000 Suns
,”
Rev. Sci. Instrum.
,
86
(
12
), p.
8
. 10.1063/1.4936976
21.
Bush
,
H. E.
,
2019
, “
Development and Characterization of Novel Reduction-Oxidation Active Materials for Two-Step Solar Thermochemical Cycles
,”
Ph.D. Dissertation
,
Georgia Institute of Technology
,
Atlanta, GA
.
22.
Dai
,
G.-L.
,
Xia
,
X.-L.
, and
Hou
,
G.-F.
,
2014
, “
Transmission Performances of Solar Windows Exposed to Concentrated Sunlight
,”
Sol. Energy
,
103
, pp.
125
133
. 10.1016/j.solener.2014.01.036
23.
Yong
,
S.
,
Fu-Qiang
,
W.
,
Xin-Lin
,
X.
,
He-Ping
,
T.
, and
Ying-Chun
,
L.
,
2011
, “
,”
Int. J. Hydrogen Energy
,
36
(
19
), pp.
12148
12158
. 10.1016/j.ijhydene.2011.07.013
24.
Mecit
,
A. M.
, and
Miller
,
F.
, “
Optical Analysis of a Window for Solar Receivers Using the Monte Carlo Ray Trace Method
,”
Proceedings of the ASME 2013 7th International Conference on Energy Sustainability Collocated With the ASME 2013 Heat Transfer Summer Conference and the ASME 2013 11th International Conference on Fuel Cell Science, Engineering and Technology, American Society of Mechanical Engineers
,
Minneapolis, MN
,
July 14–19
, p.
V001T011A003
.
25.
Müller
,
R.
, and
Steinfeld
,
A.
,
2007
, “
Band-Approximated Radiative Heat Transfer Analysis of a Solar Chemical Reactor for the Thermal Dissociation of Zinc Oxide
,”
Sol. Energy
,
81
(
10
), pp.
1285
1294
. 10.1016/j.solener.2006.12.006
26.
Cunsolo
,
S.
,
Oliviero
,
M.
,
Harris
,
W. M.
,
Andreozzi
,
A.
,
Bianco
,
N.
,
Chiu
,
W. K. S.
, and
Naso
,
V.
,
2015
, “
Monte Carlo Determination of Radiative Properties of Metal Foams: Comparison Between Idealized and Real Cell Structures
,”
Int. J. Therm. Sci.
,
87
, pp.
94
102
. 10.1016/j.ijthermalsci.2014.08.006
27.
Wei
,
G.
,
Huang
,
P.
,
Xu
,
C.
,
Chen
,
L.
,
Ju
,
X.
, and
Du
,
X.
,
2017
, “
Experimental Study on the Radiative Properties of Open-Cell Porous Ceramics
,”
Sol. Energy
,
149
, pp.
13
19
. 10.1016/j.solener.2017.04.002
28.
Cui
,
F. Q.
,
He
,
Y. L.
,
Cheng
,
Z. D.
,
Li
,
D.
, and
Tao
,
Y. B.
,
2012
, “
Numerical Simulations of the Solar Transmission Process for a Pressurized Volumetric Receiver
,”
Energy
,
46
(
1
), pp.
618
628
. 10.1016/j.energy.2012.07.044
29.
Li
,
Y.
,
Xia
,
X.-L.
,
Ai
,
Q.
,
Sun
,
C.
, and
Tan
,
H.-P.
,
2018
, “
Pore-Level Determination of Spectral Reflection Behaviors of High-Porosity Metal Foam Sheets
,”
Infrared Phys. Technol.
,
89
, pp.
77
87
. 10.1016/j.infrared.2017.12.016
30.
Akolkar
,
A.
, and
Petrasch
,
J.
,
2011
, “
Tomography Based Pore-Level Optimization of Radiative Transfer in Porous Media
,”
Int. J. Heat Mass Transfer
,
54
(
23
), pp.
4775
4783
. 10.1016/j.ijheatmasstransfer.2011.06.017
31.
Marti
,
J.
,
Roesle
,
M.
, and
Steinfeld
,
A.
,
2014
, “
Combined Experimental-Numerical Approach to Determine Radiation Properties of Particle Suspensions
,”
ASME J. Heat Transfer
,
136
(
9
), p.
092701
. 10.1115/1.4027768
32.
Marti
,
J.
,
Roesle
,
M.
, and
Steinfeld
,
A.
,
2014
, “
Experimental Determination of the Radiative Properties of Particle Suspensions for High-Temperature Solar Receiver Applications
,”
Heat Transfer Eng.
,
35
(
3
), pp.
272
280
. 10.1080/01457632.2013.825173
33.
Tien
,
C. L.
,
1988
, “
Thermal Radiation in Packed and Fluidized Beds
,”
ASME J. Heat Transfer
,
110
(
4b
), pp.
1230
1242
. 10.1115/1.3250623