<figure><img src="image_1.png" /></figure>

HAL
open science

# Coupled surface-to-surface radiative and conductive heat transfer with an immersed volume method

Rémi Gerard, Elie Hachem, Giulia Lissoni, Aurélien Larcher

<figure><img src="image_2.png" /></figure>

## To cite this version:

Rémi Gerard, Elie Hachem, Giulia Lissoni, Aurélien Larcher. Coupled surface-to-surface radiative and conductive heat transfer with an immersed volume method. 2024. hal-04769425

HAL Id: hal-04769425

https://hal.science/hal-04769425v1

Preprint submitted on 6 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coupled surface-to-surface radiative and conductive heat transfer with an immersed volume method

R. Gerardª, E. Hachemª, G. Lissoniᵇ, A. Larcherª,

ª MINES ParisTech, PSL - Research University, CEMEF - Centre de mise en forme des matériaux, CNRS UMR 7635, CS 10207 rue Claude Daunesse 06904 Sophia Antipolis Cedex, France

ᵇSciences Computers Consultants, 10 rue du plateau des Glières, 42000 Saint-Etienne, France

## Abstract

This study presents a novel approach of surface-to-surface (S2S) radiation in the context of an immersed volume method. Radiating facets are reconstructed from the solid-liquid interface. Obstructed view factors are accurately computed while energy balance is ensured through flux correction. The resulting radiative flux is then coupled with a volumic thermal tranport equation as a volumetric source term. This paper presents pure radiation and radiation-conduction simulations which agree with analytical results.

**Keywords:** Surface radiation, View factors, LevelSet, Obstruction, Conduction

## 1. INTRODUCTION

Radiation is an important part of thermal transfer, and one that plays a dominant role in a number of industrial applications that involve locally or globally high temperatures, such as in furnace vitrification [1, 2, 3], charcoal combustion chambers [4, 5, 6], steel ingot casting [7, 8, 9], thermal shielding [10, 11, 12], additive manufacturing [13, 14, 15] or photovoltaic central light collection [16, 17, 18]. Due to its strong reliance on directionality and its important local variability, the radiative transfer equation (RTE) that describes radiation at any given point in the system, is a complex integro-differential equation for which the direct numerical resolution is very costly. This has urged the scientific community to develop a number of radiation models over the years. While some of these methods, such as the P1 model [19, 20, 4] or the Rosseland approach [10, 21, 22] average out the contribution of radiation as a boundary condition, some others, like the DOM [23, 24, 25], the DTRM [25, 26, 27] or the Monte-Carlo ray tracing [28, 29, 30] instead account for the various trajectories of the heating rays. A powerful and adaptable modelling, the Surface-to-Surface (S2S) or radiosity method, has been recently gaining a lot of traction, providing a good compromise between computational cost and accuracy [31, 32, 33, 34, 35, 36, 37]. This method allows to account for the directionality of radiation and strong local variations when the radiating surfaces are diffusive and approximated as grey bodies.

The radiosity method has been recorded in thermal transfer handbooks since 1965 [38]. It has then been picked up in the 80s by computer graphic scientists, in order to describe how light radiates from surfaces to surfaces for correctly rendering illuminated images [39, 40, 41]. Since light is a form of radiation, the method is indeed useful to both domains of research.

The S2S method is applied in a domain where all the surfaces involved in radiation are split into a discrete set of *n* elementary small surfaces. The core principle of the S2S method then relies on the following fundamental operation: evaluating how the radiation leaving each of the *n* surfaces is distributed to hit all the other surfaces. This is called the view factor matrix computation. This process can be adjusted to the problem according to a number of specifications: the *n* surfaces taken for radiation computation can be

Email addresses: remi.gerard@mines-paristech.fr (R. Gerard), elie.Hachem@mines-paristech.fr (E. Hachem), glissoni@scconsultants.com (G. Lissoni), aurelien.larcher@mines-paristech.fr (A. Larcher)
directly extracted from the surface elements of the mesh. But this set can be refined if a greater accuracy is required; or facets can also be gathered into clusters if, on the contrary, speed of execution is deemed more important than accuracy [39]. These choices can be made at a local level, and indeed some parts of the surface mesh can be refined where most of the radiation occurs while some others can be clustered in areas where radiation plays a less dominant role. Furthermore, this operation is only dependent on the geometry of the system: it can be precomputed, and as long as the geometry does not change during simulation, the computed values can be used as-is for all timesteps. This adaptability and simplicity constitutes the most interesting point of the radiosity method.

Computing view factors is key to the speed and accuracy of the S2S method, and many S2S variants differ by the way they process that step. Hottel and Sarofim [42] developed the so-called "crossed-string method" which allows for an analytical computation of 2D view factors. Walton [43] transformed the double suface integral computation of view factors into a single linear integral, including the complicated case of partially obstructed view factors. Cohen and Greenberg [44] devised a method called the hemi-cube technique, in which they installed a square grid surrounding a given facet and then projected each other facet onto it. This both allowed for an efficient computation of obstruction and turned the view factor computation from that between two arbitrarily-shaped facets into that between one arbitrarily-shaped facet and one square surface parallel to the facet, simplifying the computation. Franck et al. [45] or Woop et al [46] have instead used a Monte-Carlo method, massively launching rays in random directions from each radiating surface, and the fraction of rays hitting the opposite surface out of the total number of launched rays gives a good approximation of the view factors between them.

Though it can accurately handle very detailed geometries, or radiation between moving or deforming items, the S2S method tends to grow rather cumbersome in these cases. Modelling radiation between moving or deforming items, for which the mesh varies over time, is very costly since the view factor matrix has to be recomputed each time the mesh is modified. Besides, in order to accurately model precise geometries, a very refined mesh for the radiating items needs to be provided, which makes the process more computationally costly. A more general approach would be to compute S2S radiation within the framework of an immersed mesh. Indeed, since immersed meshes do not have to be recomputed when items in the system move, we can save the computation time otherwise required to recompute the mesh in case of body-fitted approaches. We do have to recompute the view factors still, but this stands for a good efficiency improvement nonetheless. Furthermore, the immersed method only requires a moderately refined mesh for the immersed items, and still allows for an accurate S2S modelling through the reconstruction of virtual radiating facets at the interface. This allows us to lighten the load of the simulation data even more. This papers investigates the matter of radiation-conduction problems with an increased focus on obstruction handling. The main contribution of this paper is thus the devising of S2S model adapted to the use within an immersed mesh context that can prove accurate even with moderately rough meshes and be used as a solid foundation for later fast modellings of radiating moving items.

The paper is composed of the following sections. Section 1 will review the numerical framework used, and recall the main points of the immersed mesh approach. Section 2 will describe the S2S model used, its working assumptions, and its numerical implementation. Section 3 will detail various test cases where our model is compared to analytically-expected results. A brief conclusion will sum up the main contributions of this paper.

## 2. NUMERICAL FRAMEWORK

### 2.1. *Interface description*

The immersed approaches, also called monolithic approaches, impose the use of an appropriate constitutive equation describing both the fluid and the solid domains. Their use can be seen for example in [47]. This offers a great flexibility to deal with different shapes or to change easily the physical properties for each immersed structure. Therefore, we start by computing the signed distance function (level set) of a given geometry to each node of the mesh. Using the zero isovalue of this function, we can easily identify the fluid-solid interface. Consequently, different parts are immersed in a larger domain of different constitutive
laws. This part will only briefly explain the immersed approach process, but more details can be found in [48].

At any point $\mathbf{x}$ of the computational domain $\Omega$, the level-set function $\alpha$ corresponds to the signed distance from $\Gamma_{im}$. In turn, the interface $\Gamma_{im}$ is given by the zero isovalue of the function $\alpha$:

$$ 
\begin{aligned}
\alpha(\mathbf{x}) &= \pm d(\mathbf{x}, \Gamma_{im}), \mathbf{x} \in \Omega, \\
\Gamma_{im} &= \{\mathbf{x}, \alpha(\mathbf{x}) = 0\}.
\end{aligned} 
\quad (1) 
$$

In this paper, the following sign convention is used: $\alpha \ge 0$ inside the solid domain defined by the interface $\Gamma_{im}$ and $\alpha \le 0$ outside this domain.

As explained, the signed distance function is used to localize the interface of the immersed structure but it is also used to initialize the desirable properties on both sides of the latter. Indeed, for the elements crossed by the level-set functions, fluid-solid mixtures are used to determine the element effective properties. A Heaviside function $H(\alpha)$ is then defined as follow:

$$ H(\alpha) = \begin{cases} 1 & \text{if } \alpha > 0 \\ 0 & \text{if } \alpha < 0 \end{cases} \quad (2) $$

In the numerical approximation of the problem we will consider a partition of the computational domain $\Omega$ with a finite element mesh made by a collection of element $\{K\}$.

The Heaviside function can be smoothed to obtain a better continuity at the interface [?] using the following expression:

$$ H_{\varepsilon}(\alpha) = \begin{cases} 1 & \text{if } \alpha > \varepsilon \\ \frac{1}{2} \left( 1 + \frac{\alpha}{\varepsilon} + \frac{1}{\pi} \sin \left( \frac{\pi \alpha}{\varepsilon} \right) \right) & \text{if } |\alpha| \le \varepsilon \\ 0 & \text{if } \alpha < -\varepsilon \end{cases} \quad (3) $$

where $\varepsilon$ is a small parameter such that $\varepsilon = O(h_{im})$, known as the interface thickness, and $h_{im}$ is the mesh size in the normal direction to the interface. In the vicinity of the interface, it can be computed using the following expression:

