Atin Kumar

Jérôme Vicente

Jean-Vincent Daurelle

Yann Favennec

Benoit Rousseau
email: benoit.rousseau@univ-nantes.fr

A modified zonal method to solve coupled conduction-radiation physics within highly porous large scale digitized cellular porous materials

Keywords: conductive-radiative heat transfer, coupling, modified zonal method, parallelization method

Due to their distinct textures, porous cellular materials have been of interest to many engineering applications. Energy conversion through such materials, especially at high temperatures, is governed by naturally coupled conduction-radiation physics. Therefore, a thorough understanding of the conduction-radiation behavior within these materials is important to properly design and optimize these materials. Several numerical methods developed over the last few decades allow to solve and study the influence of different textural parameters on coupled conductive-radiative heat transfers. These numerical methods require 3D digitized images of the sample of interest obtained either using 3D imagery technique or generated using numerical algorithms. To better represent the complex geometry, the 3D image of the sample requires high spatial resolution, and for a sample which is representative of involved physics might contain hundreds of millions of voxels which is complex to solve and computationally expensive.

To cope with this issue, we developed a parallelized discrete scale numerical approach using cell centered Finite Volume Method (FVM) and deterministic ray tracing to solve coupled heat transfer within highly porous large scale complex cellular materials. At the heart of this method rests a decomposition approach based on modified zonal method which significantly reduces the coupling efforts, memory requirements, and computation time. The results of two test cases presented in this study are cross-verified with those presented in literature and computed using Star-CCM+ software. Finally, the method is applied to virtual fibrous sample. These results present the ability of the solver to consider different boundary conditions, large temperature gradient across boundaries, and dealing with large samples consisting hundreds of million of voxels while providing accurate results. We also present and discuss the sensitivity of other associated parameters on the final results.

Introduction

Cellular materials [START_REF] Scheffler | Cellular ceramics: structure, manufacturing, properties and applications[END_REF] such as open-cell foams, manufactured lattices, fibers felts exhibit interesting properties such as high strength-to-weight ratio, high surface-to-volume ratio, good flow-mixing capacity which make them suitable for several engineering applications. Certain high temperature applications where cellular porous materials are used include solar power power plants [START_REF] Guévelou | Representative elementary volumes required to characterize the normal spectral emittance of silicon carbide foams used as volumetric solar absorbers[END_REF][START_REF] Wu | Experimental investigations of porous materials in high temperature thermal energy storage systems[END_REF][START_REF] Fend | Porous materials as open volumetric solar receivers: Experimental determination of thermophysical and heat transfer properties[END_REF], heterogeneous combustion [START_REF] Terracciano | Thermal and acoustic performance of al2o3, mgo-zro2, and sic porous media in a flow-stabilized heterogeneous combustor[END_REF], thermal protection systems [START_REF] Helber | Microstructure and gas-surface interaction studies of a low-density carbon-bonded carbon fiber composite in atmospheric entry plasmas[END_REF][START_REF] Borner | High temperature permeability of fibrous materials using direct simulation Monte Carlo[END_REF][START_REF] Banerji | Tomography-based radiative characterisation of decomposing carbonaceous heat shield materials[END_REF][START_REF] Paglia | Carbon-phenolic ablative materials for re-entry space vehicles: plasma wind tunnel test and finite element modeling[END_REF], etc. Note that, for porous samples having mean characteristic size of their constituents greater than a few microns or much greater than the size of thermal radiation can also be refereed as macroporous materials [START_REF] Taine | Statistical modelling of radiative transfer within non-beerian effective phases of macroporous media[END_REF][START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF].

From a thermal viewpoint, the effective heat transfer within macroporous materials which operate at high temperatures (i.e., >800 K) may involve all three modes of heat transfer namely conduction, convection, and radiation. Such heat transfer and its key quantity Effective Thermal Conductivity (ETC) is generally influenced by textural parameters such as shape, size, orientation, and spatial distribution of constituents [START_REF] Rousseau | Textural parameters influencing the radiative properties of a semitransparent porous media[END_REF], the porosity of the material [START_REF] Sans | Solving transient coupled conductive and radiative transfers in porous media with a monte carlo method: Characterization of thermal conductivity of foams using a numerical flash method[END_REF] as well as on the operational conditions [START_REF] Mora-Monteros | Effective conductivity of porous ceramics in a radiative environment[END_REF]. Moreover, at high temperatures, different modes of heat transfer are generally coupled i.e., heat transfer by one mode might be affected by other mode making their effective contribution non-linear [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF][START_REF] Badri | Conductive-radiative heat transfer within sic-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis[END_REF]. Due to this coupling effect, accurate estimation of all three modes of heat transfer, both individually as well as coupled, within the volume is still a topic of discussion within academic and industrial research teams. This knowledge is crucial for designing and optimizing heat transfer through such complicated porous materials. In this article, only macroporous materials having applications in vacuum will be considered i.e., convective heat transfer will be neglected, therefore, highlighting coupling effects between conduction and radiation only as in [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF][START_REF] Badri | Conductive-radiative heat transfer within sic-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis[END_REF].

Generally, designing macroporous media with an optimized ETC involves following steps: material elaboration, testing, and optimization [START_REF] Guévelou | Representative elementary volumes required to characterize the normal spectral emittance of silicon carbide foams used as volumetric solar absorbers[END_REF][START_REF] Heisel | Digital design and 3d printing of innovative sic architectures for high temperature volumetric solar receivers[END_REF]. Till date, several methods have been developed to determine coupled conductive-radiative heat transfer while developing the innovative porous materials. The classical approach which is based on experimental methods [START_REF] Coquard | Experimental investigations of the coupled conductive and radiative heat transfer in metallic/ceramic foams[END_REF] consists of several steps (i.e., material elaboration, estimation of required properties either through experimental methods or using a relevant inverse methods) and can be tedious because it is time consuming and economically expensive. This is often not compatible with the industrial mass production of a range of materials (having volume approximately a few cubic meters). Also, while manufacturing, it is generally challenging to tune a particular parameter affecting the heat transfer. Therefore, it is extremely difficult to determine the impact of a particular textural parameter on coupled heat transfer and in a competitive global context, this traditional approach turns out to be economically and temporally inefficient.

An alternative to experimental approach is to develop a numerical approach based on modeling of heat transfer physics using a numerical method. In this approach it is possible to explore the effect of a set of textural parameters on coupled conductive-radiative heat transfer within a digitized porous sample. This approach became popular with the advancement of X-ray micro-tomography [START_REF] Baruchel | X-ray tomography in material science[END_REF]. X-ray micro-tomography provides a stack of binary images of a sample which is then arranged together to generate a 3D (threedimensional) image of the sample. The resulting 3D image can further be used for detailed investigation of textural properties of the macroporous sample [START_REF] Brun | Geometrical measurement of real foams from 3d images[END_REF] and to perform numerical simulation on it [START_REF] Vignoles | Numerical study of effective heat conductivities of foams by coupled conduction and radiation[END_REF]. Furthermore, an in-depth analysis of spatial organisation of the sample allows to subsequently develop specific algorithms for generating virtual 3D structures with prescribed textural features [START_REF] Guévelou | Representative elementary volumes required to characterize the normal spectral emittance of silicon carbide foams used as volumetric solar absorbers[END_REF] which are more economic and less time consuming. Note that, the basic structural unit of 3D digitized image is a cubic cell called voxel and its size depends on the size of constituents to better capture the geometric details. Therefore, the final 3D image required for coupled heat transfer simulation which should be representative of the entire volume, also referred to as Representative Elementary Volume (REV), might consist of hundreds of millions of voxels [START_REF] Panerai | Micro-tomography based analysis of thermal conductivity, diffusivity and oxidation behavior of rigid and flexible fibrous insulators[END_REF].

From thermal modeling point of view, the approaches to solve coupled conductive-radiative heat transfer can be broadly classified into two categories: continuous scale approaches and discrete scale approaches. Continuous scale approaches consider the computational domain homogeneous and require homogenized properties such as effective solid thermal conductivity, volumetric extinction coefficient, scattering phase function, etc. to solve the heat transfer physics. Therefore, continuous scale approaches are also referred to as homogenized approaches. In contrary to continuous scale approach, discrete scale approach uses the actual geometry of the sample and requires local material properties such as thermal conductivity and optical properties of constituents. The discrete scale approaches allow to model heat transfer differently in different phases e.g., radiation in void phase and conduction in solid phase, and are preferred over continuous scale approach when local information is needed.

Furthermore, discrete scale approaches can be further classified into two subcategories: deterministic and stochastic. Using deterministic approach, Talukdar et al. proposed a methodology based on FVM integrated with block-off region approach to model conduction and radiation heat transfer within honeycomb structures [START_REF] Talukdar | Modelling of conduction-radiation in a porous medium with blocked-off region approach[END_REF]. Luo et al. developed a method combining a discrete ordinate ray tracing technique to simulate radiation transfer and a FVM scheme for thermal conduction heat transfer to compute the effective thermal conductivity of porous samples [START_REF] Luo | Thermal radiation effect on the effective thermal conductivity of porous material by pore level numerical simulation[END_REF][START_REF] Luo | Characteristics of effective thermal conductivity of porous materials considering thermal radiation: A pore-level analysis[END_REF]. Recently, some of the authors proposed a parallelized vectorial Finite Element Method (FEM) coupled with discrete ordinate method to model coupled conduction-radiation heat transfers within virtual SiC-based cellular foams [START_REF] Badri | Conductive-radiative heat transfer within sic-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis[END_REF] and within squared lattice open-cell foam [START_REF] Ouchtout | Finite element framework for modeling conducto-radiative transfers within heterogeneous media at both discrete and continuous scales[END_REF]. In their work, the domain decomposition approach was an extension of the well-known Schwarz method [START_REF] Dolean | An introduction to domain decomposition methods: algorithms, theory, and parallel implementation[END_REF] which has been successfully extended to a finite elements radiation solver in the context of high performance computation [START_REF] Badri | Vectorial finite elements for solving the radiative transfer equation[END_REF][START_REF] Badri | High performance computation of radiative transfer equation using the finite element method[END_REF][START_REF] Badri | Preconditioned krylov subspace methods for solving radiative transfer problems with scattering and reflection[END_REF][START_REF] Favennec | Ad hoc angular discretization of the radiative transfer equation[END_REF][START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF]. In parallel to deterministic approaches, Vignoles introduced a stochastic approach based on Monte-Carlo method to model coupled conductive-radiative heat transfer within fibrous media. The solid conduction was modeled using Brownian walkers within opaque solid domain and radiative heat exchange was modeled using ray tracing approach [START_REF] Vignoles | A hybrid random walk method for the simulation of coupled conduction and linearized radiation transfer at local scale in porous media with opaque solid phases[END_REF][START_REF] Vignoles | Computation of the conducto-radiative effective heat conductivity of porous media defined by triply periodic minimal surfaces[END_REF]. This strategy was applied to a digitized fibrous sample [START_REF] Vignoles | A hybrid random walk method for the simulation of coupled conduction and linearized radiation transfer at local scale in porous media with opaque solid phases[END_REF] having approximately 1 × 10 6 voxels and then was extended to tomographied cellular ceramics with a maximum size of 125 × 10 6 voxels [START_REF] Vignoles | Numerical study of effective heat conductivities of foams by coupled conduction and radiation[END_REF]. Another Monte-Carlo strategy based on an integral formulation of the temperature [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single monte carlo algorithm[END_REF] was proposed by Caliot et al. [START_REF] Caliot | Combined conductiveradiative heat transfer analysis in complex geometry using the monte carlo method[END_REF] to estimate the temperature profiles of numerical porous samples having twelve stacked Kelvin cells. Recently, Sans et al. [START_REF] Sans | Solving transient coupled conductive and radiative transfers in porous media with a monte carlo method: Characterization of thermal conductivity of foams using a numerical flash method[END_REF] developed a Monte-Carlo method to solve transient coupled conductive-radiative heat transfer within porous foams. Fan et al. [START_REF] Fan | Numerical investigations of the coupled conductive-radiative heat transfer in alumina ceramics[END_REF] proposed a method based on FVM and mesh agglomeration Monte-Carlo ray tracing method to solve coupled conductive-radiative heat transfer within open-cell foam having upto 1.6 × 10 6 mesh elements. Perraudin et al. [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF] developed a hybrid numerical scheme based on Monte Carlo Ray Tracing and cellcentered FVM to determine the conduction-radiation coupling effects within a structured porous geometry.

Till date, several local scale approaches have been developed but the main drawback is the computing time which could be significantly high if a huge 3D image (e.g., 1000×1000×1000 number of voxels) needs to be simulated either due to large size of REV or high spatial resolution required to capture geometric details of complex sample [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF][START_REF] Gonneau | Modelling of heat transfer within heterogeneous media by brownian walkers[END_REF]. Above mentioned numerical approaches are generally applied to simple geometries having total number of unknowns in the order of few millions. However, for samples having fine details such as fibrous materials, the REV for coupled conductive-radiative problem might contain hundreds of millions of unknowns. Therefore, we require a reliable and robust numerical solver which must not only be able to solve coupled conductive-radiative heat transfer physics with high accuracy but must also be able to solve the physics on geometries having huge number of unknowns within realistic time frames.

To cope up with above mentioned issues, this article proposes a parallelized local scale numerical approach to solve coupled conductive-radiative heat transfer. Technically, the proposed approach is similar to the one developed by Perraudin et al. [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF] i.e., FVM for solid conduction and ray tracing for radiative exchange computation with two main difference. First and foremost is the way conduction and radiation are coupled. Perraudin et al. solved the radiative exchange between each triangle of solid-void interface. Knowing that the radiative exchange may take place between two surfaces placed far apart, the radiative exchange factor matrix is a square matrix of side N where N is the number of triangles in the domain. This matrix is huge and memory intensive. Also, they mentioned that the size increases quadratically with the size of the domain. Hence, if this matrix does not fit into the Random Access Memory (RAM), this might increase the computation time significantly as the required value of radiative exchange factor needs to be read from a hard drive. Secondly, they used Monte-Carlo ray tracing to compute radiative exchange factor. The process is quasi-stochastic because the point of emission of rays is chosen stochastically, however, the reflection and absorption is modeled in a deterministic manner. Since the choice of the point of emission is not made deterministically, their method may lead to a stochastic error and requires a high number of rays as well as a high number of surfaces from which the rays are emitted.

In this work, we use a decomposition approach based on modified zonal method to divide a large volume into smaller subvolumes and deterministic ray tracing to overcome the above challenges. This paper presents the detailed methodology, mathematical model, and its numerical modeling. We will also discuss the constraints while solving conduction-radiation physics using modified zonal method along with two test cases: structured porous geometry as in [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF] and a kelvin cell foam as in [START_REF] Badri | Conductive-radiative heat transfer within sic-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis[END_REF] to cross-verify the numerical approach. The results obtained using our method will be compared against those published in literature and obtained using the commercial software Star-CCM+. Finally, the method is applied to highly porous virtual fibrous sample.

Methodology

The decomposition approach used in this work divides a large volume into smaller subvolumes. Such an approach was previously used by Hottel and Cohen [START_REF] Hottel | Radiant heat exchange in a gas-filled enclosure: Allowance for nonuniformity of gas temperature[END_REF] to solve radiation heat transfer within an absorbing, emitting, nonscattering gray gas with constant absorption coefficient and is known as zonal method. The zonal method approximates the spatial domain by discretizing it into finite, isothermal subvolumes and surface area zones. Using the zonal method, energy balance is performed between two zones using precalculated direct or total exchange areas. Furthermore, Larsen and Howell [START_REF] Larsen | The exchange factor method: an alternative basis for zonal analysis of radiating enclosures[END_REF] modified this method to express the energy balance in terms of exchange factor which is a more physical and measurable quantity. The modified zonal method has been use repeatedly to solve radiative heat transfer [START_REF] Cui | A new inverse approach for the equivalent gray radiative property of a non-gray medium using a modified zonal method and the complex-variable-differentiation method[END_REF]. In this work, the method is further extended to solve coupled conduction-radiation physics. Using modified zonal method, the conductive-radiative heat transfer problem is solved in four steps as detailed below:

Step 1: The entire volume is divided into smaller parallelepipedic subvolumes of almost equal size as shown in Figure 1. For the sake of simplicity and better representation, Figure 1 shows only a few subvolumes and for each subvolume, only a few voxels. However, in reality, each volume may have hundreds (or thousands) of subvolumes and each subvolume may consist of millions of voxels and surfacic triangles representing the solid void interface. Furthermore, using modified zonal method, each subvolume in itself is now a small domain consisting of a structured grid of hexahedrons also referred as voxels. From now on, terms such as 'global level', 'subvolume level' and 'voxel level' will be used to define different levels of the geometry. Also, each subvolume consists of both solid and void voxels. For the sake of brevity, a few number of voxels are shown but in reality, each subvolume may contain millions of voxels.

Step 2: For each subvolume, the effective solid thermal conductivity tensor (λ) is computed using voxelized FVM. The method for the same is presented in section 4.2.

Step 3: The radiative exchange factor is computed between each triangle within the volume and is stored between each subvolume (rather than between each triangle). This leads to a radiative exchange factor matrix of side N where N is the number of subvolumes (rather than number of triangles as in [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF]). The method for the same is presented in section 4.3.

Step 4: Once the effective solid thermal conductivity tensor (λ) of each subvolume and radiative exchange factor (F) between each subvolume is computed, the coupling is performed at the level of subvolumes. Note that each subvolume is considered as isothermal similar to the case of modified zonal method [START_REF] Hottel | Radiant exchange in a gas-fired enclosure: allowance for nonuniformity of gas temperature[END_REF] as discussed in the beginning of this section. Furthermore, while dividing the entire volume into subvolumes, additional precaution is required to be sure that each subvolume is connected to the boundaries so that heat conduction can take place. Also for effective radiative heat transfer, the subvolumes should be small enough such that radiative exchange can take place among the subvolumes which would otherwise be absorbed within the same subvolume and will lead to the extinction of the radiative heat transfer. The effect of subdivision will be discussed in detail in results and discussion.

To solve the coupled conductive-radiative heat transfer within porous media using modified zonal method, the following assumptions are made:

• The geometry consists of only two phases: solid and void. Both phases are assumed to be at rest which means that the heat is transferred by either conduction or radiation [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF]. Similar assumption was made in several works as in [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF][START_REF] Vignoles | A hybrid random walk method for the simulation of coupled conduction and linearized radiation transfer at local scale in porous media with opaque solid phases[END_REF][START_REF] Dauvois | Modélisation du transfert thermique couplé conductif et radiatif au sein de milieux fibreux portés à haute température[END_REF].

• The solid phase is assumed to be optically smooth (i.e., only specular reflection can occur) and opaque (i.e., no radiative transfer within the solid). All the radiative energy is absorbed at the surface itself which is then transmitted through conduction within the solid.

• The void phase is assumed to be transparent i.e., no absorption of radiation takes place within void phase as well as no heat transfer through convection.

• Solid plates are added near the hot and cold boundaries which behave as emitting boundaries at a fixed temperature. Emission is allowed towards the medium only.

• Using modified zonal method, the entire volume is subdivided into parallelepipedic subvolumes of almost equal size. Each subvolume is considered as isothermal with constant properties.

• The optical properties used for the computation are assumed to be independent of temperature.

Governing equations and boundary conditions

At the level of subvolumes, the coupled conductive-radiative heat transfer can be estimated using the coupled Heat Conduction Equation (HCE) and can be written as

∇ • (q cond (x) + q rad (x)) = 0 ∀x ∈ Ω ⊂ R 3 , (1)
where q is the heat flux, subscripts cond and rad denote conduction and radiation, respectively, x = (x, y, z) is the spatial location, and Ω = (Ω s ∪ Ω v ) denotes the open bounded set with subscripts s and v denoting solid and void, respectively. The volume integration of equation (1) yields the following form:

V ∇ • q cond (x) dV + V ∇ • q rad (x) dV = 0 . (2)
On applying the divergence theorem, also known as Gauss's theorem or Ostrogradsky's theorem, which converts the volume integral into area integral, the equation (2) becomes

A q cond (x) • n dA + A q rad (x) • n dA = 0 , (3)
where n is the normal direction vector. In equation (3), the first term is the contribution of solid conduction. This term can also be explained as the heat flow rate through conduction across a surface of area 'A'. In continuous form, heat flow rate across a surface, through conduction in steady-state can be written as

Q cond (x) = A q cond (x) • n dA = - A λ(x)∇T (x) • n dA , ( 4
)
where λ is the orthotropic thermal conductivity tensor of each subvolume and T is temperature.

In equation (3), the second term is associated to the radiative heat transfer. The radiative heat flow rate is estimated by taking the emission and absorption coming from all the surrounding objects into account. While using the ray tracing approach, the directional dependent absorption term can be reformulated in terms of the radiative exchange factor. The net radiative heat flow rate, adapted from equation (8.14) in [START_REF] Modest | Radiative Heat Transfer[END_REF], can be represented as:

Q rad (x) = A q rad (x) • n dA = ϵ(x)n 2 v σ B A(x)T 4 (x) - A ϵ(x ′ )σ B T 4 (x ′ )dF dA ′ →dA dA ′ - V 4κ P (x ′′ )σ B T 4 (x ′′ )dF dV ′′ →dA dV ′′ , (5)
where ϵ is the emissivity, x and x ′ denote emitting and reflecting surfaces, respectively, n v is the refractive index of void phase, σ B is the Stefan-Boltzmann constant, A denotes the surface area, κ P (x ′′ ) is the local Planck's mean absorption coefficient of the medium at x ′′ , dF dA ′ →dA denotes the radiative exchange factor between surface element dA ′ and dA (where dA ′ = dA(x ′ )), dF dV ′′ →dA denotes the radiative exchange factor between volume element dV ′′ and surface element dA (where dV ′′ = dV (x ′′ )). In the absence of participating medium and considering only surface-to-surface exchange, equation ( 5) can be simplified to

Q rad (x) = ϵ(x)n 2 v σ B A(x)T 4 (x) - A ϵ(x ′ )σ B T 4 (x ′ )dF dA ′ →dA dA ′ . (6)
Here, F dA ′ →dA is the fraction of the total energy emitted by a surface dA ′ that is absorbed by dA, either directly or after any number of reflections. F dA ′ →dA is also called radiative exchange factor. Also, it is to be recalled here that the radiation travels within the void phase only and as soon as it is absorbed by the surface, the energy within the solid is transferred by solid conduction. Hence, equation ( 6) is solved only in void phase. To solve coupled conduction-radiation physics using above equations, following boundary conditions are required:

• Dirichlet boundary condition for equation ( 4) is prescribed at hot and cold boundaries as

T (x = 0) = T hot , T (x = L x ) = T cold , (7)
where x-direction is the direction of heat flow and L x is the length of domain in x-direction.

• Neumann boundary condition for equation ( 4) is prescribed on side faces to avoid any heat loss as

∇T • n(y = 0 & y = L y ) = 0 , ∇T • n(z = 0 & z = L z ) = 0 . (8)
• The emitting boundary condition is prescribed at fixed temperature boundaries (i.e., at hot and cold boundaries) which emit the radiation within the domain at fixed temperature.

• On four sides walls, either mirror (with ρ = 1, where ρ is the reflectivity) or vacuum (surface at T = 0 K also known as absorbing) boundary condition is prescribed. In case of mirror boundary condition, the radiation is reflected back within the domain conserving the total energy of the system. However, in case of vacuum boundary condition, the radiation that strikes any of the four side walls does not reenter, resulting in loss of total energy of the system.

• Also, the reflecting boundary condition is applied at the interface with reflectivity (ρ) = 1 -ϵ.

Numerical modeling

Solving coupled conductive-radiative heat transfer at the level of subvolumes

Using the FVM for a 2D structured grid (see Figure 2), equation ( 4) can be written in discretized form as

Q cond,(i) = (λA∇T ) W i -(λA∇T ) E i +(λA∇T ) S i -(λA∇T ) N i +(λA∇T ) F i -(λA∇T ) B i ∀ i = 1, 2, . . . , N , (9)
where E, W , S, N , F , and B denote the East, West, North, South, Front, and Back coordinates, respectively. The discretised equation has a clear physical interpretation. Equation [START_REF] Paglia | Carbon-phenolic ablative materials for re-entry space vehicles: plasma wind tunnel test and finite element modeling[END_REF] states that the Q cond entering from the west, south, and front face minus Q cond leaving the north, east, and back face are equal at steady state, i.e., it constitutes a balance equation for Q cond over the control volume. Now using a finite difference stencils formula, equation ( 9) can be written as where n denotes the number of iterations, and w denotes the equivalent thermal conductivity weight. Different conductivity weights can be calculated as

Q cond,(i) = w T i T n+1 i -w E T n E -w W T n W -w N T n N -w S T n S -w F T n F -w B T n B ∀ i = 1, 2, . . . , N , (10)
w E i = 2/ ∆x i λ i + ∆x E λ E = w E , w W i = 2/ ∆x i λ i + ∆x W λ W = w W , (11)
w N i = 2/ ∆y i λ i + ∆y N λ N = w N , w S i = 2/ ∆y i λ i + ∆y S λ S = w S , (12)
w F i = 2/ ∆y i λ i + ∆y F λ F = w F , w B i = 2/ ∆y i λ i + ∆y B λ B = w B , (13)
w T i = w E + w W + w N + w S +w F + w B . (14)
Here, λ denotes the thermal conductivity of the corresponding subvolume (in a particular direction denoted using indices 'E', 'W ', 'N ' 'S', 'F ', and 'B'), and ∆x, ∆y, and ∆z denote the size of the subvolume in x, y, and z-directions, respectively. Although, the above discretization is shown in two dimensions, an expansion to three dimensions is straight forward.

Furthermore, to solve equation ( 6) numerically, it is necessary to write it into discretized form. Assuming that the entire volume is divided into N number of subvolumes, the equation ( 6) will be transformed into the following equation

Q rad,(i) = ϵ i n 2 v σ B T 4 i A i - N j=1 j̸ =i ϵ j σ B T 4 j A j F j→i ∀ i = 1, 2, . . . , N , (15)
where A denotes the surface area of solid-void interface within a subvolume. It is worth mentioning that Q rad,(i) depends on the temperature of all surfaces but F j→i is either independent (in case of gray surfaces) or only depends on the temperature of emitting surfaces (in case of nongray surfaces). Furthermore, within a closed system, the energy remains conserved i.e., by the second law of thermodynamics, the net radiative exchange between two surfaces maintained at equal temperatures is equal to zero. Therefore, the reciprocity relation and summation rule also holds for radiative exchange factor. Note that the summation rule is true only in case of the closed enclosures without any opening. By applying the reciprocity and summation rule, equation ( 15) can be written as

Q rad,(i) = N j=1 j̸ =i ϵ i n 2 v σ B A i F i→j T 4 i -T 4 j ∀ i = 1, 2, . . . , N . (16)
The above equation is nonlinear and needs to be linearized to solve the system of equations using numerical method. The (T 4 i -T 4 j ) term is linearized by taking the factor as

T 4 i -T 4 j = (T 2 i + T 2 j )(T i + T j )(T i -T j ) ∀ i, j = 1, 2, . . . , N . (17)
Note that the linearization of T 3 equivalent term is not based on any assumption. It is just a mathematical simplification unlike other studies where [START_REF] Sans-Laurent | Modeling the flash method by using a conductoradiative monte-carlo method: Application to porous media[END_REF][START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams with a flash method at high temperature[END_REF]. This mathematical simplification shown in equation ( 17) allows to consider large temperature gradient between hot and cold boundaries and results of which will be presented in section 5. Furthermore, using equation ( 16) and ( 17), we get

T 4 i -T 4 j ≈ 4T 3 ref (T i -T j ), where T ref = (T i +T j )/2 is assumed
Q rad,(i) = N j=1 j̸ =i h n r,i→j F i→j T n+1 i -T n j ∀ i = 1, 2, . . . , N , (18)
where

h n r,i→j = ϵ i n 2 v σ B σA i T 2 i + T 2 j n (T i + T j )
n is a constant, superscript n denotes the number of iterations, and A denotes the total surface area of all triangles within a subvolume. While calculating h r,(i)→(j) , T i and T j are assumed to be known quantities. Generally, computation of F i→j is long and can be memory intensive if one records the exchange between each triangle [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF]. It would then be nearly impossible to consider a large mesh with hundreds of millions of triangles (or surfaces). Using modified zonal method, we calculate the radiative exchange among all the triangles within the volume but store the radiative exchange factor between subvolumes, as in case of modified zonal method [START_REF] Larsen | The exchange factor method: an alternative basis for zonal analysis of radiating enclosures[END_REF]. This will be discussed in detail in section 4.3.

Finally, the coupling can be performed over a simplified geometry having subvolumes equal to number of subdivisions. The modified zonal method makes the coupling simple and efficient. The discretized form of coupled HCE involves calculating T i ∈ Ω ⊂ R 3 which can be written as

w T i + h n r,i→i (1 -F i→i ) T n+1 i = w E T n E + w W T n W + w N T n N + w S T n S +w F T n F + w B T n B - N j=1 j̸ =i h n r,i→j F i→j T n j ∀ i = 1, 2, . . . , N . (19)
Equation ( 19) is a linear equation and explicit for T n+1 i written on the left hand side of the equation. Solving this equation at the level of subvolume requires the knowledge of equivalent thermal conductivity tensor (λ i ) of each subvolume (cf. section 4.2), and the radiative exchange factor (F i→j ) between subvolumes (cf. section 4.3).

Note that the coupled heat conduction equation is greatly resolved using modified zonal method and introducing radiative exchange factor for directional dependent term. Once the effective solid thermal conductivity tensor and radiative exchange factor are computed, the system can be assembled and coupling can be performed by solving the set of linear equations iteratively using any numerical method. A schematic showing the coupled system using resistance network and possible boundary conditions is shown in Figure 3.

But before going further, a special attention is required in case of coupling. Firstly, equation [START_REF] Brun | Geometrical measurement of real foams from 3d images[END_REF] written for each subvolume leads to a square matrix of N × N , where N denotes the number of subvolumes. But unlike the pure conduction problem, the resulting matrix is dense. Such a matrix requires additional care while solving numerically, otherwise, it may lead to an unexpected divergence of the system. In this work, thanks to a small size of matrix, a direct solver developed in an open source library Eigen [START_REF] Guennebaud | Eigen v3[END_REF] was used. The matrix structure in case of coupling is represented as

AX = B , (20)
where A = A C + A R with A C and A R representing the contribution form solid conduction and radiation, respectively, X is the unknown temperature vector, and B is the known boundary condition vector. These matrices and vectors are further expanded as

A C =     w t 1 w S 1 • • • 0 w N 1 w t 1 • • • 0 . . . . . . . . . . . . 0 • • • w N N w T N     , A R =     h r,(1→1) (1 -F 1→1 ) h r,(1→2) F 1→2 • • • h r,(1→N ) F 1→N h r,(2→1) F 2→1 h r,(2→2) (1 -F 2→2 ) • • • h r,(2→N ) F 2→N . . . . . . . . . . . . h r,(N →1) F N →1 • • • h r,(N →N -1) F N →N -1 h r,(N →N ) (1 -F N →N )     X =     T 1 T 1 . . . T N     , B = T hot     h r,(1→B1) F 1→B1 + h r,(1→B1) F 1→B1 h r,(2→B1) F 2→B1 + h r,(2→B1) F 2→B1 . . . w E N + h r,(N →B1) F N →B1 + h r,(N →B1) F N →B1     +T cold     w W 1 + h r,(1→B2) F 1→B2 + h r,(1→B2) F 1→B2 w W 1 + h r,(2→B2) F 2→B2 + h r,(2→B2) F 2→B2 . . . h r,(N →B2) F N →B2 + h r,(N →B2) F N →B2     ,
where subscript N denotes the number of subvolumes. Secondly, even though the coupled HCE is greatly reduced (see equation ( 19)), the dense matrix is solved and the solution is converged, the problem remains non-linear. This is due to the fact that the equivalent T 3 term is assumed to be constant which is calculated using the initial temperature or temperature obtained in previous iteration. Therefore, another non-linear loop needs to be setup. Overall, while solving the coupled conduction-radiation problem, temperature is initialized for each subvolume and constants (w and h) are calculated. Then, the coupled HCE is solved until convergence. The resulting temperature field is compared with the previous temperature field. If the variation is less than a predefined threshold (ε), the solution is accepted. Otherwise, the temperature field is updated and constants are recalculated to resolve the coupled HCE. This process is repeated until convergence is reached. Algorithm 1 presents the two loops iterative algorithm to solve the coupled conduction-radiation problem. T m = T n 14 end while

Computing effective thermal conductivity tensor of each subvolume

In order to determine the orthotropic thermal conductivity tensor of each subvolume, the heat conduction problem needs to be solved at voxel level. Note that, at the level of voxels, equation (4) persists. However, the heat conduction problem is to be solved within the solid domain only because the void phase is assumed to be vacuum i.e., λ v = 0. At the level of voxels, the heat flow rate can be written as

Q cond (x) = A q cond (x) • n dA = - A λ s (x)∇T (x) • n dA = 0 , (21)
where λ s is the solid phase orthotropic thermal conductivity tensor as in equation ( 4) except that equation ( 21) is solved for the solid voxels only, therefore, the subscript 's' is used. At voxel level, we assume that λ x = λ y = λ z = λ s so, the above equation can be re-written as

Q cond (x) = - A λ s ∇T (x) • n dA = 0 . (22)
Furthermore, before solving equation [START_REF] Talukdar | Modelling of conduction-radiation in a porous medium with blocked-off region approach[END_REF], it is necessary to find the connected volume i.e., the solid part should be connected to the boundaries on which the temperature gradient is applied. This is necessary because there is no heat transfer within hanging control volume, definition of initial temperature in those control volumes may lead to wrong estimation of flux. Hence, unconnected volumes need to be removed before computing the temperature field. Once the connected volumes are filtered, the discretized form of the above equation is solved to calculate T ∈ Ω s ⊂ R and can be written as

w T i T n+1 i = w E T n E + w W T n W + w N T n N + w S T n S +w F T n F + w B T n B ∀ i = 1, 2, . . . , N . (23)
The above equation is exactly the same as of equation [START_REF] Taine | Statistical modelling of radiative transfer within non-beerian effective phases of macroporous media[END_REF] except that the equation ( 23) is to be solved for solid voxels only. This leads to a seven diagonal sparse metrix which is then solved using Generalized Minimal Residual Method (GMRES) method developed in PARALUTION, an open source library [START_REF] Labs | Paralution v1.1.0[END_REF].

Once the temperature field is determined, the conductive heat flux across the boundaries (of known temperature) can be calculated using the well-known Fourier's law. For a case of one dimensional heat transfer along e x -direction, the effective conductive heat flux (q eff cond,x ) can be written as

q eff cond,x = q eff cond • e x = λ eff cond,x T hot -T cold L x , (24)
where e x = (1, 0, 0) T denotes the direction vector along x-direction, the subscript x denotes the direction of heat flow, T hot and T cold are the temperature at hot and cold boundaries, respectively, L x is the length of the domain along the direction of heat flow and λ eff cond is the effective solid conductivity of the material which is to be determined. Furthermore, conductive heat flux towards x-direction can be determined as

q cond,x (x) = q cond (x) • e x = -λ s ∇T (x) • e x = -λ s ∂T (x) ∂x . ( 25
)
Volume averaged conductive heat flux in x-direction q cond,x can be calculated using the following relation

q cond,x = 1 V Vs q cond,x (x) dV = - 1 V Vs λ s ∂T ∂x dV . (26)
Since the total conductive heat flux along the direction of heat flow remains conserved, from equation ( 24) and ( 26), we could deduce the directional effective thermal conductivity due to conduction only as

λ eff cond,x = L x T hot -T cold q cond,x = L x T hot -T cold 1 V Vs λ s ∂T ∂x dV . (27)
Similarly, λ eff cond is computed in all three directions to compute effective thermal conductivity tensor of each subvolume.

The computation of effective solid thermal conductivity tensor of each subvolume is independent to each other, therefore, parallelization was quite simple due to embarrassingly parallel problem. Once the entire volume is divided into subvolumes and each subvolume is checked for connectivity as discussed in previous section, each subvolume is then passed to different processors. Within each processor, the matrix is formed for respective subvolume and is stored in its memory. Once all the matrices (equal to the number of processors available for computation) are formed, PARALUTION is called. The PARALUTION solves the matrix using multiple processors available for computation using open MPI support. Then, the same step is repeated for remaining subvolumes and for other directions.

Computing radiative exchange factor between each subvolume

In order to compute the radiative exchange factor, ray tracing approach is used. Using this approach, a surfacic mesh which consists of triangles is generated using fast marching cube algorithms initially developed in iMorph software [START_REF] Brun | Imorph: A 3d morphological tool to fully analyse all kind of cellular materials[END_REF][START_REF] Rousseau | Evolution of the homogenized volumetric radiative properties of a family of α-SiC foams with growing nominal pore diameter[END_REF]. Rays are emitted from the center of each triangle towards the void phase in prescribed number of directions within the hemisphere. Each ray carries a fixed amount of energy which is normalized by the surface area of triangle, and rays are traced as they travel within the medium until they are fully absorbed or leave the system. A schematic of ray tracing is shown in Figure 4. Note that for conciseness, one triangle per subvolume is shown but as mentioned in section 2, each subvolume might consist millions of number of triangles.

As rays travel within the volume, at each interaction, the amount of energy reflected is computed. For this purpose, the value of reflectivity can be fixed (in case of gray surface) or can be computed using the Fresnel's relation from the knowledge of optical properties. Rest of the energy after reflection is considered to be absorbed by the surface. The energy balance is performed by assuming the amount of energy carried by an emitted ray from surface placed in subvolume (i) denoted by E e,(i) = A i N s,(i) /A total,(i) . Here A (i) denotes the area of the triangle and N s,(i) denotes the number of triangles within a subvolume i. A fraction of the energy will be reflected by the surface placed in subvolume j is denotes by E e,(j) = ρE e,(i) , and finally, the energy absorbed by surface placed in subvolume j denoted by E a,(j) = E e,(i) -E e,(j) . This is to be recalled that in this work, we only consider specular reflection and the value of reflectivity is either provided to the code or computed from the knowledge of optical properties (n -ik).

The ray tracing approach used in this work is based on a fast voxel traversal algorithm [START_REF] Amanatides | A fast voxel traversal algorithm for ray tracing[END_REF]. The method is simple and consist of two steps: Initialization and transversal increment. In initialization step, the emitting point of the ray and the step value in each direction (1 for increment and -1 for decrement) are initialized which can be denoted as ∆x = [∆x, ∆y, ∆z]. In the next step, the time t(x) at which the ray crosses the first boundary in x, y, and z-directions is recorded. The minimum of the three time steps (t) defines the minimum distance which can be traveled along the ray within the current voxel without leaving it. Finally, t × ∆x is computed which indicates the distance which can be traveled along a direction. Detailed algorithm of the method can be found in [START_REF] Amanatides | A fast voxel traversal algorithm for ray tracing[END_REF].

Then, the radiative exchange factor is computed between subvolumes. Each subvolume can consist of number of triangles and each triangle is allowed to emit several rays in all the directions within the hemisphere. The radiative exchange factor is calculated by recording the energy emitted by all the triangles within a subvolume and a fraction of which is absorbed by all the triangles in other subvolumes. The radiative exchange factor between subvolumes can be written as

F (i)→(j) =
total energy absorbed by all the surfaces within subvolume j either directly or after a number and types of reflections total energy emitted by all the surfaces of subvolume i .

In this work, the radiative exchange factor is computed between each subvolume and also between subvolumes and boundary walls. Unlike the solid conduction, where exchange with and only with neighboring subvolume is possible, the radiative exchange can take place with subvolumes placed far away. It is to be noted that the ray tracing procedure is completely separated from the computation of effective solid thermal conductivity tensor and the coupling can be performed separately.

The biggest advantage of ray tracing is that it is easily parallelizable because each ray is independent of others. This makes ray tracing method well suited for large scale complex geometries consisting of hundreds of millions of triangles and from each triangle, hundreds of rays are emitted which leads to manage billions of rays. Also, multiple reflections may further increase the complexity and hence, the computation time. The parallelization strategy which was initially developed in "genMat" for the Monte Carlo Ray Tracing (MCRT) computation is adapted in this work [START_REF] Guévelou | Representative elementary volumes required to characterize the normal spectral emittance of silicon carbide foams used as volumetric solar absorbers[END_REF][START_REF] Guévelou | A simple expression for the normal spectral emittance of opencell foams composed of optically thick and smooth struts[END_REF]. Rays are initialized by defining the emitting point, direction, and initial energy, and they are stored in a stack. Each ray is then accounted by a different processor and its movement is traced as it travels within the domain until it leaves the system or is fully absorbed. The process is repeated for all the rays. Finally, the energy corresponding to each triangle within the domain is accounted to compute the radiative exchange factor.

Numerical Experiments

In this section, we present and discuss the results of the two different test cases obtained using the numerical tool developed in this work. The first test case deals with the cross-verification of numerical results against literature results and the second test case presents the ability of the numerical approach to deal with different boundary conditions, huge temperature gradient over the boundaries, as well as effect of other input parameters on final results.

To describe and analyze the behavior of numerical experiments, the following quantities are used:

• Plane averaged temperature T x=x , i.e., is the mean temperature of any yz-plane located at x = x :

T x=x = Ω T (x)1 [x=x] dx Ω 1 [x=x] dx . (29)
• Similarly, plane averaged conductive heat flux (q cond ) in x-direction can be determined as

q cond,x=x = -λ s Ωs T + -T - ∆x 1 [x=x] dx Ω 1 [x=x] dx , (30)
where T -and T + are the temperatures on the left and right faces of the control volume, respectively, ∆x is the length of the control volume in x-direction.

• The plane averaged radiative flux (q rad,x=x (i)→(j) ) in a direction is calculated by taking the directional component of it. So, plane averaged radiative heat flux in x-direction can be calculated as

q rad,x=xi→j = σ B n 2 v ϵF i→j T 4 i -T 4 j • e x ∀ i, j = 1, 2, . . . , N , (31)
where e x = (1, 0, 0) T denotes the direction vector along x-direction.

• Temperature difference :

∆T = T hot -T cold , (32)
where T hot and T cold are the temperature at hot and cold boundaries, respectively.

• Reference temperature i.e., the mean of the temperatures applied at boundaries :

T ref = T hot + T cold 2 . ( 33
)
• Non-dimensional distance in x-direction :

X * = x/L x , (34)
where L x is the total length of the domain in x-direction.

• Effective coupled thermal conductivity:

λ eff coup = L x T hot -T cold q coup,x , (35)
where q coup,x = q cond,x + q rad,x denotes the volume averaged coupled heat flux. q cond,x is calculated using equation ( 26) and to calculate q rad,x , the radiative heat flux is calculated as

q rad,xi→j (x) = q rad,xi→j (x) • e x = σ B n 2 v ϵF i→j T 4 i -T 4 j • e x ∀ i, j = 1, 2, . . . , N . (36)
Then, volume averaged radiative heat flux in x-direction q rad,x can be calculated as

q rad,x = 1 V Vv q rad,x (x) dV . (37)

Test case 1: Cross-bar geometry

The objective of this test case is to cross-verify the numerical method developed in this work against the results published in literature. It must be noted that there are only a handful of numerical studies related to numerical modeling of coupled heat transfer which can be used for the cross-verification. As the first test case, we considered a simple geometry consisting of twelve square cross-sectional bars. The geometry is proposed by Perraudin et al. [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF] to study pure conduction, pure radiation, and conduction-radiation coupling. In this section, we detail and discuss the results in terms of effect of number of subdivisions as well as the effect of conduction-to-radiation parameter on coupled heat transfer.

Geometry and input parameters for test case 1

The geometry was reconstructed using C++ code to produce a stack of binary images and is shown in Figure 5 with applied boundary conditions. This geometry is composed of two phases namely solid (shown in black) and void (surrounded by solid phase) and both having equal volume fractions (0.5). To study coupling, two physical boundary walls were added at hot and cold sides also shown in Figure 5. The domain measures 1 × 1×1 cm 3 and is composed of 32 voxels in each direction giving total number o voxels ≈33 × 10 3 and total number of triangles ≈15 × 10 3 . As in reference [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF], we considered that the solid phase is opaque and optically smooth with solid phase thermal conductivity (λ s ) equal to 5 × 10 -3 W m -1 K -1 and void phase (which is replaced by a conducting, non-participating media) thermal conductivity (λ v ) equal to 5 × 10 -2 W m -1 K -1 giving λ v /λ s = 10. Note that Perraudin et al. studied λ v /λ s varying from 2.1 × 10 -4 -4800 to study the effect of solid and void phase thermal conductivity on coupled conductive-radiative heat transfer. However, plots of plane-averaged temperature is given for λ v /λ s = 10, only which are used for the cross-verification of this test case. For the simulation, 1314 number of rays per triangle were emitted and a total of 16 (3.5 GHz) processors were used.

Results of test case 1

Effect of number of divisions: using the modified zonal method, the coupling was performed for different number of subvolumes. The reference temperature was quite high in this case i.e., T ref = 2800 K with ∆T = 100 K. At such high temperature, it is expected that contribution of radiation can be significant. Therefore, before studying the effect of subdivisions on coupled heat transfer, it is necessary to check the effect of number of subdivisions on solid conduction only.

To study the effect of pure conduction, the entire geometry was divided into number of subvolumes and for each subvolume, the effective solid thermal conductivity tensor (λ) was calculated. Then, the HCE was solved again at the level of subvolumes and the plane-averaged temperature profile (T ) for different number of divisions is plotted in Figure 6. Note that the number of points at which the temperature is computed are equal to the number of divisions in the direction of heat flow (plus two more points, one at each boundary). The temperature profiles show a good agreement with the literature. In case of N div = 4 (which denotes the number of divisions in each direction), a slight mismatch can be seen between two temperature profiles. This is due to the resolution as four subvolumes are not enough to adapt the temperature gradient accurately. This mismatch reduces by increasing the number of divisions. It can be seen for N div = 8 and 16 that temperature profiles are in good agreement. Furthermore, it is worth mentioning that the λ eff cond was computed to be the same in all three cases i.e., 2 × 10 -2 W m -1 K -1 . Therefore, it is concluded that in this case, the divisions of geometry does not have any effect on solid conduction. Next, the effect of number of subdivisions on conduction-radiation coupling is performed and the results for different number of divisions is plotted in Figure 7.

N div = 4 N div = 8 N div = 16
Figure 7 shows that the results follow a similar trend and are in good agreement with literature. A slight difference in temperature profile can be noticed in case of coupling. At this stage, it is necessary to recall the methodologies used by Perraudin et al. [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF]. As discussed in section 1, they accounted the radiative exchange factor between each surface giving more localized control of heat transfer. However, in this work the radiative exchange factor is accounted between subvolumes which could be the reason of slight mismatch in the two results.

Furthermore, it can be seen that the average temperature increases near the cold boundary as high intensity radiation coming from the hot region will be absorbed here, and vice versa occurs near the hot boundary. This is a common phenomenon in case of semi-transparent media [START_REF] Viskanta | Boundary layer in thermal radiation absorbing and emitting media[END_REF]. Figure 8 shows the difference between plane-averaged temperature profiles of pure conduction (cf. Figure 6) and coupled heat transfer (cf. Figure 7). The figure confirms that the average temperature gain near the cold boundary (at X * = 0) is higher than the temperature loss near the hot boundary (at X * = 1) and also confirms that the profile is asymmetric as shown in various works [START_REF] Perraudin | Numerical quantification of coupling effects for radiation-conduction heat transfer in participating macroporous media: Investigation of a model geometry[END_REF][START_REF] Dauvois | Modélisation du transfert thermique couplé conductif et radiatif au sein de milieux fibreux portés à haute température[END_REF][START_REF] Viskanta | Boundary layer in thermal radiation absorbing and emitting media[END_REF]. Effect of Stark number or conduction-to-radiation parameter: the conduction-to-radiation parameter (N ) [START_REF] Viskanta | Boundary layer in thermal radiation absorbing and emitting media[END_REF] which is also known as Stark number as discussed in [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Siegel | Thermal radiation heat transfer[END_REF] is defined as

N = κλ eff cond 4σ B n 2 T 3 ref , ( 38
)
where κ is the effective volumetric absorption coefficient and n is the effective refractive index of the medium.

Figure 9 shows the plane-averaged temperature profile (T ) for different Stark numbers calculated by varying λ eff cond .

At high value of the Stark number (N = 0.1), the variation between temperature profiles for pure conduction and coupled heat transfer is negligible and it seems like both the profiles are overlapping. It means that the relative contribution of conduction is much higher than the radiative one. As the Stark number decreases, the relative contribution of radiation increases which is clearly reflecting in the temperature profile in Figure 9. As expected, with the increase in radiative contribution, the temperature increases near the cold boundary and decreases near the hot boundary. The effect is more significant in case of low value of Stark number. Viskanta et al. [START_REF] Viskanta | Boundary layer in thermal radiation absorbing and emitting media[END_REF] shows the similar trend for the case of coupled conductive-radiative heat transfer within an absorbing/emitting medium placed between two diffuse, non-black, infinite, isothermal parallel plates separated by a finite distance. However, for this test case, the temperature profile is slightly different than the one obtained by Viskanta et al. This could be explained by the fact that the geometries and temperature used in this test case is different than the one used by Viskanta et al. Nonetheless, on the basis of results obtained corresponding to different values of Stark number, it can be concluded that the numerical method developed in this work is able to capture the coupled physics accurately.

Test case 2: Kelvin-cell foam

The main objective of this test case is to present the ability of the solver to consider two different boundary conditions: mirror and vacuum, as well as to handle strong temperature gradient over the boundaries. Also, in this test case, we further increased the complexity both in terms of geometry and number of voxels (or unknowns). In this case, a network of Kelvin-cells was considered. A Kelvin-cell which is also known as tetrakaidecahedron is a polygon composed of fourteen faces: eight hexagonal and six quadrilateral, and twenty four edges. Due to structural simplicity and easy generation of samples with different numbers of cells and porosity, Kelvin-cell network has been the subjects of various studies including [START_REF] Badri | Conductive-radiative heat transfer within sic-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis[END_REF][START_REF] Sans-Laurent | Modeling the flash method by using a conductoradiative monte-carlo method: Application to porous media[END_REF][START_REF] Pabst | Young's modulus and thermal conductivity of closed-cell, open-cell and inverse ceramic foams -model-based predictions, cross-property predictions and numerical calculations[END_REF][START_REF] Caliot | Combined conductive-radiative heat transfer analysis in complex geometry using the Monte Carlo method[END_REF].

Geometry and input parameters for test case 2

The Kelvin-cell foam used for the computation is generated in "genMat" as in [START_REF] Badri | Discrete-scale numerical radiative transfer analysis of porous participating media[END_REF] and is shown in Figure 10 with applied boundary conditions to solve coupled heat transfer. The geometry measures 2.5×2.5×2.5 cm 3 and is composed of 250 voxels in each direction giving total number of voxels equal to ≈15.6 × 10 6 and number of triangles equal to 3 × 10 6 . The foam has a porosity of 90 %. Similar to previous test cases, the geometry was placed between two walls whose faces were maintained at fixed temperatures (T hot = 1800 K and T cold = 1200 K) to generate strong temperature gradient. We assumed that the solid phase is opaque and optically smooth with solid phase thermal conductivity (λ s ) equal to 0.3 W m -1 K -1 , void phase is considered to be non-conductive with thermal conductivity (λ v ) equal to 0 W m -1 K -1 , and a fixed total emissivity (ϵ) equal to 0.9. For this simulation, 96 (3.2 GHz) processors were used.

Results of test case 2

Effect of boundary condition: the plane-averaged temperature profile of a Kelvin-cell foam for both vacuum and mirror boundary conditions (on side walls) is shown in Figure 11. The number of points at which the temperature is computed is equal to the number of divisions in the direction of heat flow (plus two more points, one at each boundary).

For vacuum boundary condition on sides, the temperature in the middle goes lower than those fixed at hot side (T =1800 K) and cold side (T =1200 K). This occurs due to the loss of heat resulting from radiation or energy leak from the system. In case of vacuum boundary condition, the side walls are set free i.e., any radiation emitted from the geometry falling on any side walls will be directed to the space outside of the geometry and will be lost forever. This is similar to the case of black walls at absolute zero. In other words, it can be considered that the radiation from the material is being exchanged with black walls maintained at absolute zero temperature. It is evident that heat conduction through solid maintains the continuous flow of heat within the system otherwise, entire energy will be leaked out to the outer space through radiation. Also, due to the heat loss to the outer space, the deviation of T profile from the linear profile is significantly high with more significant loses close to the hot boundary.

For mirror boundary condition on sides, the temperature does not go beyond those fixed at hot and cold sides because the system remains closed i.e., no energy leaves the system. Further, we noticed asymmetric behavior of temperature profile as in previous case. The average temperature increases near the cold boundary as high intensity radiation coming from the hot region is absorbed here and vice-versa occurs near the hot boundary as explained by Viskanta et al. [START_REF] Viskanta | Boundary layer in thermal radiation absorbing and emitting media[END_REF].

The results obtained using our approach were compared with the results produced using a commercial software STAR-CCM+, also presented in Figure 11. STAR-CCM+ is also based on FVM and uses ray tracing to calculate the view factor instead of radiative exchange factor as computed in this work. The difference is that multiple reflections are not considered in case of view factor computation. For the computation in STAR-CCM+, the geometry was remeshed and the coupling was performed at local scale i.e., at the level of triangles. From Figure 11, it is evident that the results are in good agreement with those of STAR-CCM+, hence cross-validating our approach. Furthermore, in case of vacuum boundary conditions on side, an oscillatory behavior can be observed in Figure 11 (see location 0.1, 0.3, ..., 0.9). It signifies that the center parts of Kelvin-cells leak more radiation. Qualitatively, in case of radiative heat transfer within Kelvin-cell foam (similar to the geometry used in this test case), similar behavior was also observed in [START_REF] Badri | Discrete-scale numerical radiative transfer analysis of porous participating media[END_REF][START_REF] Badri | Efficient finite element strategies for solving the radiative transfer equation[END_REF]. Also, it is worth mentioning that the total computation time required to solve each case was less than 10 min [START_REF] Kumar | A numerical method based on domain decomposition to solve coupled conduction-radiation physics using parallel computing within large porous media[END_REF].

Effect of number of directions: the number of directions (or rays) used for exchange factor computation may have an effect on correct estimation of radiative heat transfer [START_REF] Badri | Discrete-scale numerical radiative transfer analysis of porous participating media[END_REF]. The results obtained using number of directions (N dir ) equal to 253 (angular sector = 10°) and 1169 (angular sector = 5°) are also shown in Figure 11. In case of Kelvin-cell foam, the effect was found to be negligible for both types of boundary conditions. However, the computation time increased more than double by increasing the number of rays from 253 (and total number of initial rays equal to 7.6 × 10 8 ) to 1169 (and total number of initial rays equal to 3.5 × 10 9 ). Hence, it can be concluded that 253 directions are enough for this test case.

Effect of number of divisions: similar to previous test case, the effect of number of divisions was also checked for Kelvin-cell foam. The effective solid thermal conductivity (λ eff cond ) values computed for different number of divisions is plotted in Figure 12. It can be seen that in case of porous media with one phase conduction, the number of divisions has significant effect. A maximum difference of up to 15 % between 1×1×1 and 7×7×7 subdivisions, as shown in bars can be seen in Figure 12.

Also, the total time taken for simulation for each case is also plotted by dots in Figure 12. The time first decreases with the increase in number of divisions. As mentioned in section 4, to solve effective solid thermal conductivity tensor of each subvolume, each subvolume is assigned to a single processor. With the increase in number of subdivision, size of each subvolume reduces which means less number of voxels per subvolume which leads to smaller matrix, thereby requiring less time to solve per subvolume. However, on further increasing the number of divisions (e.g., 6 × 6 × 6), total number of subvolumes increases, thereby increasing the time to read the geometry and building the matrix which increases the overall computation time. Furthermore, to study the effect of number of divisions (or size of the subvolume) on λ eff coup , two cases were selected: 5×5×5 and 7×7×7 corresponding to the maximum and minimum values of λ eff cond . T profile is plotted in Figure 13. Here, it is evident that number of divisions does not have any significant effect on coupled heat transfer, for both types of boundary conditions.

The mean heat flow rate throughout the thickness is calculated and tabulated in Table 1. The results obtained using STAR-CCM+ and our method show a very good agreement and confirm that the number of divisions do not affect the coupled heat transfer in this case. Finally, λ eff coup comes out to be approximately 0.31 W m -1 K -1 for all the cases. It should be noted that the λ eff coup /λ eff cond ≈6.2 at T ref = 1500 K and ∆T =600 K. This means that the contribution of radiation in coupled heat transfer is significantly high and a small variation in λ eff cond does not affect the coupled heat transfer. From the above discussion, it is clear that the effect of number of divisions (or size of the subvolume) on coupled heat transfer is negligible when the edge length of subvolume is comparable to the extinction mean free path of photons (Λ). The Λ for the sample was computed to be 334 µm using the Radiative Distribution Function Identification (RDFI) method [START_REF] Rousseau | Evolution of the homogenized volumetric radiative properties of a family of α-SiC foams with growing nominal pore diameter[END_REF] implemented in "genMat" software. For both the simulations performed using the number of divisions equal to 5×5×5 and 7×7×7, the side length of subvolumes was calculate to be 500 µm and 357 µm, respectively, and the ratio of side length of subvolume to the mean free path of photons (Λ) was calculated to be 1.5 and 1.0, respectively. Therefore, it can be concluded that a small variation in the ratio of side length of subvolume to the mean free path of photons (Λ) is acceptable as it does not make any effect on λ eff cond .

1 × 1 × 1 2 × 2 × 2 3 × 3 × 3 4 × 4 ×

Type

Test case 3: Virtual fibrous material

Finally, coupled conductive-radiative heat transfer is solved within highly porous virtual fibrous media. Unlike cellular porous media (e.g., Kelvin-cell foam), fibrous materials have fibers characterized by small diameter (varying around 10-50 µm which is approximately 10-100 times smaller than the ligament's size of foams) and uneven distribution. Also, the REV for such media is comparatively bigger than cellular materials which makes the estimation of coupled heat transfer more complex and challenging.

Geometry and input parameters for test case 3

In this case, two different geometries are considered. Geometries are generated using an in-house algorithm implemented in "genMat" software which allows to generate virtual fibrous materials samples with controlled parameters. The fibers are allowed to interpenetrate and their orientation is controlled. For random sample, the elevation/out-of-plane (yz-plane) angle was allowed to vary between 0-90°, however, in case of oriented media, the elevation angle is controlled using normal distribution with mean value fixed to 0°and standard deviation equal to 30°. Note that the azimuth/in-plane angle (ϕ) was uniformly varied from 0-360°. The geometry measures 2×2×2 mm 3 or 500×500×500 voxels giving total number of voxels equal to 125 × 10 6 and number of triangles ≈66 × 10 6 . Fibers are 200 mm long and 20 mm in diameter. Similar to previous test cases, the geometry was placed between two walls whose faces were maintained at fixed temperatures (T hot = 2000 K and T cold = 1500 K) to generate strong temperature gradient. We assumed that the solid phase is opaque and optically smooth with solid phase thermal conductivity (λ s ) equal to 2 W m -1 K -1 and void phase is considered to be non-conductive with thermal conductivity (λ v ) equal to 0 W m -1 K -1 . In this case, optical properties (n -ik) of the material is used to compute local emissivity (ϵ) and fixed to 3-i0.001. For this simulation, only 16 (3.5 GHz) processors were used.

Results of test case 3

Effect of orientation -Random vs. Oriented fibrous media: Plane-averaged temperature profile for both samples i.e., with randomly distributed and oriented fibers; with mirror boundary condition is shown in Figure 16. Like in previous case, the asymmetric behavior of temperature profile is observed. In addition to the plane-averaged temperature profile, the plane-averaged flux for conduction (c.f., equation ( 30)) and radiation (c.f., equation [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF]) is plotted in Figure 17. It can be seen that the conductive flux is high at the inlet (at solid plate) which first decreases suddenly, then it increases gradually within the medium until it reaches the maximum at the end (at solid plate). In contrast, the radiative flux behaves exactly the opposite. Radiative flux is small at the inlet and at the outlet, and decreases within the medium as it moves from hot boundary towards the cold boundary. Consequently, the total flux through the thickness remains almost constant and is also shown in Figure 17. It was also found that the total flux value is high in case of random media, which states that the heat transfer is high, and giving a high value of effective coupled thermal conductivity. The λ eff coup for random and oriented media comes out to be 0.306 W m -1 K -1 and 0.263 W m -1 K -1 , respectively. This states that the oriented fibrous sample is better insulator than the random fibrous sample. Effect on number of directions: Furthermore, the effect of number of directions was also checked. 841 (angular sector = 6°) and 361 number of directions (angular sector = 9°) per triangle were chosen and the results obtained are shown in Figure 16. It can be seen that the effect of number of directions is negligible in this case. However, the computation time increases by a factor of two while using more number of rays. Here, it can be concluded that 361 directions are enough for the computation.

Effect of number of divisions: Figure 18 shows the effect of N div on λ eff cond . It should be noted that in case of randomly distributed fibers, λ eff cond,x ≈ λ eff cond,y ≈ λ eff cond,z . However, in case of oriented media, λ eff cond,x < λ eff cond,y ≈ λ eff cond,z . This confirms that the fibers are uniformly distributed in all directions in case of random media whereas, majority of fibers are aligned in yz-plane in case of oriented media giving higher value of λ eff cond in y-and z-directions than in x-direction. Furthermore, as expected, N div has an effect on λ eff cond . Here, we will consider N div = 1, and resulting value of λ eff cond as reference. In case of random media, on varying N div to 5, 10, and 20, λ eff cond reduces by 7.0 %, 15.8 %, and 28.1 % of the reference value, respectively. Similarly, in case of oriented media, on varying N dir to 5, 10, and 20, λ eff cond reduces by 5.1 %, 10.1 %, and 19.5 % of the reference value, respectively. It is clear that the variation increases with the increase in N div .

Table 2 shows the effect of number of divisions on λ eff coup . As discussed in the last subsection, λ eff cond decreases with increase in number of divisions, due to the lack of connectivity. The results for λ eff coup also show the same trend. In case of random media, on increasing the number of divisions from 5 to 10, λ eff cond decreases by 5.5 % and λ eff coup decreases by 7.8 %. Similarly, in case of oriented media, on increasing the number of divisions from 5 to 10, λ eff cond decreases by 14.6 % and λ eff coup decreases by 17.5 %. From the above discussion, it is clear that number of subdivisions (or size of the subvolume) may have a significant effect on λ eff coup as shown in Table 2 but the effect is rather coming from λ eff cond due to loss of connectivity between solid subvolumes.

For both the computations performed using number of divisions equal to 5×5×5 and 10×10×10, the side length of the subvolumes was calculated to be 400 µm and 200 µm, respectively. Also, mean free path of photon (Λ) was calculated to be 62 µm for random sample and 68 µm for oriented sample. Then the ratio of side length of subvolume to the mean free path of photons (Λ) was calculated to be ≈6 and ≈3, respectively. Since the ratio is close to one for N div equal to 10×10×10, the N div equal to 10×10×10 (or size of subvolumes equal to 200 µm) is a good choice over 5×5×5 (or size of subvolumes equal to 400 µm). However, using N div equal to 20×20×20 would give ratio equal to 1.5 but it leads to low λ eff cond due to lack of connectivity between solid subvolumes as discussed in Figure 2, therefore, it is crucial to choose an appropriate number of divisions for the computation to obtain better results.

Type

Random

Conclusion

A numerical approach based on modified zonal method is proposed in this article to solve coupled conductive-radiative heat transfer at local scale within large scale digitized porous materials. To solve conduction-radiation coupling, the volume is divided into subvolumes and for each subvolume, effective thermal conductivity tensor (λ) and radiative exchange factor (F) between subvolumes are computed, and are used to perform coupling at the level of subvolumes. Note that these two quantities i.e., λ and F need to be computed only once (knowing that the optical properties are independent of temperature) and can be computed a priori, thereby, reducing the solving time at different temperatures and coupling efforts drastically. In addition to that, modified zonal method developed in this work leads to an embarrassingly parallel system as each subvolume is independent to each other i.e., no communication needed between subvolumes while computing λ which simplifies the code and further reduces the computation time.

The results of two test cases shown in this article confirm that the method is able to capture the physics and is able to produce accurate results. Finally, the method was applied to the two samples of virtual fibrous media namely random and oriented. We first confirmed that the method is not sensitive to the number of divisions considering that the size of subvolumes is equivalent to the mean free path length of extinction of thermal radiation. We also checked the effect of boundary conditions on temperature distribution within the media. It was found that boundary conditions significantly affect the temperature distribution within the media. However, the effective coupled thermal conductivity for the same geometry remains the same in case of vacuum as well as mirror boundary conditions. In case of virtual fibrous samples, results exhibit the strong impact of orientation and confirms that the sample with oriented fibers is better insulator of heat than that of random fibers.

Another advantage of the solver is that it is able to consider large temperature gradients over the boundary. The mathematical simplification used for the linearization allows to consider large temperature gradient between hot and cold boundaries, and provides accurate results. Furthermore, the method is pure deterministic in nature for solving both conduction and radiation physics, therefore, the results are free from any stochastic errors. Last but not the least, the excessive number of directions or rays used for the computation of radiative exchange factor does not affect the coupled heat transfer significantly as long as the total number of rays are very high (e.g., in the order of 1 × 10 8 to 1 × 10 9 ).

Overall, the method is able to solve coupled conductive-radiative heat transfer at local scale within large scale complex geometries having hundreds of millions of unknowns while dealing with different boundary conditions and large temperature gradients over the boundary. Knowing that the modified zonal method developed in this work to solve coupled conduction-radiation physics reduces the coupled system of equations drastically, we aim to further extend the capability of the solver to simulate heat transfer within "realistic" micrometric sized, heterogeneous, fibrous media having constituents sizing less than 10 µm which produces a 3D digitalized image of the size of REV consisting of approximately 1000×1000×1000 voxels.

Figure 1 :

 1 Figure1: A 2D representation of the subdivision of a volume into subvolumes. (i) and (j) denote the identity of subvolumes. Also, each subvolume consists of both solid and void voxels. For the sake of brevity, a few number of voxels are shown but in reality, each subvolume may contain millions of voxels.

Figure 2 :

 2 Figure 2: A control volume (i) surrounded by six neighbors to which the conductive heat energy can be exchanged. Here, the terms i denote the identity of the control volume at which the conductive heat balance is computed and E, W , N , S, F , and B denote the East, West, North, South, Front, and Back coordinates, respectively.

Figure 3 :

 3 Figure3: 2D Schematic of coupled system having 4×4 subvolumes only, assembled using the resistance network with possible boundary conditions required to solve coupled conduction-radiation physics. Boundary conditions written in black are those required to solve conductive heat transfer and in orange are those required to solve radiative heat transfer.

Algorithm 1 : 7 8 9 n = n + 1 10 end while 11 Calculate 12 m = m + 1 13

 1789111121 Two step iterative algorithm for coupled conduction-radiation modeling.Input: geometry, number of subdivisions (Nx, Ny, Nz), physical parameters: λ, F, ϵ, σB, . . . , numerical parameters: ε =1 × 10 -12 , δ1 = δ2 =1 × 10 12 , m = n = 0, . . . , and Boundary conditions 1 Initialize temperature T m (x ∈ Ω): Linearly distributed between T hot and T cold boundaries (using linear interpolation) 2 while δ1 > ε ▷ Outer loop 3 Solve coupled HCE (equation (19)) with given boundary conditions to determine → T n+1 (x ∈ Ω) Calculate residual δ2 = ∥T n+1 -T n ∥∞ residual δ1 = ∥T m -T n ∥∞

Figure 4 :

 4 Figure 4: Schematic of ray tracing procedure to compute radiative exchange factor. Emitted ray from triangle in subvolume (i) is shown in orange and the reflected ray from triangle in subvolume (j) is shown in blue. Note: For conciseness, only two subvolumes and only one surface per subvolume is shown.

Figure 5 :

 5 Figure 5: Geometry with added walls at hot and cold sides, and boundary conditions.

Figure 6 :

 6 Figure 6: Effect on number of divisions on solid conduction. Here, N div denotes the number of divisions in each directions.

Figure 7 :

 7 Figure 7: Effect of number of divisions on coupled heat transfer.

Figure 8 :

 8 Figure 8: ∆T profile for N div = 16.

Figure 9 :

 9 Figure 9: Effect of the Stark number on T in case of coupled heat transfer.

Figure 10 :

 10 Figure 10: Geometry of a Kelvin-cell foam and boundary conditions.

Figure 11 :

 11 Figure 11: Plane-averaged temperature (T ) profile for coupled conductive-radiative heat transfer computation for vacuum boundary condition (top) and mirror boundary condition (bottom) applied on sides. The images showing the temperature distribution within Kelvin-cell foam are generated in STAR-CCM+ .

Figure 12 :

 12 Figure 12: Effect of number of divisions on λ eff cond .

Figure 13 :

 13 Figure 13: Effect of number of subdivisions on coupled conductive-radiative heat transfer for vacuum boundary condition (top) and mirror boundary condition (bottom) applied on sides.

Figure 14 :

 14 Figure 14: Geometry of virtual fibrous sample and boundary conditions.

  (a) Randomly distributed fibers (b) Oriented fibers with θ ± σ = 0°±30°F igure 15: Cross-sectional (yz-plane) view of two samples used for the cross-validation showing the angular distribution of fibers. Longer cross-section on right side image states that the fibers are more aligned along the yz-plane. The cross-section of both the images is 2×2 mm 2 .

Figure 16 :

 16 Figure 16: Plane-averaged temperature profile for random fibrous sample (top) and oriented fibrous sample (bottom).

4 X

 4 * [-]

Figure 17 :

 17 Figure 17: Plane-averaged flux profile for random fibrous sample (top) and oriented fibrous sample (bottom).

Figure 18 :

 18 Figure 18: Effect of number of subdivisions on λ eff cond for random fibrous sample (left) and for oriented fibrous sample (right).

Table 1 :

 1 Effect of number of divisions. N dir was fixed to 253.

		Vacuum BC	Mirror BC
	N div	5×5×5	7×7×7	5×5×5	7×7×7
	Mean heat flow rate [W]	0.445	0.446	0.445	0.447

Table 2 :

 2 Effect of number of divisions on λ eff coup . N dir was fixed to 361.

	media	Oriented media

Acknowledgement

The author Atin Kumar would like to thank Mersen, Gennevilliers, France and National Association for Research and Technology (ANRT) for funding this research during his Ph.D. under the Industrial Agreement of Training through Research (CIFRE). The author also thanks Dr. Sophia Haussener, EPFL, Switzerland for fruitful and valuable discussions while developing the approach.

Competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Gérer et implanter des formations ligneuses sur une ferme : Que voyez-vous sur la ferme expérimentale INRAE ? De la gestion à l'implantation d'arbres sur la ferme expérimentale 1

• • •

Angles de vue (voir photos ci-dessous) 1