$$ h_{im} = \max_{j,l \in K} \nabla \alpha \cdot \mathbf{x}^{jl}, \quad (4) $$

where $\mathbf{x}^{jl} = \mathbf{x}^l - \mathbf{x}^j$ and $K$ is the mesh element under consideration. According to the chosen approximations, the Heaviside function is then approximated using linear interpolations (P1) between fluid and solid properties or a piecewise constant interpolation (P0).

Finally, we combine this approach with an anisotropic mesh adaptation algorithm to ensure an accurate capturing of the discontinuities at the fluid-solid interface. However, the levelset function intersects the mesh element arbitrarily. It is possible then to overtake the discontinuity appearing at the interface by using mesh adaptation and regularization. The regularization parameter can be seen as the thickness or the resolution of the interface. It is shown that using local adaptivity, stretched elements at the interface are obtained and one can then choose a priori the resolution of the thickness to be very small, which leads to very sharp interfaces, favorable for simulating fluid-structure conduction problems.

# 3. SURFACE TO SURFACE RADIATION

## 3.1. S2S hypotheses

We will now recall the assumptions under which the S2S model is considered valid.

1. Surface to surface radiation, as its name indicates, handles radiation as a purely surface phenomenon. This means radiation is considered to be utterly absorbed into the bulk of the radiating elements, and emitted, within a characteristic length of penetration that is small in regard of the dimension of the
system. Such radiating elements are said to be opaque. Most of the production of radiation occurs within 20 nm of the surface [49]. For radiation absorption, Wien's law tells us the wavelength $\lambda_m$ at which a blackbody emits its maximum amount of radiation is dependent on its surface temperature T [50]:

$$ \lambda_m T = 2897,9 \ \mu m \cdot K \qquad (5) $$

Then, 95% of the total amount of radiant energy emitted by the blackbody stands in a wavelength band of 0.5 $\lambda_m - 5 \lambda_m$. For most industrial applications, ambient temperature can vary between 300 K and 2500 K, which amounts to a wavelength band of 0.6 - 50 $\mu m$ (from visible to far infrared wavelengths). One can infer from [51] that for most metals, complete absorption of visible and infrared radiation occurs within 0, 1 $\mu m$ of the surface. For other opaque (in the visible range) materials such as silicon nitride ceramics [52] or zirconium-yttrium ceramics [53] skin thickness for visible and near infrared wavelengths is of the order of 100 $\mu m$. This is a negligible thickness for industrial products for which typical detail size is that of 1 cm. Most liquids and gases are not opaque to radiation however, but the S2S model only considers thermal radiation between solid surfaces.

2. This brings the next hypothesis: the medium between the solid, radiating surfaces is considered transparent to radiation. Obviously this is true when the medium is vacuum, but this also holds true for diatomic gases who do not possess a polar moment ($N_2$, $O_2$, etc) [54]. This means for example dry air is actually transparent to radiation. The medium does not absorb, nor reflects any radiation, but fully transmits it. It is possible to extend the S2S method to the case of non-transparent media, in which case it is called the zonal method [55].

3. Radiating surfaces are considered diffuse. This means that an incoming ray of energy will be reflected isotropically, with no direction of reflection carrying a greater fraction of the incoming energy than the others. In addition, blackbody radiation is considered to be emitted isotropically. This amounts to considering that emissivity and reflectance are independent of direction. This assumption is needed since obtaining detailed direction-dependant radiation measures is fairly difficult. In reality, these values do vary a lot with direction: see for example [56] for tables of reflectance for various materials. However, in closed systems, the multiple reflections that occur between all surfaces tend to even out any direction-dependant behaviour. This argument makes the diffuse surfaces hypothesis reasonable.

4. The last hypothesis we make is also the strongest one. We assume emissivity, absorptance and reflectivity are not only independent of direction (which is a reasonable assumption to make, as previously explained), but are in fact constant. In reality, these values depend on a great deal many parameters, many of which constantly evolve during the industrial process: temperature, wavelength of radiation emission, surface rugosity, oxydation status, and even surface stress ([57] offers a good review of this situation). Not only is it nearly impossible to have access in real time to the surface state of all the radiating items during the process, but also the relationship between local emissivity and all these variables is, most often, not known for a given material. Because of this, we decide to take the simplification of constant radiative values. In particular, this implies considering emissivity does not depend on emitted wavelength. This amounts to make the "grey body assumption", where the graph of radiant energy depending on wavelength for any given surface is identical to that of a blackbody at the same temperature, multiplied by the constant emissivity standing between 0 and 1.

Then, to compute immersed radiative transfer, we need more than an implicit definition of the interface $\Gamma_1$. The first step of the S2S algorithm is to reconstruct a list of explicit facets out of the immersed mesh to determine where radiation effectively takes place.

### 3.2. Facet list construction

S2S methods used in the context of body-fitted meshes usually have a direct access to the list of facets (that is, lines in 2D and triangles in 3D) that contributes to radiation. Indeed, these facets are the border elements of the mesh of the domain and of the mesh of the solid items. When used in combination with our immersed mesh, all participating facets are not border elements: the facets on the solid interface $\Gamma_1$ actually belong to bulk elements, so there is no easy access to the list of elements containing the participating facets.
Moreover, in immersed meshes, participating facets do not always (and in fact, rarely do) coincide with existing element facets.

<figure><img src="image_4.png" /><figcaption>Figure 1: 2d intersections of elements by the solid interface</figcaption></figure>

The first step is thus to determine all the elements that are intersected by the solid interface $\Gamma_1$. Figure 3.2 shows the multiple cases that can be encountered in 2D. For each element of the mesh, the solver evaluates the levelset at each vertex of the elements. Several cases then arise:

1. The levelset takes a 0 value at a single vertex of the element, and takes either all-negative or all-positive values at all the other vertices (case 1 in Figure 3.2). In such a case the solver does not add anything to the list of participating facets.

2. The levelset takes a 0 value at exactly two vertices of the element (case 2 in Figure 3.2). This means one of the edges of the element is indeed on the interface $\Gamma_1$, and this edge is added to the list of participating facets.

3. The levelset takes a 0 value at one vertex, and two values of different signs at the other two (case 3 in Figure 3.2). In that case the solver creates a virtual vertex at the position of the intersection between the solid interface and the edge of the element between the two vertices at which the levelset takes a nonzero value. It then creates a virtual facet anchored between the 0-value vertex and the virtual vertex, and adds it to the list.

4. The levelset takes three nonzero values, not all of them being of the same sign (case 4 in Figure 3.2). This means the solid interface crosses the element right through the middle. Similarly to what happens in case 3, the solver creates two virtual vertices and then a virtual facet between them. This facet is then added to the list.
<figure><img src="image_5.png" /><figcaption>Figure 2: Relative position of elementary surfaces $dS_i$ and $dS_j$.</figcaption></figure>

5. The levelset function takes all-negative or all-positive values at all vertices. This means the solid interface does not cross the element at all.

The outer shell of the domain $\Omega$ is also defined by a levelset function, and the facets of the border of the domain that contribute to radiation are also computed by interface detection. We could have simply retrieved the flat elements composing the border of the domain, but detecting them by levelset as well ensures maximum generability of the method. This way it is easy to modify the dimensions of the outer radiating object (ie, a room's walls) on the fly without having to remesh the domain.

### 3.3. View factor computation

The next step is the computation of the so-called view factor matrix. The intuitive definition of the view factor $F_{ij}$ is that $F_{ij}$ is equal to the proportion of the total radiation emitted by facet $F_i$ that will strike facet $F_j$. Its mathematical definition is as follows: Let there be two facets $F_i$ and $F_j$ of respective surfaces $S_i$ and $S_j$, normals $\vec{n}_i$ and $\vec{n}_j$, separated by a vector $\vec{r}_{ij}$. These normals are such that the angles between $\vec{n}_i$ and $\vec{r}_{ij}$ and $\vec{n}_j$ and $\vec{r}_{ij}$ are, respectively, $\Theta_i$ and $\Theta_j$. This setup is displayed with elementary 2D surfaces Figure 2. Then:

$$F_{ij} = \frac{1}{S_i} \int_{S_i} \int_{S_j} \frac{\cos(\Theta_i) \cdot \cos(\Theta_j)}{r_{ij}^2} dS_j dS_i \quad (6)$$

This adimensional number depends on geometrical parameters only, and needs to be calculated only once at the first time step of the simulation (except if at some point any radiating item moves inside the system). Filling the view factor matrix is, however, an arduous process. The view factor is quite a computationnaly expensive term to calculate, because of its double integral nature. Moreover, the radiosity matrix is rarely sparse, so a system involving $N$ facets will have to calculate $O(N^2)$ terms. As a result, finding various ways to accelerate the computation of the view factors is in and of itself a dense topic of research [58, 59, 46, 60]. First, three properties of the view factor allow for a more efficient computation. Given a closed system containing $N$ facets, for any given facets $F_i$ and $F_j$ of respective surfaces $S_i$ and $S_j$:

$$F_{ii} = 0 \qquad \text{no self-radiation} \qquad (7)$$

$$S_i F_{ij} = S_j F_{ji} \qquad \text{reciprocity of radiation} \qquad (8)$$

$$\sum_{j=1}^{N} F_{ij} = 1 \qquad \text{conservation of energy} \qquad (9)$$

The first advantage of these formulae is that thanks to Equation 6, one only has to compute half of all the view factors and deduce the other half through reciprocity of radiation. Moreover, for any one give facet $F_i$, we only consider the facets $F_j$ such that $\vec{n}_i$ points towards $F_j$ and $\vec{n}_j$ points towards $F_i$. For all other cases, this means at least one of the facets is "turning its back" to the other, and as such no radiation should be considered and the view factor is set to 0. After this, several ways exist to compute the view
<figure><img src="image_6.png" /><figcaption>Figure 3: Crossed strings method to determine the view factor between two surfaces in 2D</figcaption></figure>

factors in an efficient way. For 2D geometries, we have already mentioned the crossed-string method [42], that provides an analytical answer to the problem. Given two surfaces (even if of curved shape) A₁ and A₂, of temperature T₁ and T₂, formed of the lines (ab) and (cd), then:

$$A_1 F_{12} = A_2 F_{21} = \frac{ad + bc}{2(ac + bd)} \qquad (10)$$

Figure 3 explains the process. This allows quick computations of 2D view factors. As far as 3D geometries are concerned, several approaches are possible. Monte-Carlo ray tracing [33, 46, 60] or the hemicube method [31, 59] are some of the existing methods, but this paper will focus on a semidirect integration approach [43, 31]. Indeed, Equation 6 is a double surface integral. A direct numeric integration is possible, but quite costly. A more efficient approach is to turn the integral into a double line integration over the edges of each facet, thanks to the Stokes theorem [43]. We then numerically integrate only over the edges of the facet, which is less expensive. We split the contours of the facets into nᵢ and nⱼ, and have the following:

$$F_{ij} = \frac{1}{4\pi S_i} \sum_{k=1}^{n_i} \sum_{l=1}^{n_j} \int_{e_k^i} \int_{e_l^j} \ln(d(e_k^i, e_l^j)) u_k v_l du_k dv_l \qquad (11)$$

where eᵢⱼ is the k-th segment making up the discretized contour of Fᵢ, d(eᵢⱼ, eⱼⱼ) is the Euclidean distance between the two segments being integrated and vᵢ and vⱼ are the integration variables. A Gauss quadrature is then used to compute the value of the logarithm function. As one can see, a problem arises when two facets sharing an edge are analyzed, as the logarithm will be evaluated at 0. A different formula is used in this case [61]. For two elementary surfaces sharing an edge dl:

$$d^2 F_{ij} = \frac{dl^2}{2} (3 - \ln(dx^2)) \qquad (12)$$

### 3.4. Obstruction handling

The last parameter to consider when computing the view factors between two facets is the possibility that a third facet is located between these two and blocking off, totally or partially, emitted radiation from one to the other. Obstruction handling is very important, because if properly carried out, it allows to assign 0 to a number of view factors without having to compute them. This also makes the view factor matrix more sparse. Finally, is it crucial for conservation of energy as a whole that radiation be not computed where none should be received. Obstruction handling is a costly operation, because at first approximation, for each of the n² view factors Fᵢⱼ in a system of n facets, there are O(n) possibly obstructing facets to be
<figure><img src="image_7.png" /><figcaption>Figure 4: Obstruction check for $F_k$ between $F_i$ and $F_j$. $B_k$ is not within the projected tetrahedra and thus $F_k$ does not obstruct $F_{ij}$.</figcaption></figure>

checked. That means view factor computation as a whole takes now $O(n^3)$ operations. As a result, many workarounds are deployed to try and diminish the required number of operations, and some papers focus specifically on this point [60].

1. For closed systems of a convex shape, the facets that make up the border of the system cannot be in the way of any possible couple of other facets. For nonconvex systems, one can still compute, for any given border facet $B_i$, the plane in which $B_i$ belongs, and then if all other facets belong in only one of the two half-spaces defined by this plane (that is, whether all other facets are "above" or "below" $B_i$). Such facets can be removed from the list of facet checked for obstruction.

2. Starting from facet $F_i$ and to look if a possible facet $F_k$ is in the way between $F_i$ and $F_j$, if facet $F_k$ does obstruct, then there is a solid item of nonzero width in the way. This means there is facet $F_l$ on the other side of the obstructing object, and of a normal opposed to that of $F_k$, that will obstruct too. For this reason one only has to check for the $F_k$ that have a normal in the same sense as $F_i$. If $F_k$ has a normal opposite to that of $F_i$ (an operation that is much less computationally costly than checking for the specific relative geometric positions of each facet), then it can be dismissed at once, knowing at some point in the list facet $F_l$ will be encountered.

3. Facet $F_k$ must lie in the portion of space between the planes containing $F_i$ and $F_j$.

The biggest simplification made with obstruction, however, is the following: obstruction is considered to be a boolean value. Either facet $F_k$ completely obstructs the pair $F_i, F_j$, and we have $F_{ij} = 0$; or it does not obstruct at all, and we proceed to check for $F_{k+1}$. Theoretically, one should account for partial obstructions, and devise a corrected view factor accordingly. However, such an operation is extremely computationally costly for a modest benefit in accuracy if the facets are small enough. When the mesh is refined enough, the resulting facets tend to be very small before the objects immersed in the mesh, so that in practice all obstructions are all-or-nothing.

Obstruction is defined by the following: when considering facets $F_i, F_j$ and $F_k$, the barycentres $B_i, B_j$ and $B_k$ are computed. Then, vectors from $B_i$ to each of the vertices of $F_j$ are computed. This defines a triangle in 2D, its edges being the two such drawn vectors and facet $B_j$ itself; or a tetrahedra in 3D, its edges being the three such drawn vectors and the edges of facet $B_j$ itself. Then, we check whether barycentre $B_k$ is inside this triangle or tetrahedra. If it is not, we then construct the vectors from $B_j$ to the vertices of $F_i$, define again a triangle or tetrahedra, and check again whether $B_k$ is within it or not. If the answer is still no, we then consider facet $F_k$ to be too offset from the path of radiation from $B_i$ to $B_j$ to obstruct it. $F_k$ is considered non obstructing, and we move on to check for $F_{k+1}$. Figure 4 sums up this process.
### 3.5. Flux computation

Once the view factor matrix has been computed, added with the temperature and emissivity of all the facets, we have all the data we need to compute the radiation emitted by each facet. A useful intermediate value is introduced, which is called radiosity, noted $J$. It stands for the sum of all surface radiation leaving the facet, regardless of radiation that is received by it. This means that facet $F_j$ receives from $F_i$ (of surface $S_i$) a total amount of radiation equal to $S_i F_{ij} J_i$. In return, $F_i$ receives from $F_j$ an amount of radiation equal to $S_j F_{ji} J_j$. Thanks to the rule of reciprocity (Equation 8) we can rewrite this as $F_i$ receiving $S_i F_{ij} J_j$. Radiosity $J_i$ of facet $F_i$ with temperature $T_i$ and emissivity $\epsilon_i$ is thus implicitly defined to be equal to the sum of the flux emitted by the facet by black body radiation and the reflected part of all radiation received from other facets.

$$S_i J_i = S_i \epsilon_i \sigma T_i^4 + (1 - \epsilon_i) \sum_j F_{ji} J_j \quad (13)$$

Equation 13 shows that the vector of radiosities can be deduced from a matrix inversion with some reformulating:

$$J_i - (1 - \epsilon_i) \sum_j F_{ij} J_j = \epsilon_i \sigma T_i^4 \Leftrightarrow \mathbf{M}\mathbf{J} = \mathbf{T} \quad (14)$$

where matrix $\mathbf{M}$ and vector $\mathbf{T}$ are defined as follows:

$$M_{ij} = \begin{cases} 1 & \text{if } i = j \\ -(1 - \epsilon_i)F_{ij} & \text{if } i \neq j \end{cases} \quad (15)$$

$$T_i = \sigma \epsilon_i T_i^4 \quad (16)$$

As a result, we simply need to invert matrix $\mathbf{M}$. It is not very sparse: common orders of magnitude for the ratio of nonzero values are between 40% and 60%. That means a direct inversion would be quite computationally costly. However, it also is a symmetric positive, diagonally-dominant matrix. These properties make it a good candidate for a GMRES resolution method, which under these circumstances turns computational time from $O(N^2)$ to $O(N)$ operations for a matrix of order $N$. Once the matrix is inverted, we obtain the radiosities vector. A final formula allows us to extract the net outgoing surface radiative flux from the radiosity:

$$\phi_i = \frac{\epsilon_i}{1 - \epsilon_i} (\sigma T_i^4 - J_i) \quad (17)$$

A final check is made at this stage. Due to energy conservation, we should have $\sum_i S_i \phi_i = 0$. It happens often enough that the sum of the fluxes is actually a non-negligible quantity, due to inaccuracies in view factor computation or in matrix inversion. To ensure no energy is either created or lost, we define a vector of secondary fluxes through a very simple correction:

$$\phi'_i = \phi_i - \frac{1}{S_i N} \sum_j S_j \phi_j \quad (18)$$

### 3.6. Coupling with thermal transfer

Now that we have computed the list of the net outgoing surface radiative flux for all the facets in the system, let us see how we can plug this data in the general resolution of thermal transfer. We use a convection-diffusion (CD) equation, which is a volumic equation that solves thermal transport for a given P1 element.

$$\rho c_p \frac{\partial T}{\partial t} + \rho c_p \vec{v} \cdot \vec{\nabla}(T) - \lambda \nabla \cdot (\vec{\nabla} T) = Q \quad (19)$$

where:

* $\lambda \nabla \cdot (\vec{\nabla} T)$ stands for the diffusive term. $\lambda$ is the conductivity in the given element.
* $\rho c_p \vec{v} \cdot \vec{\nabla}(T)$ stands for the convective term. $v$ is the speed of the fluid (if any) at the boundaries of the element. $\rho$ and $c_p$ are the volumic mass and the specific heat capacity.

* $Q$ stands for the source term. It accounts for any input (or output) of energy, of which amount has already been calculated earlier in the time step or is constant with regard to the temperature (like energy produced by radioactivity).

All values in the above equation are volumic in nature. That explains why radiation does not have a term directly associated to it: in the context of the S2S modelisation, it is an inherently surface phenomenon. In order to account for the energy gained or lost by the element through radiation, we compute beforehand how much net absolute energy is brought by radiation to the edges of the cell. We then divide that absolute amount by the volume of the element, so that we obtain an equivalent volumic radiative terme $Q_r$. That term is then added to Equation 19 as a source term. The underlying physical hypothesis is that we essentially consider radiation to be stationary when compared to convection or diffusion. Indeed, at each new time step, we take as input the field of temperatures calculated at the end of the previous time step, we first solve the radiation with the S2S method using that temperature, and we then use the radiation field as an input for Equation 19 after the aforementioned conversion. This is justified, however, by the fact that the characteristic time of radiation is $1/c$, much smaller than the characteristic times of convection or conduction.

# 4. NUMERICAL RESULTS

We present here the results of the simulations produced by our solver coded as explains beforehand. We conduct comparisons with analytical test cases. We first recall the set of assumptions under which the S2S modelling takes place in this paper:

1. All solid are fully opaque to radiation.

2. Emissivity is a constant for a given surface.

3. Radiating surfaces are diffuse.

4. The medium between radiating surfaces is fully transparent to radiation.

## 4.1. Simple radiation

The first step is to verify our solver performance in cases of purely radiative thermal transfer. To this end, we consider several simple geometries, both in 2D and in 3D, for which analytical solutions can be found.

### 4.1.1. 2D test cases: concentric cylinders

We will evaluate radiation between two concentric circles separated by a nonparticipating medium. Since theory tells us the net outgoing surface radiative flux should be equal all along the surface of each of these two circles for symmetry reasons, we can make detailed analyses of the accuracy of our algorithm. The two circles have a respective radius of $r_{in} = 0.4$ and $r_{out} = 0.6$ with an emissivity $\epsilon_{in} = \epsilon_{out} = 0.5$ for both surfaces. We set $T_{in} = 500$ K and $T_{out} = 1000$ K. The mesh is uniform with a mesh size of 0.02, which amounts to a total of 8756 elements, and out of these 476 radiating facets. Figure 5 displays the setup.

In such a case, the only physical phenomenon that is considered is thermal radiation between two concentric circles at a fixed temperature. When steady state is reached, the net exchanged radiation is null (since the system is closed). The total amount of radiation $\Phi_{tot,i \to o}$ leaving the inner circle and absorbed by the outer circle is (see [50]) with $\sigma$ standing for the Stefan constant:

$$ \Phi_{tot,i \to o} = \frac{2\pi r_{out} \sigma (T_{in}^4 - T_{out}^4)}{(1/\epsilon_{in} \cdot r_{out}/r_{in} + 1/r_{out} - 1)} \quad (20) $$

The amount of radiation going from the outer circle to the inner circle is the opposite of this value. Dividing by the respective surface area of each circle gives the analytical value of the surface flux to be expected on each point of the surface of each respective circle, $\phi_{s,i \to o}$ and $\phi_{s,o \to i}$:
<figure><img src="image_8.png" /><figcaption>Figure 5: Pure radiation test case in concentric circles. The leftmost figure displays the enforced temperatures on the boundaries. The central figure displays the mesh employed. The rightmost figure displays the computed radiative fluxes at the radiating surfaces.</figcaption></figure>

$$ \left\{ \begin{aligned} \phi_{s,i \to o} &= \frac{\Phi_{tot,i \to o}}{2\pi r_{in}} = \frac{r_{out}\sigma(T_{in}^4 - T_{out}^4)}{r_{in}(1/\epsilon_{in} \cdot r_{out}/r_{in} + 1/r_{out} - 1)} \\ \phi_{s,o \to i} &= \frac{\Phi_{tot,o \to i}}{2\pi r_{out}} = \frac{r_{out}\sigma(T_{in}^4 - T_{out}^4)}{r_{out}(1/\epsilon_{in} \cdot r_{out}/r_{in} + 1/r_{out} - 1)} \end{aligned} \right. \qquad (21) $$

We then compare the range of outgoing surface radiative flux obtained over the inner and outer circles, and compare with the analytical values we deduce from Equation 21. First, we compute the average surface flux, obtained by integrating the computed local surface flux values over the circumference of each circle. These average values should match the analytical surface flux value. Then, in order to determine dispersion, for each circle we look up the local surface flux value of the mesh point where this value is farthest away from the expected value out of all mesh points on that circle. We then compute the error between that value and the expected surface value. This gives us the maximum error committed by the solver at any given point on each circle. Table 1 sums up the results. We can see that the obtained results are within 1.2% of the theoretical values.

<table><caption>Table 1: Comparison of net output radiative flux with analytical solution for concentric circles pure radiation case</caption><thead><tr><th>Source</th><th>Inner circle surface flux (W·m<sup>-1</sup>)</th><th>Outer circle surface flux (W·m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Analytical result</td><td>-19935</td><td>13290</td></tr><tr><td>Present work</td><td>-19877 (-0.3%)</td><td>13401 (+0.8%)</td></tr><tr><td>Maximum error value</td><td>-19876 (-0.3%)</td><td>13445 (+1.2%)</td></tr></tbody></table>

### 4.1.2. 2D test cases: square within a square

We compute pure radiation between two homothetic squares separated by a nonparticipating medium. The outer square has an edge size of 1 and is set to $T_{out} = 1000$ K. The inner square has an edge size of 0.5 and is set to $T_{in} = 500$ K. Emissivity is set at $\epsilon = 0.5$ for both squares. The meshing is homogeneous with a mesh size of 0.04. The number of elements is 2314, with 161 facets participating in radiation. The setup can be seen on Figure 6. In this simple configuration, we can not find an analytical solution describing the surface radiative flux exchanged at every point of each of the two squares. However, we can find a solution
<figure><img src="image_9.png" /><figcaption>Figure 6: Pure radiation test case in homothetic squares</figcaption></figure>

to the total radiative flux exchanged by each square towards the other. For this, we will consider each side of the squares to be a single emitting facet, for a total of 8 facets in the system. Analytical formulations of the view factors in such a case can be processed by hand with the use of the crossed strings methods as per Equation 10 and displayed on Figure 3.

The view factor matrix ($F_{ij}$) is displayed in Equation 22. From here on, we then invert by hand the 8x8 radiosity matrix $M$ as seen in Equation 15, and from this we deduce the radiosity vector $J$ and finally the flux $\phi$ of all radiating facets as per Equation 17.

$$ (F_{ij}) = \begin{pmatrix} 0 & 1 - \frac{\sqrt{10}}{4} & 1 - \frac{\sqrt{10}}{4} & \frac{\sqrt{10}-3}{2} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & \frac{\sqrt{10}-\sqrt{2}}{4} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & 0 \\ 1 - \frac{\sqrt{10}}{4} & 0 & \frac{\sqrt{10}-3}{2} & 1 - \frac{\sqrt{10}}{4} & \frac{\sqrt{10}-\sqrt{2}}{4} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} \\ 1 - \frac{\sqrt{10}}{4} & \frac{\sqrt{10}-3}{2} & 0 & 1 - \frac{\sqrt{10}}{4} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & \frac{\sqrt{10}-\sqrt{2}}{4} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} \\ \frac{\sqrt{10}-3}{2} & 1 - \frac{\sqrt{10}}{4} & 1 - \frac{\sqrt{10}}{4} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{4} & \frac{\sqrt{10}-\sqrt{2}}{4} \\ \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & \frac{\sqrt{10}-\sqrt{2}}{2} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & 0 & 0 & 0 \\ \frac{\sqrt{10}-\sqrt{2}}{2} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & 0 & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & 0 & 0 & 0 \\ \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & 0 & \frac{\sqrt{10}-\sqrt{2}}{2} & \frac{1-\frac{\sqrt{10}-\sqrt{2}}{2}}{2} & \frac{\sqrt{10}-\sqrt{2}}{2} & 0 & 0 & 0 \end{pmatrix} \quad (22) $$

The analytical calculation gives us the expected total radiative flux integrated over each of the four facets (that is, each side) of each square, but cannot give us the expected surface radiative flux at any given point of the geometry. For this reason, we can only compare the surface radiative flux of our computed solution integrated over the four facets of each square to the analytical expected value. We thus compare global surface flux to global surface flux, and not surface values as was the case for the concentric circles.
The results are indicated on Table 2. Once again the computed solution stands very close to the analytical solution —though we cannot get as many accuracy criteria as in the concentric circles case.

<table><caption>Table 2: Comparison of net output radiative flux with analytical solution for concentric circles pure radiation case</caption><thead><tr><th>Source</th><th>Inner square integrated flux (W)</th><th>Outer square integrated flux (W)</th></tr></thead><tbody><tr><td>Analytical result</td><td>-42528</td><td>42528</td></tr><tr><td>Present work</td><td>-42139 (-0.9%)</td><td>42140 (-0.9%)</td></tr></tbody></table>

### 4.1.3. 3D test cases: cube within a cube

We now present a test case in pure radiation for the following 3D geometry: a cube of edge size 0.3 within a cube (which is a homothetic transformation of the first one) of edge size 0.6. The inner cube is set at temperature $T_{in} = 1000$ K and the outer cube at $T_{out} = 500$ K. Emissiviy is set at $\epsilon = 0.5$ for both cubes. The solver computes a total of 164 participating facets. The geometry is shown Figure 7. We measure the total net surface radiative flux integrated on respectively the inner cube and the outer cube. For this simple case, in order to compute the analytical solution of the problem, we need to compute 3D view factors between the facets of the two cubes. For this we will consider each cube to be composed of only 6 radiating facets: this means we will consider each side of each cube to be a single unique radiating facet. This way, our system is composed of 12 radiating facets, and thus we only need to set up a view factor matrix of size 12. In 3D, we can not use the crossed strings method to compute the view factor. However the "cube within cube" view factor 3D integrals of Equation 6 have been analytically solved in the literature: for example in Howell [62]. From these values we can once again manually set up and invert the radiosity matrix, and thus obtain the analytical radiative flux received by each side of each cube.

Our computation gives us the radiative flux received by each of the 164 facets in the mesh, and we then integrate the value from all the facets composing each individual cube side, so that we obtain the total computed flux received by each side of each cube. We then compare these values to those devised analytically. Table 3 shows the obtained results. The simulated fluxes are very close from the analytical fluxes.

<table><caption>Table 3: Comparison of net output radiative flux with analytical solution for cube within cube pure radiation case</caption><thead><tr><th>Source</th><th>Inner cube flux (W·m<sup>-2</sup>)</th><th>Outer cube flux (W·m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Analytical result</td><td>-21747</td><td>9665</td></tr><tr><td>Present work</td><td>21755</td><td>9669</td></tr><tr><td>Relative difference (%)</td><td>0.04</td><td>0.04</td></tr></tbody></table>

### 4.2. Conduction-radiation

Now that we are confident with our radiation results for the considered geometries, we can simulate a multiphysics heat transfer problem. Once again we will first provide 2D examples before moving on to 3D test cases.

#### 4.2.1. 2D test cases: infinitely long rods

We will resolve thermal conduction, coupled with radiation as a source term, on the following geometry: two parallel, infinitely long solid rods of thickness $e_{rod} = 1$ separated by a gas channel of thickness $e_{gas} = 3$. There is no fluid circulation, meaning the convection term of Equation 19 is zero. The bottom side of the bottom rod sees a boundary condition $T_1 = 1000$ K while the top side of the top rod sees $T_4 = 0$ K.
<figure><img src="image_10.png" /><figcaption>Figure 7: Pure radiation test case in cube within cube</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Mesh used for the infinitely long rods. It has been refined around the radiating surfaces.</figcaption></figure>

Conduction operates through the width of the bottom rod, then through the air, and finally through the width of the top rod. The gas itself is considered transparent to radiation. Meanwhile, radiation also occurs between the top side of the bottom rod and the bottom side of the top rod. The chosen values for emissivity and conductivity are arbitrary: we take $\epsilon = 0.5$ for all radiating surfaces, $\lambda_{rod} = 25 \text{ W} \cdot \text{K}^{-1} \cdot \text{m}^{-1}$ for the conductivity in the steel and $\lambda_{gas} = 10 \text{ W} \cdot \text{K}^{-1} \cdot \text{m}^{-1}$ for the conductivity in the air. Volumic masses are $\rho_{rod} = 7700 \text{ kg} \cdot \text{m}^{-3}$ and $\rho_{gas} = 1 \text{ kg} \cdot \text{m}^{-3}$ respectively, and specific heat capacity are $cp, rod$ and $c_{p,gas} = 1000 \text{ J} \cdot \text{kg}^{-1} \cdot \text{K}^{-1}$ respectively. Remeshing takes place around the radiating facets, as explained in Part 1. We can afford a coarser mesh in the bulk of each domain, where only conduction takes place. We start from a uniform mesh of 37816 elements, reduced to 5642 elements after remeshing, with 1004 radiating facets. The setup is displayed on Figure 9, while the mesh can be seen on Figure 8. It amounts to a 1D-thermal transfer problem for which an analytical solution can be devised in steady state. In particular, the temperatures $T_2$ and $T_3$ can be measured on the simulation after convergence and compared to the analytical solution. The total surface radiative flux crossing the air zone between the rods can also be retrieved and compared to the analytical solution. Figure 11 shows thermal propagation and radiative flux in the rods at various timesteps. When the computation begins, only the bottom surface is at a temperature other than OK. In theory this means the radiating surfaces should output an absolute radiative flux of $0 \text{W} \cdot \text{m}^{-2}$. Due to numerical inaccuracies, very small non-zero fluxes are computed instead, though they do not adversely impact thermal propagation. After 10 time steps we can see the radiating surfaces have reached a nonzero temperature and start radiating flux towards the other. After 50 time steps, steady state is reached, and we can use the temperatures obtained at that stage for analytical comparison.

In order to devise the analytical solution to the problem, we find it convenient to use the electricity analogy. Figure 10 shows a graphical representation. Our 1D system is composed of four different potential points: $T_1, T_2, T_3$ and $T_4$. Out of these, $T_1$ and $T_4$ are known. $T_1$ and $T_2$, and $T_3$ and $T_4$ respectively, are distant from a thickness $e_{rod}$ of conductivity $\lambda_{rod}$, so the thermal conduction resistance between these
<figure><img src="image_12.png" /><figcaption>Figure 9: Conduction–radiation test case. The mesh has been refined around the radiating surfaces.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Electrical analogy of thermal transfer through the rods</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11: Thermal propagation and radiative flux in the rods. The first line shows the initial setup. The second line shows the state after 10 timesteps. The third line shows the state after 50 timesteps, when steady state is reached.</figcaption></figure>

potentials is $R_{cond,1→2}$ and $R_{cond,3→4}$ respectively. Since the walls facing each other are of infinite surface (in order to ensure the physics occurs in 1D), surface values will be considered. All terms of surface normally appearing in the formulae will be replaced by a set surface value of $1m^2$. The expressions of the resistances are thus defined as:

$$R_{cond,1→2} = R_{cond,3→4} = \frac{e_{rod}}{\lambda_{rod}} \qquad (23)$$

These resistances are expressed in $K \cdot m^2 \cdot W^{-1}$, which is the inverse of $W \cdot K^{-1} \cdot m^{-2}$, which is indeed the unit of a thermal conductivity. Then, $T_2$ and $T_4$ are connected by two parallel branches. The branch accounting for the conduction thermal flux holds a thermal conduction resistance $R_{cond,2→3}$ such as:

$$R_{cond,2→3} = \frac{e_{gas}}{\lambda_{gas}} \qquad (24)$$

And the branch accounting for the radiation thermal flux holds a thermal radiation resistance $R_{rad,2→3}$ the expression of which can be deducted as follows. According to [50], the thermal surface exchanged radiative flux $\phi$ between two infinite, parallel walls of respective temperature $T_1$ and $T_2$ and emissivity $\epsilon_1$ and $\epsilon_2$ is

$$\phi = \frac{\sigma(T_1^4 - T_2^4)}{\frac{1}{\epsilon_1} + \frac{1}{\epsilon_2} - 1} \qquad (25)$$

which can be rewritten as the expression of a potential difference divided by a thermal resistance as follows:

$$\phi = \frac{(T_1 - T_2)}{\frac{\frac{1}{\epsilon_1} + \frac{1}{\epsilon_2} - 1}{\sigma(T_1 + T_2)(T_1^2 + T_2^2)}} \qquad (26)$$

We thus obtain a thermal radiation resistance the expression of which depends on both $T_1$ and $T_2$:

$$R_{rad,1→2}(T_1, T_2) = \frac{\frac{1}{\epsilon_1} + \frac{1}{\epsilon_2} - 1}{\sigma(T_1 + T_2)(T_1^2 + T_2^2)} \qquad (27)$$

We can observe that its value decreases as $T_1$ and $T_2$ increase, when the difference between the two is fixed. This makes sense, as the difference between the fourth power of the two increases, making the radiative exchange more important and thus the resistance smaller.
Continuing the electricity analogy, we can define an equivalent thermal resistance $R_{eq}$ from $T_1$ to $T_4$ that takes into account the four previously defined resistances. Here, we take $\epsilon_1 = \epsilon_2 = \epsilon_{rod}$ (since the gas is transparent to radiation, its emissivity does not intervene).

$$R_{eq} = R_{cond,1 \to 2} + \frac{1}{\frac{1}{R_{cond,2 \to 3}} + \frac{1}{R_{rad,2 \to 3}}} + R_{cond,3 \to 4} \quad (28)$$

When steady state is reached, since there is no local source or sink of flux, total thermal flux $\Phi_{tot}$ is constant throughout the system. We can thus write Ohm's law from one end of the system at $x = 0$ at $T = T_1$ to the other at $x = 5$ and $T = T_4$:

$$\Delta T = T_1 - T_4 = R_{eq} \Phi_{tot} \Leftrightarrow \frac{T_1 - T_4}{R_{eq}} - \Phi_{tot} = 0 \quad (29)$$

However, $\Phi_{tot}$ is not the only unknown here, since the expression of $R_{eq}$ includes $R_{rad,2 \to 3}$, which itself depends upon both $T_2$ and $T_3$. We will thus express $T_2$ and $T_3$ from $\Phi_{tot}$ by writing Ohm's law between $x = 0$ at $T = T_1$ and $x = 1$ at $T = T_2$, and $x = 4$ at $T = T_3$ and $x = 5$ at $T = T_4$, respectively. We have seen that the only thermal resistance between $T_1$ and $T_2$ is the thermal conduction resistance, and its expression is composed of known constants:

$$\Delta T = T_1 - T_2 = \Phi_{tot} R_{cond,1 \to 2} = \Phi_{tot} \frac{e_{rod}}{\lambda_{rod}} \quad (30)$$

From which we can express $T_2$ with the only unknown $\Phi_{tot}$:

$$T_2 = T_1 - \Phi_{tot} \frac{e_{rod}}{\lambda_{rod}} \quad (31)$$

And similarly, we obtain

$$T_3 = T_4 + \Phi_{tot} \frac{e_{rod}}{\lambda_{rod}} \quad (32)$$

This allows us to express $R_{rad,2 \to 3}$ as a function of $\Phi_{tot}$:

$$R_{rad,2 \to 3} = \frac{\frac{1}{\epsilon_{rod}} + \frac{1}{\epsilon_{rod}} - 1}{\sigma(T_2 + T_3)(T_2^2 + T_3^2)} \quad (33)$$

$$R_{rad,2 \to 3} = \frac{\frac{1}{\epsilon_{rod}} + \frac{1}{\epsilon_{rod}} - 1}{\sigma\left((T_1 - \Phi_{tot}\frac{e_{rod}}{\lambda_{rod}}) + (T_4 + \Phi_{tot}\frac{e_{rod}}{\lambda_{rod}})\right)\left((T_1 - \Phi_{tot}\frac{e_{rod}}{\lambda_{rod}})^2 + (T_4 + \Phi_{tot}\frac{e_{rod}}{\lambda_{rod}})^2\right)} \quad (34)$$

We can now return to Equation 29 to find an equation with only one unknown, being $\Phi_{tot}$:

$$\frac{T_1 - T_4}{R_{cond,1 \to 2} + \frac{1}{\frac{1}{R_{cond,2 \to 3}} + \frac{1}{R_{rad,2 \to 3}}} + R_{cond,3 \to 4}} - \Phi_{tot} = 0 \quad (35)$$

$$\frac{T_1 - T_4}{\frac{e_{rod}}{\lambda_{rod}} + \frac{1}{\frac{1}{\frac{e_{gas}}{\lambda_{gas}} + \frac{1}{R_{rad,2 \to 3}(\Phi_{tot})}}} + \frac{e_{rod}}{\lambda_{rod}}} - \Phi_{tot} = 0 \quad (36)$$

We will not develop the equation further for readability reasons. However, Equation 35 can be numerically solved to find the value of $\Phi_{tot}$ setting it to 0. Once $\Phi_{tot}$ is determined, we can express the numerical values of $T_2$ and $T_3$. From that point, the temperature at any point in the system can be written as a piecewise
linear function, which is a simple linear interpolation between ($x_1 = 0, T_1$) and ($x_2 = 1, T_2$), ($x_2 = 1, T_2$) and ($x_3 = 4, T_3$), and ($x_3 = 4, T_3$) and ($x_4 = 5, T_4$) respectively.

$$ T(x) = \begin{cases} T_1 + (T_2 - T_1) \frac{x - x_1}{x_2 - x_1} & \text{if } x \in [x_1; x_2] \\ T_2 + (T_3 - T_2) \frac{x - x_2}{x_3 - x_2} & \text{if } x \in [x_2; x_3] \\ T_3 + (T_4 - T_3) \frac{x - x_3}{x_4 - x_3} & \text{if } x \in [x_3; x_4] \end{cases} \qquad (37) $$

This, at last, can be plotted against the computed value for the temperature taken at a vertical slice right in the middle of Figure 13. The results are shown Table 4. Since the simulation cannot really depict a 1D test case, the measured values vary inside an interval depending on where on the horizontal axis is the measure done. We take the range of values observed within the interval of width 20 as depicted on Figure 9, and we indicate it in the table.

<table><caption>Table 4: Comparison of temperatures with analytical solution for infinitely long rods conduction–radiation test case</caption><thead><tr><th>Source</th><th><i>T</i><sub>2</sub>(<i>K</i>)</th><th><i>T</i><sub>3</sub>(<i>K</i>)</th></tr></thead><tbody><tr><td>Analytical result</td><td>729</td><td>271</td></tr><tr><td>Present work</td><td>729-732</td><td>267-271</td></tr><tr><td>Relative difference (%)</td><td>0.4</td><td>1.4</td></tr></tbody></table>

The results are in very good accordance with the analytical predictions.

We can also plot the temperature variation along the width of the rod, both with analytical values and with our computed solution. This way we can assess the accuracy of the solution relatively to the position inside the rod. Figure 15 shows the temperature value inside the rod from bottom to top, depending on the distance $x$ from the bottom only since the problem is 1D. This figure also displays the absolute temperature error between the computed solution and the analytical one. The worst error is of around 6 K at width $x = 3.90$, right before the interface between air and upper steel part. This errors amounts to 2% of the analytical temperature at this point, which remains a fairly good accuracy.

### 4.2.2. 2D test cases: concentric cylinders

The next geometry we assess is that of two infinitely long, concentric hollow cylindrical rods. They can be described as two solid cylinders separated by a layer of gas. Since they are infinitely long, no physics is dependent upon the $z$ coordinate, meaning the physics that takes place can be described in 2D. The inner cylinder is comprised between radii $r_1 = 1.2$ and $r_2 = 3$, and the outer cylinder between radii $r_3 = 6$ and $r_4 = 7.5$. The geometry is constructed as follows. A global domain shaped like a 2D slice of a hollow cylinder comprised between radii $r_1$ and $r_4$ is defined and meshed uniformly. Essentially, the domain is comprised between two concentric circles. We then immerse into that domain a mesh of a 2D slice of a hollow cylinder comprised between radii $r_2$ and $r_3$. Once again this amounts to defining a domain comprised between two concentric circles. This immersion allows us to define three different spatial domains: between radii $r_1$ and $r_2$, between radii $r_2$ and $r_3$, and between radii $r_3$ and $r_4$. The mesh starts out from a uniform mesh of 8056 elements, and is then refined around the 0-isovalue of the levelset of the immersed cylinder slice down to 3966 elements. This amounts to refining the mesh around the $r_3$ and $r_4$ circles, since they form the implicit boundaries of the immersed slice. These circles are where radiation will take place. Figure 13 displays the setup.

We let thermal propagation run until steady state is reached after 50 timesteps of $2.10^4$ s each. Figure 14 shows the situation both for temperature and for radiative flux at different points of the simulation. As was the case for the rods, we can see a non-zero, but negligible, radiative flux at the two radiating circles at the beginning of the simulation.
<figure><img src="image_15.png" /><figcaption>Figure 12: Progression of the temperature inside the rod from bottom ($x=0$) to top ($x=5$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13: Immersion of a 2D cylinder slice in the domain. The two white circles stand for the 0-isovalue of the levelset defining the slice. The immersed item covers all the space between them. The mesh has been refined at the interface.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 14: Thermal propagation and radiative flux in the rods. The first line shows the initial setup. The second line shows the state after 10 timesteps. The third line shows the state after 50 timesteps, when steady state is reached.</figcaption></figure>

Even if the refinement does not seem very circle-shaped around the interfaces, the process of extracting virtual facets ensures the radiation is computed around a proper set of facets forming the correct circles. The values of emissivity, conductivity, specific heat and volumic mass are the same as in the infinite rods case. The inner side of the inner cylinder standing at $r_1$ is set to $T_1 = 1000$ K while the outer side of the outer cylinder standing at $r_4$ is set to $T_4 = 0$ K. Once again, conduction occurs through the thickness of the inner cylinder (standing for an arbitrary metal) between $r_1$ and $r_2$. Then, it occurs again between $r_2$ and $r_3$ in the thickness between the two cylinders (standing as an arbitrary radiation-transparent gas), plus radiation occurs between the outer surface of the inner cylinder at $r_2$ and the inner surface of the outer cylinder at $r_3$. Conduction occurs again in the thickness of the outer cylinder between $r_3$ and $r_4$. The physical parameters are $\epsilon = 0.5$ for all radiating surfaces, $c_{p,steel} = 25 \text{ W} \cdot \text{K}^{-1} \cdot \text{m}^{-1}$ for the conductivity in the cylinders and $c_{p,air} = 10 \text{ W} \cdot \text{K}^{-1} \cdot \text{m}^{-1}$ for the conductivity in the medium between the cylinders. We then let the system stabilize towards a steady state in temperature. We can devise an analytical solution for the thermal situation of the problem, exactly how we did so for the infinitely long rods. The only difference relies in the expression of the thermal resistances. Figure 10 still describes the electrical analogy of thermal transfer in this system. We take a set length of $1m$ for the cylinders, in order to make the length disappear from the equations. Thus, we have the following values for the various thermal conduction resistances:

$$R_{cond,1 \to 2} = \frac{\ln(\frac{r_2}{r_1})}{2\pi\lambda_{rod}} \qquad (38)$$

$$R_{cond,2 \to 3} = \frac{\ln(\frac{r_3}{r_2})}{2\pi\lambda_{gas}} \qquad (39)$$

$$R_{cond,3 \to 4} = \frac{\ln(\frac{r_4}{r_3})}{2\pi\lambda_{rod}} \qquad (40)$$
<figure><img src="image_18.png" /><figcaption>Figure 15: Progression of the temperature inside the rod from bottom (x=0) to top (x=5).</figcaption></figure>

And the value of the thermal radiation resistance is [50]:

$$R_{rad,2 \to 3} = \frac{\frac{1}{\epsilon_{rod}} \frac{r_3}{r_2} + \frac{1}{\epsilon_{rod}} - 1}{2\pi r_3 \sigma (T_2 + T_3)(T_2^2 + T_3^2)} \quad (41)$$

These changes done, we can once again compare the analytical solution to the computed solution at the radiating interfaces. The results are indicated on Table 5 and show the computed solution matches very well with the analytical solution.

<table><caption>Table 5: Comparison of temperatures and surface radiative flux with analytical solution for concentric cylinders conduction-radiation test case</caption><thead><tr><th>Source</th><th><i>T</i><sub>2</sub>(K)</th><th><i>T</i><sub>3</sub>(K)</th></tr></thead><tbody><tr><td>Analytical result</td><td>547</td><td>110</td></tr><tr><td>Present work</td><td>544-547</td><td>110-112</td></tr><tr><td>Relative difference (%)</td><td>0.7</td><td>1.6</td></tr></tbody></table>

We can also plot the temperature along with a radial progression through the cylinders, both for the analytical and computed solutions. This can be seen on Figure 15. Once again, we can see that our computed solution matches very well with the analytical solution.

### 4.2.3. 3D test cases: concentric spheres

Moving on to 3D geometries, we model conduction and radiation between four concentric spheres. The whole 3D domain is comprised between two concentric spheres $S_1$ and $S_4$ of respective radii $r_1 = 4$ and $r_4 = 10$. Two surface meshes, one of a sphere $S_2$ of radius $r_2 = 6$ and one of a sphere $S_3$ of radius $r_3 = 8$,
<figure><img src="image_19.png" /><figcaption>Figure 16: Immersion of two spheres in the spherical domain. The mesh has been refined around these two inner spheres. Conduction occurs in the bulk of the system, and radiation occurs between the two immersed spheres.</figcaption></figure>

are immersed in the domain. The space between $S_1$ and $S_2$, and $S_3$ and $S_4$ respectively, is filled with an arbitrary solid material, through which only conduction occurs. The space between $S_2$ and $S_3$ is filled with an arbitrary, immobile gas, transparent to radiation. This means both conduction through that gas and radiation between the two spheres occur at that stage. The setup is displayed on Figure 16. Starting from a uniform mesh with a mesh size of $h_{init} = 0.3$ and 635377 elements, we refine around the immersed inner spheres so as to achieve a mesh size around the interfaces $h_{final} = 0.03$ and 286434 elements. The outer sphere surface is set to $T_1 = 600$ K while the inner sphere surface is set to $T_4 = 0$ K. At the beginning of the simulation, temperature is equal to 0 K except for the outer sphere. We let thermal transfer heat up the volume starting from the outer sphere towards the inner sphere and wait for thermal equilibrium to be reached.

The equations are fairly similar to those of the 2D concentric cylinders test case. The expression of thermal conduction resistance $R_{cond,1→2}$ between two concentric spheres of respective radius $r_1$ and $r_2$ (with $r_1 < r_2$) separated by a solid medium of constant and homogeneous thermal conductivity $\lambda$ amounts to:

$$R_{cond,1→2} = \frac{1}{4\pi\lambda}\left(\frac{1}{r_1} - \frac{1}{r_2}\right) \qquad (42)$$

The thermal radiation resistance between these same spheres, of respective surface emissivity $\epsilon_1$ and $\epsilon_2$, is close to the expression for concentric circles, Equation 41:

$$R_{rad,1→2} = \frac{\frac{1}{\epsilon_1 r_2} + \frac{1}{\epsilon_2} - 1}{4\pi r_2^2 \sigma (T_1 + T_2)(T_1^2 + T_2^2)} \qquad (43)$$

These changes aside, the same computations can be made than for the 2D test cases, and similarly, we can compare the achieved temperature at certain specific points between the analytical solution and the computed solution. We look at the temperature at the surfaces of spheres 2 and 3. Table 6 shows the gathered data.
<figure><img src="image_20.png" /><figcaption>Figure 17: Progression of the temperature according to radius starting from the inner sphere surface at $r = 4$ to the outer sphere surface at $r = 10$</figcaption></figure>

<table><caption>Table 6: Comparison of temperatures and surface radiative flux with analytical solution for concentric spheres conduction–radiation test case</caption><thead><tr><td>Source</td><td>$T_2(K)$</td><td>$T_3(K)$</td></tr></thead><tbody><tr><td>Analytical result</td><td>355</td><td>491</td></tr><tr><td>Present work</td><td>338-344</td><td>490-493</td></tr><tr><td>Maximum relative difference (%)</td><td>-4.8</td><td>0.5</td></tr></tbody></table>

We can see that the results are in good agreement with the analytical values, though slightly less than for the 2D cases. And once again we can plot the temperature at a fixed radius and compare it to the analytical temperature. In theory, temperature should depend on radius only. The experimental data we use for comparison must not depend too much on local differences and mesh imperfections. This means we can not just draw an arbitrary radius segment through the sphere, gather the temperature along that segment, and plot it against the analytical value. Instead, we gather the temperature along with six different radius segments, and we then take the mean of all these to obtain a representative result from our computed data. For two of these six segments, we choose the following segments that run along the *x* axis: from ($r_1, 0, 0$) to ($r_4, 0, 0$), and from ($-r_4, 0, 0$) to ($-r_1, 0, 0$). We proceed similarly to obtain four other segments from the *y* and *z* axes. At last, the data is plotted against analytical value and shown Figure 17.

Once again we can see the computed temperature matches very closely with the analytical temperature. This test case gives us confidence in our ability to model 3D geometries with reasonable accuracy.

## 5. CONCLUSIONS

An immersed volume approach has been applied to conduction–radiation simulation cases. Radiating facets have been precisely extracted from the implicit nature of the boundaries between the different elements
in the mesh. This paper has provided a very general method of determining participating facets in the mesh,
and of handling obstruction while maintaining energy balance. A major point is the integration of a surface
phenomenon, radiation, into the resolution of a volumetric conduction equation, not only for the border
elements but also for the bulk elements of the mesh. The results have proven to be in very good accordance
with analytical solutions. An interesting future development would be try and apply this method to moving
objects inside an immersed mesh, for example ingots circulating through a furnace. The immersed volume
method lends itself very well to modelling of such moving objects, since it does not call for a remesh after
each movement of the objects. Another venue of progress could be to enhance the computation time of the
solver. This will be addressed in a future paper.

# 6. ACKNOWLEDGEMENTS

The authors would like to express their gratefulness to Thales Alenia Space France for having funded
the Ph.D thesis the results of which are presented in this paper.

## References

[1] F.-T. Lentes, N. Siedow, Three-dimensional radiative heat transfert in glass cooling processes, Glass Science and Technology : Glastechnische Berichte 335 (1999) 188-196.

[2] G. Thömmes, R. Pinnau, M. Seaïd, T. Götz, A. Klar, Numerical methods and optimal cotro for glass cooling processes, Transport Theory and Statistical Physics 31 (2002) 513-529.

[3] D. Clever, J. Lang, Optimal control of radiative heat transfer in glass cooling with restrictions on the temperature gradient, Optimal Control: Applications and Methods 33 (2012) 157-175.

[4] G. Krishnamoorthy, A computationally efficient p1 radiation model for modern combustion systems utilizing pre-conditioned conjugate gradient methods, Applied Thermal Engineering (june 2017) 197-206.

[5] L. Wu, X. Xu, W. Zhou, Y. Su, X. Li, Heat transfer analysis of blast furnace stave, International Journal of Heat and Mass Transfer 51 (jun. 2008) 2824-2833.

[6] V. Cuervo-Piñera, D. Cifrián-Riesgro, P.-D. Nguyen, V. Battaglia, M. Fantuzzi, A. D. Rocca, M. Ageno, A. Rensgard, C. Wang, J. Niska, T. Ekman, C. Rein, W. Adler, Blast furnace gas based combustion systems in steel reheating furnaces, Energy Procedia 120 (august 2017) 357-364.

[7] B. Thomas, I. Samarasekera, J. Brimacombe, Mathematical model of the thermal processing of steel ingots : Part i. heat flow model, Metallurgical and Material Transactions B (mar. 1987) 119-130.

[8] J. Wei, H. Zhang, L. Zheng, C. Wang, B. Zhao, Modeling and improvement of silicon ingot directional solidification for industrial production systems, Solar Energy Materials and Solar Cells 93 (september 2009) 1531-1539.

[9] Z. Malinowski, T. Telejko, B. Hadała, Influence of heat transfer boundary conditions on the temperature field of the continuous casting ingot, Archives of Metallurgy and Materials 57 (march 2012).

[10] S. L. Foll, Modelisation du couplage conduction/rayonnement dans les systèmes de protection thermique soumis a de très hauts niveaux de températures, Ph.D. thesis, INSTITUT NATIONAL DES SCIENCES APPLIQUEES DE LYON, 2014.

[11] L. A. Dombrovsky, S. Dembele, J.-X. Wen, A simplified model for the shielding of fire thermal radiation by water mists, International Journal of Heat and Mass Transfer 96 (may 2016) 199-209.

[12] A. V. Nenarokomov, L. A. Dombrovsky, I. V. Krainova, O. M. Alifanov, S. A. Budnik, Identification of radiative heat transfer parameters in multilayer thermal insulation of spacecraft, International Journal of Numerical Methods for Heat and Fluid 27 (march 2017).

[13] S. A. Khairallah, A. T. Anderson, A. Rubenchik, W. E. King, Laser powder-bed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spattern and denuation zones, Acta Materialia 108 (april 2016) 36-45.

[14] Z. Yan, W. Liu, Z. Tang, X. Liu, N. Zhang, M. Li, H. Zhang, Review on thermal analysis in laser-based additive manufacturing, Optics and Laser Technology 106 (october 2018) 427-441.

[15] P. Michaleris, Modeling metal deposition in heat transfer analyses of additive manufacturing processes, Finite Elements in Analysis and Design 86 (september 2014) 51-60.

[16] K. Reddy, N. S. Khumar, Combined laminar natural convection and surface radiation heat transfer in a modified cavity receiver of solar parabolic dish, International Journal of Thermal Sciences 109 (jan. 2008) 1647-1657.

[17] J. Polo, J. Ballestrín, E. Carra, Sensitivity study for modelling atmospheric attenuation of solar radiation with radiative transfer models and the impact in solar tower plant production, Solar Energy 134 (september 2016) 219-227.

[18] A. Razagui, K. Abdeladim, S. Semaoui, A. H. Arab, S. Boulahchiche, Modeling the forecasted power of a photovoltaic generator using numerical weather prediction and radiative transfer models coupled with a behavioral electrical model, energy Reports 6 (february 2020) 57-62.

[19] P. Cheng, Two-dimensional radiating gas flow by a moment method, American Institute of Aeronautics and Astronautics (1964) 1662-1664.

[20] R. McClarren, C. Hauck, Robust and accurate filtered spherical harmonics expansions for radiative transfer, Journal of Computational Physics (aug. 2010) 5597-5614.
[21] M. Malek, N. Izem, M. S. Mohamed, M. Seaid, M. Wakrim, Numerical solution of rosseland model for transient thermal radiation in non-grey optically thick media using enriched basis functions, Mathematics and Computers in Simulation 180 (february 2021) 258-275.

[22] K. Prasad, K. Vajravelu, P. Datti, B. Raju, Mhd flow and heat transfer in a power-law liquid film at a porous surface in the presence of thermal radiation, Journal of Applied Fluid Mechanics 6 (2013) 385-395.

[23] W. Fiveland, Three-dimensional radiative heat-transfert solutions by the discrete-ordinates method, Journal of Thermo-physics and Heat Transfer (oct. 1988) 309-316.

[24] P. Ricca, R. Furfaro, Analytical discrete ordinate method for radiative transfer in dense vegetation canopies, Journal of Quantitative Spectroscopy and Radiative Transfer (mar. 2013) 60-69.

[25] S. M. et al, Development and comparison of the dtm, the dom and the fvm formulations for the short-pulse laser transport through a participating medium, International Journal of Heat and Mass Transfer (jan. 2006) 1820-1832.

[26] F. Lockwood, N. Shah, A new radiation solution method for incorporation in general combustion prediction procedures, in: Eighteenth International Symposium on Combustion, 1981, pp. 1405-1414.

[27] S. H. Sarvari, A new approach to solve the radiative transfer equation in plane-parallel semitransparent media with variable refractive index basde on the discrete transfer method, International Communications in Heat and Mass Transfer (nov. 2016) 54-59.

[28] J. Farmer, J. Howell, Comparison of monte carlo strategies for radiative transfer in participating media, Advances in Heat Transfer (jun. 1998) 333-429.

[29] A. Kersch, W. Morokoff, A. Schuster, Radiative heat transfer with quasi-monte carlo methods, Transport Theory and Statistical Physics (1994) 1001-1021.

[30] S. Gratyi, D. Levin, A. Walker, Rassvet: Backward monte carlo radiative transfer in spherical-shell planetary atmospheres, Icarus (mar. 2010) 188-196.

[31] B. Gaume, Reduction d'un probleme d'auto-rayonnement par modes de Branche: Application aux echanges thermiques dans un domaine multi-enceintes, Ph.D. thesis, Universite Evry-Val d'Essonne, 2016.

[32] Q. Schmid, Numerical radiative transfer using an immerse volume method, Ph.D. thesis, Ecole Nationale Superieure des Mines de Paris, 2016.

[33] A. F. et al, Modeling of the surface-to-surface radiation exchange using a monte carlo method, Journal of Physics: Conference Series (2016).

[34] S. Bindick, M. Stiebler, M. Krafczyk, Fast kd-tree based hierarchical radiosity for radiative heat transport problems, Internation Journal for Numerical Methods in Engineering 86 (september 2014) 1082-1100.

[35] C. Sevart, T. Bergman, Evolutionay design method for a conducting solid cooled by combined free convection and radiation, Journal of Heat Transfer (january 2021).

[36] T. G. Revanna, S. S. Katte, Effects of axial conduction and internal radiation on thermodynamic optimization of a tubular space radiator, Thermal Science and Engineering Progress 22 (may 2021).

[37] K. De, H. Wu, G. Huang, X. Xu, Y. Liu, Condensation-free radiant cooling with double-skin infrared-transparent membranes, Building and Environment 193 (april 2021).

[38] E. M. Sparrow, Advances in heat transfer, volume 2, Elsevier, 1965.

[39] P. Hanrahan, D. Salzman, L. Aupperle, A rapid hierarchical radiosity algorithm, Computer Graphics 25 (1991).

[40] C. M. Goral, K. E. Torrance, D. P. Greenberg, B. Battaile, Modeling the interaction of light between diffuse surfaces, in: SIGGRAPH '84: Proceedings of the 11th annual conference on Computer graphics and interactive techniques, 1984, pp. 213-222.

[41] D. R. Baum, H. E. Rushmeire, J. Winget, Improving radiosity solutions through the use of analytically determined form-factors, in: Proceedings of the 16st Annual Conference on Computer Graphics and Interactive Techniques, 1989, pp. 325-334.

[42] E. Eckert, Radiative transfer, hc hottel and af sarofim, mcgraw-hill book company, new york, 1967. 52 pages, AIChE Journal 15 (1969) 794-796.

[43] G. Walton, Calculation of obstructed view factors by adaptive integration, National Institute of Standard and Technology Internal Report (nov. 2002).

[44] M. Cohen, M. F. Greenberg, The hemi-cube: a radiosity solution for complex environments, in: SIGGRAPH '85 conference proceedings, volume 19, July 1985, pp. 31-39.

[45] A. Franck, W. Heidemann, K. Spindler, Modeling of the surface-to-surface radiation exchange using a monte carlo method, Journal of Physics: Conference Series 745 (2016).

[46] C. B. S. Woop, A.T. Afra, Stbvh: A spatial-temporal bvh for efficient multi-segment motion blur, Proceedings of High Performance Graphics (jun. 2017).

[47] E. Hachem, S. Feghali, R. Codina, T. Coupez, Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation, International Journal for Numerical Methods in Engineering 94 (2013) 805-825.

[48] C. Bahbah, M. Khalloufi, A. Larcher, Y. Mesri, T. Coupez, R. Valette, E. Hachem, Conservative and adaptive level-set method for the simulation of two-fluid flows, Computers and Fluids 191 (2019).

[49] B. Py, Rayonnement thermique des métaux microrugueux ou dispersés, International Journal of Heat and Mass Transfer 10 (1967) 1735-1736.

[50] D. Marchio, P. Reboux, Introduction aux transferts thermiques, Presses des Mines, 2008.

[51] M. A. Ordal, L. L. Long, R. J. Bell, S. E. Bell, R. R. Bell, R. W. A. Jr, C. A. Ward, Optical properties of the metals al, co, cu, au, fe, pb, ni, pd, pt, ag, ti, and w in the infrared and far infrared, APPLIED OPTICS 22 (1983) 1099-1120.

[52] V. Prokopets, I. Shaikevitch, L. Robur, Optical properties of silicon-nitride-based ceramics with molybdenum and aluminum admixtures, Ukrainian Journal of Physics 50 (2005).
[53] C. Liebert, Emittance and absorptance of nasa ceramic thermal barrier coating system, NASA technical paper 1190 (1978).

[54] D. Law, D. Zhou, Developing Solid Oral Dosage Forms (Second Edition) - Chapter 3 - Solid-State Characterization and Techniques, Pharmaceutical Theory and Practice, 2017.

[55] H. Ebrahimi, A. Zamaniyan, J. S. Mohammadzadeh, A. A. Khalili, Zonal modeling of radiative heat transfer in industrial furnaces using simplified model for exchange area calculation, Applied Mathematical Modelling 37 (2013) 8004–8015.

[56] M. Polyanskiy, Refractive index.info, refractive index database, https://refractiveindex.info/, 2021.

[57] T. E. Ariceta, Infrared spectral emissivity studies on metals and materials for solar thermal applications, Ph.D. thesis, Euskal Herriko Uniberstitatea, 2016.

[58] P. Schroder, P. Hanrahan, On the form factor between two polygons, ACM SIGGRAPH Computer Graphics (sept. 1993) 163–164.

[59] M. Coehn, D. Greenberg, The hemi-cube, a radiosity solution for complex environments, ACM SIGGRAPH Computer Graphics (jul. 1985).

[60] S. Mazumder, Methods to accelerate ray tracing in the monte-carlo method for surface-to-surface radiation transport, Journal of Heat Transfer (feb. 2006) 945–952.

[61] A. Ambirajan, S. Venkateshan, Accurate determination of diffuse view-factors between planar surfaces, International Journal of Heat and Mass Transfer 36 (1993) 2203–2208.

[62] J. Howell, A catalog of radiation heat transfer configuration factors, http://www.thermalradiation.net/tablecon.html (1982).