Morgan Sans

Olivier Farges

Vincent Schick

Gilles Parent
email: gilles.parent@univ-lorraine.fr

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

Keywords:

This paper presents a single Monte-Carlo algorithm used to solve transient conductive and radiative heat transfers in three-dimensional porous media. The complete methodology presented step by step herein enables practical and efficient study of geometrical and multiphysical complexities. The code was validated against results obtained by commercial software, analytical and semi-analytical solutions. Computation times were found to be greatly reduced when radiative transfer is predominant compared with those obtained using a deterministic solver. This kind of approach allows a probe calculation in the frame of linear thermal transfers and is well suited for the numerical characterization of heterogeneous media. In this work a numerical flash method was reproduced and enabled us to evaluate the effective total conductivity of the equivalent homogenized medium. The influence of various parameters such as porosity, size of the unit cell, bulk conductivity of the solid phase, reference temperature and emissivity was studied for a stack of Kelvin cells. This tool enables the parametric investigation of geometric and thermal properties. The results are in good agreement with those of the literature.

Nomenclature

𝐧

Unitary outward normal vector (-)

Introduction

Ceramic and metallic foams are commonly used in many engineering applications because of their interesting thermal, mechanical and optical properties [START_REF] Sommers | Ceramics and ceramic matrix composites for heat exchangers in advanced thermal systems-A review[END_REF]. Improved manufacturing processes enable the production of porous media with a very wide range of structural parameters such as porosity, cell size distribution, geometries or materials used. These can thus be designed and used for several applications such as catalysis, volumetric solar absorbers, insulators, heat exchangers or gas-phase heat recovery. The foam structural and thermal parameters play a key role in the efficiency of the studied systems because of the coupling between the different physical phenomena involved such as conduction in the fluid/solid phase, convection, transport, radiation or chemical reaction. Their optimization necessarily requires better knowledge of the link between the porous heterogeneous medium and the coupled heat transfer mechanisms and this can be achieved by using numerical and experimental tools.

Nowadays, combined heat transfers can be solved directly on complex geometry at the sample scale, i.e. at a scale equal to a few (dozen of) size of pores, but this becomes difficult at the scale of the system itself. Therefore, it is very common for the thermal problem to be simplified by considering foam as a homogeneous medium with equivalent thermal properties, which represents the overall thermal transport through the medium and has to be characterized. In this framework, a unique Equivalent Thermal Conductivity (ETC), which can include conduction, convection and radiation has been widely considered. Interested readers can refer to [START_REF] Ranut | On the effective thermal conductivity of aluminum metal foams: Review and improvement of the available empirical and analytical models[END_REF], which gives a broad overview of the studied geometries in different scopes of validity. Analytical solutions or numerical tools have been provided for the simplest geometries [START_REF] Landauer | The electrical resistance of binary metallic mixtures[END_REF] while numerical or experimental characterization procedures have been developed to study 3D complex porous geometries. The thermal behavior of metallic foams at ambient temperature were studied by Coquard and Baillis [START_REF] Coquard | Numerical investigation of the radiative properties of polymeric foams from tomographic images[END_REF] or by R. Wulf et al. [START_REF] Wulf | Experimental and numerical determination of effective thermal conductivity of open cell FeCrAl-alloy metal foams[END_REF] who used a numerical hot guarded plate method from tomographic images. Zhao et al. [START_REF] Zhao | The temperature dependence of effective thermal conductivity of open-celled steel alloy foams[END_REF], Fend et al. [START_REF] Fend | Porous materials as open volumetric solar receivers: Experimental determination of thermophysical and heat transfer properties[END_REF] and Coquard et al. [START_REF] Coquard | Experimental investigations of the coupled conductive and radiative heat transfer in metallic/ceramic foams[END_REF] used respectively an experimental guarded hot plane, a transient plan source and a flash method. Such studies can be used to explore the validity of the theoretical models described above, validate employed characterization procedures and highlight the main thermal or geometrical properties in heat transport.

Ceramic and metallic foams are particularly interesting for hightemperature applications in which radiation transfer plays a major role. The description of the radiative behavior in such media through the identification of equivalent optical properties has mostly been studied independently of the other modes of heat transport. Authors like Tancrez and Taine [START_REF] Tancrez | Direct identification of absorption and scattering coefficients and phase function of a porous medium by a Monte Carlo technique[END_REF], Guévelou et al. [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] or Loretz et al. [START_REF] Loretz | Metallic foams: Radiative properties/comparison between different models[END_REF] proposed numerical methods based on ray-tracing approaches and provided correlations from tomographic images while authors like Hendricks and Howell [START_REF] Hendricks | Absorption/scattering coefficients and scattering phase functions in reticulated porous ceramics[END_REF] or Baillis et al. [START_REF] Baillis | Identification of spectral radiative properties of polyurethane foam from hemispherical and bi-directional transmittance and reflectance measurements[END_REF] performed a spectral estimation of equivalent radiative properties from spectrometric measurements. Characterization following an entirely numerical procedure requires efficient modeling of the coupled thermal transfers directly on the complex and 3D heterogeneous medium. Solving this problem can remain quite tricky and time-consuming because of the different nature of the two modes of heat transfer. Indeed, sufficiently finely resolved and realistic geometry and optimized solvers for each of the transport modes are required to achieve this. To our knowledge, very few studies have focused on the numerical identification of equivalent properties of foams with coupled conductive and radiative transfer. Some authors like Mendes et al. [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] have proposed using a numerical hot guarded plate applied to metallic foams at high temperature and demonstrated the limits of using the Rosseland approximation to describe the equivalent radiative conductivity. Subsequently, these same authors developed a complete numerical study [START_REF] Mendes | Detailed and simplified models for evaluation of effective thermal conductivity of open-cell porous foams at high temperatures in presence of thermal radiation[END_REF] to determine the possible bias generated by simplified radiative modeling of the medium (heterogeneous, homogeneous, Rosseland...). Empirical laws were then provided to make modeling these heterogeneous environments easier. Very recently, Vignoles and Ortona [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] proposed an original stochastic method to model both conduction and radiation in ceramic/metallic foams and fibers and identified total equivalent properties.

Experimental measurements of the conductivity of foams at high temperature are not common at all. The characterization procedure requires the use of a fast and accurate coupled model. Zhao et al. [START_REF] Zhao | The temperature dependence of effective thermal conductivity of open-celled steel alloy foams[END_REF], Coquard and Baillis [START_REF] Coquard | Modeling of heat transfer in low-density EPS foams[END_REF] and Mendes et al. [START_REF] Mendes | Experimental validation of simplified conduction-radiation models for evaluation of effective thermal conductivity of open-cell metal foams at high temperatures[END_REF] used the method involving a hot guarded plate while Coquard et al. [START_REF] Coquard | Experimental and theoretical study of the hot-wire method applied to low-density thermal insulators[END_REF] used the hotwire to study coupled heat transfers in expanded polystyrene foams. Niezgoda et al. [START_REF] Niezgoda | Modeling of time-resolved coupled radiative and conductive heat transfer in multilayer semitransparent materials up to very high temperatures[END_REF] applied the flash method to the study of a silica aerogel and other thermal insulators such as low density foams. However, due to a lack of sensitivity or to correlations between parameters, the identification of equivalent radiative properties was difficult and was performed separately. For the very first time, Coquard et al. [START_REF] Coquard | Experimental investigations of the coupled conductive and radiative heat transfer in metallic/ceramic foams[END_REF] simultaneously characterized an equivalent conductivity and equivalent absorbing and scattering coefficients for metallic and ceramic foams (Zirconia, FeCr alloy, Mullite and NiCrAl) over a temperature range between 296 K and 673 K. The authors of the present paper recently reported this kind of characterization at temperatures up to 1000 K applied to SiC and SiSiC foams taking into account a wide variety of influential parameters such as the geometry of cells, pore diameter, bulk conductivity, the reference temperature and so on [START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams at high temperature[END_REF]. Nevertheless, this kind of experimental characterization is always time-consuming and costly, especially when aiming high temperatures. Moreover, the large number of possible combinations of intrinsic parameters makes a systemic parametric exploration difficult as this requires the right samples to be available. A numerical approach is thus necessary to assist with experimental characterization and anticipate needs.

As explained above, the numerical resolution of coupled heat transfers directly on complex geometry remains tricky and requires efficient numerical tools. To our knowledge, deterministic methods are commonly used such as the Finite Element Method (FEM) for conduction and the Discrete Ordinate Method for radiation [START_REF] Wulf | Experimental and numerical determination of effective thermal conductivity of open cell FeCrAl-alloy metal foams[END_REF][START_REF] Mendes | Detailed and simplified models for evaluation of effective thermal conductivity of open-cell porous foams at high temperatures in presence of thermal radiation[END_REF][START_REF] Patel | Development of correlations for effective thermal conductivity of a tetrakaidecahedra structure in presence of combined conduction and radiation heat transfer[END_REF][START_REF] Badri | Vectorial finite elements for solving the radiative transfer equation[END_REF]. Recently, Fournier et al. [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm[END_REF] suggested the idea of solving multi-physical phenomena with a single Monte-Carlo algorithm. They also showed that this kind of algorithm can be constructed if the set of considered equations can be written as a Fredholm equation of the second kind corresponding to the linear heat transfer framework. Based on an integral formulation, a recursive (backward in time) algorithm computes the temperature at a given time and at a given location (probe calculation) of a complex 3D geometry. The different heat transfer modes are solved at once so there is no need to couple two independent solvers. As this method requires only a surface mesh, it is very practical to use and convenient for dealing with geometric and multiphysical complexities. Moreover, coupling this with tools developed by the image synthesis community (such as hierarchical grids) means performance levels can be increased by reducing computation times and through a demonstrated insensitivity to the mesh refinement [START_REF] Villefranque | A path-tracing Monte Carlo library for 3-D radiative transfer in highly resolved cloudy atmospheres[END_REF]. Caliot et al. [START_REF] Caliot | Combined conductive-radiative heat transfer analysis in complex geometry using the Monte Carlo method[END_REF] applied this algorithm in a stationary case involving a structured foam made of Kelvin cells. This approach has been validated for conductive and radiative transfers by comparison with a deterministic method. Ibarrart et al. [START_REF] Ibarrart | Combined conductive-convective-radiative heat transfer in complex geometry using the Monte Carlo method : Application to solar receivers[END_REF] extended this to consider coupled conductive, advective and radiative heat transfers while Penazzi et al. [START_REF] Penazzi | Toward the use of symbolic Monte Carlo for conduction-radiation coupling in complex geometries[END_REF] studied the case of a semi-transparent medium. Assuming the linearity of heat transfers is a not limiting factor because only a small increase of the equilibrium temperature of the medium is required in most of the characterization methods, including the flash method. Moreover, recent studies have shown that solutions can be found to manage nonlinear problems [START_REF] Tregan | Transient conducto-radiative heat transfer in a single Monte-Carlo algorithm: Handling the nonlinearity[END_REF]. Finally, the Monte-Carlo algorithm and in particular its probe calculation aspect as proposed by Fournier et al. [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm[END_REF] appears highly suitable to use in the field of thermal characterization as it provides the solution to the direct model on the complex geometry. In this paper, we propose to apply the algorithm to the particular case of the flash method and to demonstrate the high potential of this approach for the study of coupled transfers in a heterogeneous medium.

In the following, we shall first describe the thermal problem and governing equations concerned by our work. Secondly, we shall go on to describe the generic principle of the Monte-Carlo algorithm, which allows geometrical and multiphysical complexities to be dealt with in a very practical and efficient way. The complete methodology and the particular configuration of the flash method are detailed. Compared to previous studies [START_REF] Caliot | Combined conductive-radiative heat transfer analysis in complex geometry using the Monte Carlo method[END_REF][START_REF] Ibarrart | Combined conductive-convective-radiative heat transfer in complex geometry using the Monte Carlo method : Application to solar receivers[END_REF][START_REF] Penazzi | Toward the use of symbolic Monte Carlo for conduction-radiation coupling in complex geometries[END_REF], the present problem requires a transient solution. Thus, a validation case is provided for a 3D geometry which could be useful as a benchmark case for further algorithmic improvements. A comparison with a calculation performed by a commercial software based on deterministic methods also helps highlight the interests of our approach. Next the construction of the geometry of foams is presented and the identification procedure based on the numerically obtained thermograms is briefly described. Finally, a parametric study is presented. The influence of main parameters is studied including porosity, size of the cells, bulk conductivity of the foam struts, the bulk emissivity and the reference temperature on the total heat transfer.

Methods

As explained previously, the flash method is an efficient method for measuring the total equivalent conductivity of a semi-transparent medium. In addition, the transient nature of this characterization method means it is possible to separate different contributions and the authors successfully achieved this to characterize experimentally silicon carbide foams at high temperature [START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams at high temperature[END_REF]. The idea of the present paper is to replace an experimental flash method by an entirely numerical one. Although the guarded hot plate or hot wire methods appear as references and are also well suited for such an identification of thermal properties, only the flash method is discussed here to benefit from the previous experience developed and allow comparison. However, there is no limit to the extension of the methodology described in this section to other characterization techniques.

In the experimental flash method (cf. Fig. 1(a)), a sample at thermal equilibrium 𝑇 𝐼 is excited with a short duration heat flux 𝛷(𝑡) on its front face. The transient evolution of the mean surface temperature of the rear face 𝑇 (𝐱 obs , 𝑡) is measured with an infrared detector or an infrared camera. It is common to place the sample between two soleplates or to cover each face with a black coating [START_REF] Coquard | Experimental investigations of the coupled conductive and radiative heat transfer in metallic/ceramic foams[END_REF][START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams at high temperature[END_REF] Here this is necessary because of the semi-transparent nature of the sample and because of the experimental configuration used. The first soleplate absorbs the incident radiative heat flux and avoids the direct transmission of the laser excitation to the detector. The second is useful to avoid measuring volumetric emission. This means that a proper spatially averaged temperature can be measured. Finally, the porous medium is considered under vacuum, such that no natural convection or conduction through the fluid phase can occur. This choice was only made to isolate and identify conduction through the solid phase but no experimental or numerical limitations prevent such study of the influence of the fluid phase. After a presentation of the context of the study, the following section is dedicated to the description of the methodology used to provide the coupled direct model describing the heat transfer through the heterogeneous medium.

Monte-Carlo Method: General case

Fundamentally, the Monte Carlo (MC) method is a method that estimates integrals (and discrete sums) by interpreting them as expectations. Any quantity that can be expressed in an integral form can be estimated using a MC algorithm, including differential problems that are initially formulated outside the framework of stochastic processes [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF][START_REF] Kac | On distributions of certain Wiener functionals[END_REF]. In this paper, our work aimed to compute the temperature at location 𝐱 obs and at time 𝑡, by following a large number of paths crossing the complex geometry. The temperature 𝑇 (𝐱 obs , 𝑡) is then the expectation of a random variable 𝛩(𝐱 obs , 𝑡):

𝑇 (𝐱 obs , 𝑡) = [𝛩(𝐱 obs , 𝑡)] ≈ 1 𝑁 mc 𝑁 mc ∑ 𝑖=1 𝜃 𝑖 (1)
with 𝑇 (𝐱 obs , 𝑡) the Monte-Carlo estimator of 𝑇 (𝐱 obs , 𝑡), 𝛩 a random variable and 𝜃 𝑖 the ith realization of this variable, which is determined from the ith path inside the calculation domain. A path starts at the location 𝐱 obs and at the time 𝑡 where and when the temperature needs to be computed. A path stops when a known temperature such as the initial condition or a Dirichlet boundary condition is reached. The achieved temperature is retained as the weight of the ith realization of the algorithm.

Monte-Carlo Method: Conductive path

Concerning the conduction, the transient Energy Balance Equation (EBE) is:

𝜕𝑇 𝜕𝑡 = 𝛼∇ 2 𝑇 ∀𝐱 ∈  𝑠 (2)
One way to solve the given thermal problem is to use the finite difference method. In the case of a unidirectional problem, using an explicit Euler scheme (order 1) to approximate the time derivative, and a centered scheme (order 2) for the spatial derivative, the EBE (cf. Eq. ( 1)) becomes :

𝑇 (𝑥, 𝑡) -𝑇 (𝑥, 𝑡 -𝜏) 𝜏 = 𝛼 𝛿 2 [𝑇 (𝑥 + 𝛿, 𝑡 -𝜏) -2𝑇 (𝑥, 𝑡 -𝜏) + 𝑇 (𝑥 -𝛿, 𝑡 -𝜏)] (3)
with 𝛿 the spatial discretization step and 𝜏 the temporal discretization step. The temperature at position 𝑥 and time 𝑡 is thus expressed as:

𝑇 (𝑥, 𝑡) = 𝜇𝑇 (𝑥 + 𝛿, 𝑡 -𝜏) + (1 -2𝜇)𝑇 (𝑥, 𝑡 -𝜏) + 𝜇𝑇 (𝑥 -𝛿, 𝑡 -𝜏) (4)
with 𝜇 = 𝛼𝜏∕𝛿 2 . For this example, the Courant-Friedrichs-Lewis (CFL) condition requires 𝜇 ≤ 1∕2 for the explicit scheme to be stable. Following this condition, Eq. ( 4) allows a probabilistic interpretation and is used to construct propagation paths:

𝑇 (𝑥, 𝑡) = 𝑝 𝑟 𝑇 (𝑥 + 𝛿, 𝑡 -𝜏) + 𝑝 𝑐 𝑇 (𝑥, 𝑡 -𝜏) + 𝑝 𝑙 𝑇 (𝑥 -𝛿, 𝑡 -𝜏) (5)
It is then possible to construct a walk based on the corresponding 1D regular mesh of step 𝛿, with 𝑝 𝑟 , 𝑝 𝑐 and 𝑝 𝑙 the probabilities associated with a computation of the temperatures 𝑇 (𝑥 + 𝛿, 𝑡 -𝜏) (on the right), 𝑇 (𝑥, 𝑡 -𝜏) (in the center) and 𝑇 (𝑥 -𝛿, 𝑡 -𝜏) (on the left), respectively. Each move 𝛿 requires a time decrement 𝜏. If the temperature at a new location at time 𝑡 -𝜏 is unknown, the process is thus continued recursively until a known temperature such as a boundary (Dirichlet) temperature or the temperature at the initial time is reached.

The transition from this kind of walking technique to a threedimensional grid is straightforward. It was initially proposed by Haji-Sheikh and Sparrow [START_REF] Haji-Sheikh | The floating random walk and its application to Monte Carlo solutions of heat equations[END_REF] and allow to simulate a Brownian motion process. However, the previous example is described so that the purpose of the walk building process can be better understood. This type of walk is actually not practical when dealing with boundary conditions which are different from Dirichlet in complex geometry because of the derivative terms with respect to the normal. Therefore the same authors proposed a walk based on a sphere. In a spherical coordinate system, the formal solution of Eq. (1) at a given position 𝐱 (here 𝐱 is the position in 3D and has thus 3 coordinates) and at time 𝑡 is:

𝑇 (𝐱, 𝑡) = ∫ 1 𝐹 =0 ∫ 1 𝐺=0 ∫ 𝑡 𝜏=0 𝑇 (𝐱 + 𝛿𝐮, 𝑡 -𝜏)d𝐹 (𝜑)d𝐺(𝜃)d𝐻(𝜏) (6)
{ 𝐹 (𝜑) = 𝜑∕2𝜋, 𝐺(𝜃) = 1 2 (1 -cos(𝜃)) 𝐻(𝜏) = 1 + 2 ∑ ∞ 𝑝=1 (-1) 𝑝 exp(-𝑎𝑝 2 𝜋 2 𝜏∕𝛿 2 ) (7)
with 𝛿 the radius of a sphere centered in 𝐱, 𝜃 the polar angle and 𝜑 the azimuthal angle. As in the previous example, Eqs. ( 6) and ( 7) also admit a probabilistic interpretation. The quantities 𝐹 , 𝐺 and 𝐻 can be regarded as cumulative density functions of the azimuthal angle 𝜑, the polar angle 𝜃 and the time decrement 𝜏, respectively. The latter expression can be thus rewritten:

𝑇 (𝐱, 𝑡) = ∫ 4𝜋 𝑝 𝐔 (𝐮)d𝐮 ∫ ∞ 0 𝑝 𝑇 (𝜏)d𝜏 [(𝑡 -𝜏 ≤ 0) 𝑇 𝐼 + (𝑡 -𝜏 > 0) 𝑇 (𝐱 + 𝛿𝐮, 𝑡 -𝜏)] ( 8
)
with 𝑝 𝐔 (𝐮) and 𝑝 𝑇 (𝜏) the probability densities associated with the sampling of a direction 𝐮 and a time 𝜏, respectively.  is the Heaviside function. Three random numbers uniformly distributed between 0 and 1 allow us to sample a direction 𝐮 and a time 𝜏. In the case of an infinite medium, the corresponding algorithm to Eq. ( 8) is described in algorithm 1: Hence, if the initial condition is not reached, the path is placed in 𝐱 + 𝛿𝐮 and the time 𝑡 associated with the path is decremented by the quantity 𝜏. A random walk based on the construction of a sphere is generated through iteration as shown in Fig. 1(b). In the case of a finite medium, the conductive path stops if the initial condition or if a boundary where temperature is known is reached.

Algorithm

Monte-Carlo Method: radiative path

This paper only deals with the case of an opaque medium and a transparent void phase to simulate experimental conditions of our previous study [START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams at high temperature[END_REF]. However, the present method is not only limited to this type of interface condition and can be extended, for example, to the case of a semi-transparent medium or to a convective exchange.

In our case, the radiative transfer is thus limited to a surface-tosurface exchange. The coupling between conduction and radiation is thus carried out at the solid/void interface 𝜕 sf . The balance of the conductive and the radiative fluxes is written as follows: with 𝜖 the emissivity of the surface (uniform and hemispherical), 𝐧 the unitary outward vector normal to the surface, 𝜆 𝑠 the thermal conductivity of the solid, 𝐱 𝑏 the position at the solid/void boundary, 𝐯 the direction, 𝐼 0 the blackbody intensity and 𝐼 the intensity. For the sake of simplicity, radiative quantities averaged on the whole spectrum are considered. As described in [START_REF] Caliot | Combined conductive-radiative heat transfer analysis in complex geometry using the Monte Carlo method[END_REF], the wall temperature is non-linearly dependent on the blackbody intensity which means the latter must be linearized around a reference temperature 𝑇 ref :

𝐧 ⋅ 𝜆 𝑠 𝛁𝑇 = ∫ 2𝜋 |𝐯 ⋅ 𝐧| 𝜖 ( 𝐼 0 (𝐱 𝑏 ) -𝐼(𝐱 𝑏 , -𝐯) ) d𝐯 (9)
𝐼 0 (𝐱 𝑏 ) = 𝜎 𝜋 𝑇 4 (𝐱 𝑏 ) ≈ 𝜎 𝜋 𝑇 4 ref + 4 𝜎 𝜋 𝑇 3 ref (𝐱 𝑏 ) ( 𝑇 (𝐱 𝑏 ) -𝑇 ref ) (10)
with 𝜎 the Stefan-Boltzmann constant. As previously discussed, this assumption does not pose particular problems in the case of the thermal characterization methods used because of the low temperature elevation of the sample (just a few Kelvin) relative to the equilibrium temperature. As the void phase is transparent, the incident intensity at the wall 𝐼(𝐱 𝑏 , -𝐯) is equal to the intensity leaving from 𝐱 i in the direction -𝐯:

𝐼(𝐱 𝑏 , -𝐯) = 𝐼(𝐱 i , -𝐯) = 𝜖𝐼 0 (𝐱 i ) + ∫ 2𝜋 𝑝 R (𝐱 i )𝜌 ′ |𝐯 r ⋅ 𝐧 i | 𝐼(𝐱 i , -𝐯 r )d𝐯 r ( 11
)
with 𝑝 R (𝐱 i )𝜌 ′ the Bidirectional Reflectance Distribution Function, which is the product of a probability density over the reflected directions 𝐯 r and the directional-hemispherical reflectivity. 𝐧 i is the outward normal at the location 𝐱 i . According to Eq. ( 11), the incident intensity is therefore decomposed into an emission term and a reflection term integrated over the directions 𝐯 r . In this expression, the incident intensity 𝐼(𝐱 i , -𝐯 r ) can also be decomposed according to Eq. ( 11), which gives rise to nested integrals and leads to an iterative algorithm. Finally, we can express the incoming intensity as a function of an intensity emitted from an unknown position 𝐱 𝛾 . The latter intensity can reach the location 𝐱 𝑏 in the direction -𝐯 after zero, one or multiple reflections. For the sake of clarity, the corresponding nested multiple integrals are simply represented by one integral on the space of optical paths 𝛤 :

𝐼(𝐱 𝑏 , -𝐯) = ∫ 𝛤 𝑝 𝛤 (𝛾)𝐼 0 (𝐱 𝛾 )d𝛾 (12)
with 𝑝 𝛤 (𝛾) the probability density associated with an optical path 𝛾. By injecting Eqs. ( 10) and ( 12) (linearized) in Eq. ( 9), we can write:

-𝐧 ⋅ 𝜆 𝑠 𝛁𝑇 = 4𝜖𝜎𝑇 3 ref ( 𝑇 (𝐱 𝑏 ) -∫ 2𝜋 |𝐯 ⋅ 𝐧| 𝜋 d𝐯 ∫ 𝛤 𝑝 𝛤 d𝛾𝑇 (𝐱 𝛾 ) ) = ℎ rad ( 𝑇 (𝐱 𝑏 ) -∫ 2𝜋 𝑝 𝐕 (𝐯)d𝐯 ∫ 𝛤 𝑝 𝛤 d𝛾𝑇 (𝐱 𝛾 ) ) = ℎ rad ( 𝑇 (𝐱 𝑏 ) -𝑇 rad (𝐱 𝑏 ) ) (13)
with ℎ rad a radiative exchange coefficient. 𝑇 rad represents an average radiative temperature seen at the interface and due to the radiative exchanges through the void phase. The latter quantity can be estimated with the following algorithm (cf. algorithm 2), which is the corresponding algorithmic interpretation of Eqs. ( 12) and ( 13):

Algorithm 2 Sampling of a radiative path: estimate 𝑇 rad (𝐱 𝑏 )

Sample a direction of travel 𝐯 according to 𝑝 𝐕 Evaluate a first intersection 𝐱 i Sample canonically a random number

𝑟 in [0, 1] if 𝑟 ≤ 𝜖 then 𝑇 rad (𝐱 𝑏 ) = 𝑇 (𝐱 i ) (absorbed) else Sample a new direction 𝐯 r according to 𝑝 R 𝑇 rad (𝐱 𝑏 ) = 𝑇 rad (𝐱 i ) (recursivity) end if
Hence, an emission direction obeying a lambertian law is sampled when evaluating the radiative temperature. Then, until an absorption event is achieved, the path jumps through the void phase as represented in Fig. 2. At each new reflection, a new direction according to the BRDF is sampled. Each jump does not imply time decrementation because of the very high velocity of radiative transfers compared to other heat transfer modes such as conduction in the solid phase.

Monte-Carlo Method: Computation of the wall temperature

As was the case for the EBE, a finite difference scheme (order 1) requiring the introduction of an infinitesimal length 𝛿 𝑏 is used. The conductive heat flux in Eq. ( 13) is therefore approximated as:

-𝐧 ⋅ 𝜆 𝑠 𝛁𝑇 ≈ - 𝜆 𝑠 𝛿 𝑏 ( 𝑇 (𝐱 𝑏 ) -𝑇 (𝐱 𝑏 -𝛿 𝑏 𝐧) ) (14)
By injecting Eq. ( 14) in Eq. ( 13), the wall temperature is expressed:

𝑇 (𝐱 𝑏 ) = ℎ rad 𝜆 𝑠 ∕𝛿 𝑏 + ℎ rad 𝑇 rad (𝐱 𝑏 ) + 𝜆 𝑠 ∕𝛿 𝑏 𝜆 𝑠 ∕𝛿 𝑏 + ℎ rad 𝑇 (𝐱 𝑏 -𝛿 𝑏 𝐧) (15)
Hence, the temperature at the interface in Eq. ( 15) can be interpreted as the average of two temperatures weighed by two associated probabilities 𝑝 rad and 𝑝 cond :

𝑝 rad = ℎ rad 𝜆 𝑠 ∕𝛿 𝑏 + ℎ rad (16a) 𝑝 cond = 𝜆 𝑠 ∕𝛿 𝑏 𝜆 𝑠 ∕𝛿 𝑏 + ℎ rad (16b)
The corresponding Monte-Carlo algorithm is straightforward as described in algorithm 3:

Algorithm 3 Computation of the wall temperature 𝑇 (𝐱 𝑏 ) Sample canonically a random number 𝑟 in [0, 1] if 𝑟 ≤ 𝑃 𝑐𝑜𝑛𝑑 then 𝑇 (𝐱 𝑏 ) = 𝑇 (𝐱 𝑏 -𝛿 𝑏 𝐧) = 𝑇 (𝐱, 𝑡) (conductive path cf. algorithm 1) else 𝑇 (𝐱 𝑏 ) = 𝑇 rad (𝐱 𝑏 ) (radiative path cf. algorithm 2) end if

Monte-Carlo Method: Computation of the wall temperature (non diffusive wall)

Considering a given non-emitting (𝜖 = 0) wall 𝜕 𝛷 , excited by a transient heat flux 𝛷(𝑡), the balance at the interface is: As with the finite difference scheme described earlier, the wall temperature can be expressed as:

𝐧 ⋅ 𝜆 𝑠 𝛁𝑇 = 𝛷(𝑡) (17)
𝑇 (𝐱 𝑏 , 𝑡) = 𝑇 (𝐱 𝑏 -𝛿 𝑏 𝐧, 𝑡) + 𝛷(𝑡)𝛿 𝑏 ∕𝜆 s ( 18
)
This equation can be interpreted with a Monte-Carlo algorithm. If the path reaches an interface with 𝐱 𝑏 ∈ 𝜕 𝛷 , it is put back at the location 𝐱 𝑏 -𝛿 𝑏 𝐧 and the Monte-Carlo weight is increased by a source term equal to 𝛷(𝑡)𝛿 𝑏 ∕𝜆 s . This kind of procedure is repeated at each contact with the excited surface. In the same manner, if the impacted wall is nondiffusive and adiabatic (𝛷(𝑡) = 0), the path is simply put back at the location 𝐱 𝑏 -𝛿 𝑏 𝐧 and the Monte-Carlo weight is not increased.

Monte-Carlo Method: Random walk techniques

As discussed previously, the generation of the random walk requires the introduction of two arbitrary lengths 𝛿 and 𝛿 𝑏 which respectively correspond to a displacement within the solid volume and a replacement at the interface. As for the deterministic methods, the computed solution with the MC algorithm converges to the exact solution for 𝛿 ⟶ 0 and 𝛿 𝑏 ⟶ 0. However, as the calculation time increases greatly when the displacement steps are decreased, a compromise is therefore required between the computing cost and precision of the solution.

Moreover, several walking techniques can be employed to mimic the brownian motion [START_REF] Haji-Sheikh | The floating random walk and its application to Monte Carlo solutions of heat equations[END_REF]. The values of displacement steps 𝛿 and 𝛿 𝑏 can be alternately fixed or floating. These two different techniques are represented in Fig. 3. The first approach consists of setting a fixed value for displacement step 𝛿 as shown in Fig. 3 (a). At each iteration, a direction of travel 𝐮 is sampled according to 𝑝 𝐔 and the distance 𝛿 𝑤 between the considered position 𝐱 and the boundary 𝜕 s in the direction 𝐮 is computed. If 𝛿 < 𝛿 𝑤 , the associated boundary condition is processed and the wall temperature 𝑇 (𝐱 𝑏 ) is computed. If the path has to be put back in the solid, a fixed value of 𝛿 𝑏 is used. If 𝛿 > 𝛿 𝑤 , the path is simply placed in 𝐱 + 𝛿𝐮 with the associated time decrement. Hence, the chosen displacement is computed in the form of 𝛿 = min(𝛿, 𝛿 𝑤 ).

As Fig. 3 shows, the second technique is based on a floating value of the displacement step 𝛿. The idea is to compute the radius 𝛿 𝑠 corresponding to the inner circle. The displacement step is chosen in the form of 𝛿 = min(𝛿 𝑎 , 𝛿 𝑠 ) with 𝛿 𝑎 being an arbitrary value chosen to prevent the displacement being too large. Then, a direction of travel 𝐮 is sampled and the jump is performed. If the considered position 𝐱 is located between the closest position 𝐱 𝑏 to the boundary 𝜕 s and 𝐱 𝑏 -𝛿 𝑏 𝐧, the associated boundary condition is treated. If the path has to be put back in the solid, the same value of 𝛿 𝑏 is used and the path is placed in 𝐱 -𝛿 𝑏 𝐧. Finally, each floating displacement 𝛿 in the volume also implies a time decrement according to the chosen value of 𝛿.

Between these two described walking techniques, the floating random walk was selected and was used to perform our modeling. Even though the calculation of the inner sphere can be more time-consuming than a fixed walk, it means we can avoid the introduction of a bias due to a truncated sphere (represented in green in Fig. 3 (a)). The fixed walk is also suitable for small values of 𝛿 and 𝛿 𝑏 .

Monte-Carlo Method: Complete formulation of the algorithm to the calculation of the flash thermogram

All the different elements required to solve the considered problem and the construction of the coupled walk in conduction and radiation were described previously. This section is dedicated to the description of the general architecture of the algorithm as implemented. Details are given of certain specific elements to make these easier to understand. The aimed configuration to solve is represented in Fig. 4. The solid volume  s is delimited by the surface 𝜕 s , which is broken down by the surfaces 𝜕 fs , 𝜕 rf , 𝜕 𝛷 and 𝜕 p respectively corresponding to the fluid(void)/solid surfaces (diffusive walls), the rear face surface, the excited surface and the adiabatic surfaces (periodicity condition). Periodicity is also applied in the void phase. Hence, each path reaching this virtual wall (represented in blue) is simply specularly reflected.

Initialization

We aim to compute the transient mean surface temperature 𝑇 rf (𝑡) of the rear face, which is defined by:

𝑇 rf (𝑡) = 1 𝑆 ∫ 𝑆 𝑇 (𝐱 obs , 𝑡) d𝑆 ∀𝐱 obs ∈ 𝜕 rf ( 19
)
with 𝑆 the surface of the rear face. The latter expression is rewritten to allow a probabilistic interpretation:

𝑇 rf (𝑡) = 1 𝑙 𝑥 𝑙 𝑦 ∫ 𝑙 𝑥 𝑥=0 ∫ 𝑙 𝑦 𝑦=0 𝑇 (𝐱 obs , 𝑡) d𝑥d𝑦 = ∫ 𝑙 𝑥 𝑥=0 𝑝 𝑋 (𝑥)d𝑥 ∫ 𝑙 𝑦 𝑦=0 𝑝 𝑌 (𝑦)d𝑦 𝑇 (𝐱 obs , 𝑡) ( 20
)
with 𝑙 𝑥 , 𝑙 𝑦 , 𝑝 𝑋 (𝑥) = 1∕𝑙 𝑥 and 𝑝 𝑌 (𝑦) = 1∕𝑙 𝑦 respectively the two lengths of the defined rectangular surface and the two corresponding probability (𝑗) 𝛾 ∈ 𝜕 𝛷 (excited surface) for the 𝑗th times. This path 𝛾, after a given number of jumps, has an associated time 𝑡 -𝜏 (𝑗) 𝛾 and the probability density associated with sampling such a path is 𝑝 (𝑗) 𝛤 (𝐱 obs ). Hence, in the space of the paths 𝐷 𝛤 , the rear face temperature calculation can be summarized as:

𝑇 rf (𝑡) = (𝑡 -𝜏 (1) 𝛾 ≤ 0) × 𝑇 𝐼 + (𝑡 -𝜏 (1) 𝛾 > 0) × ∫ 𝐷 (1)
𝛤 𝑝 (1) 𝛤 (𝐱 (1) 𝛾 )d𝐱 ( 1)

𝛾 [ 𝑇 (𝐱 (1) 𝛾 -𝛿 𝑏 𝐧, 𝑡 -𝜏 (1) 𝛾 ) + 𝛷(𝑡 -𝜏 (1) 𝛾 )𝛿 𝑏 𝜆 s ] (21)
The present algorithm developed for describing the flash method has two important characteristics.

• A path stops only when the initial condition 𝑇 𝐼 (the only known temperature) is reached.

• As described in Section 2.5, the Monte-Carlo weight is incremented when the path reaches the excited face. After this kind of interaction, the walk continues as a conductive walk within the solid phase.

Hence, the Monte-Carlo weight associated with the 𝑖th realization of the algorithm can be simply written as:

𝜃 𝑖 = 𝑇 𝐼 + 𝐽 𝑖 ∑ 𝑗=1 (𝑡 -𝜏 (𝑗) 𝛾 > 0)
𝛷(𝑡 -𝜏 (𝑗) 𝛾 )𝛿 𝑏 𝜆 s ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝑆 𝑖𝑗 [START_REF] Badri | Vectorial finite elements for solving the radiative transfer equation[END_REF] with 𝐽 𝑖 the number of interactions with the excited face of the 𝑖th realization of the MC algorithm. Thus, the MC estimator 𝑇 rf (𝑡) at a given time 𝑡 is:

𝑇 rf (𝑡) ≈ 1 𝑁 mc 𝑁 mc ∑ 𝑖=1 𝜃 𝑖 = 1 𝑁 mc 𝑁 mc ∑ 𝑖=1 ( 𝑇 𝐼 + 𝐽 𝑖 ∑ 𝑗=1 𝑆 𝑖𝑗 ) = 𝑇 𝐼 + 1 𝑁 mc 𝑁 mc ∑ 𝑖=1
𝑆 𝑖 [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm[END_REF] with:

𝑆 𝑖 = 𝐽 𝑖 ∑ 𝑗=1
𝑆 𝑖𝑗 [START_REF] Villefranque | A path-tracing Monte Carlo library for 3-D radiative transfer in highly resolved cloudy atmospheres[END_REF] the source term for the 𝑖th realization of the MC algorithm. In other words, the temperature of the rear face at a time 𝑡 is simply the temperature at the initial condition (t = 0) increased by a given amount due to the energy the flash brings to the front face.

It should be noted that the computation is performed for a given and a unique value of 𝑡. Computing an entire thermogram thus requires time discretization in 𝑁 𝑡 values in the [0, 𝑡 max ] interval and the calculation of the temperature performed for each of these times. Efficient and correct storage of the information generated thus far would mean computing the entire thermogram using only the largest time 𝑡 max could be carried out with great rapidity but the results obtained for different times but computed from the same paths would not be statistically independent. Hence, a Monte-Carlo calculation composed of 𝑁 mc paths for each value of the time interval is preferable and thus was carried out in this work.

Complete formulation of the algorithm

The complete formulation of the algorithm for one realization of the algorithm is fully described in Fig. 5. It summarizes all the previous described steps of the procedure. The interactions with walls (identified by a surface mesh) are added and ensure the transition between conductive and radiative paths.

• The starting point is to sample the starting position on the rear face and to set the starting time (Input, cf. Section 2.7.1). Then, the Monte-Carlo weight and estimator of the source are set at 0.

• The loop starts with a test which consists in determining if the considered position is on a boundary (TB). If this test is negative, the conductive path is performed (C1-C5, cf.Section 2.2). Otherwise, the algorithm determines which boundary is reached (TB1 and TB2). If the path is located on the excited face, the estimator of the source 𝑆 is increased (B1 and B2, cf.Section 2.5).

• If the path is located on a diffusive wall, a Bernoulli test is performed (B4-B5) to determine whether the path follows a conductive path in the solid or a radiative path through the void phase (cf. Section 2.4).

• If the path follows a conductive path, it is simply put back (B3) within the solid phase at a distance 𝛿 𝑏 from the boundary.

Otherwise, the walk continues as a radiative path (R1-R5, cf. Section 2.3)

• The path is stopped when the associated time becomes negative (cf. Section 2.7.2). The resulting Monte-Carlo weight is stored for this realization of the algorithm and a new path is initialized.

Validation on combined conductive and radiative transfers: 3D heterogeneous geometry, benchmark configuration

The purpose of this section is to validate the methodology described above as a direct model for its application in the framework of an inversion procedure to identify equivalent thermophysical properties. The interest of the method is highlighted, a step by step development of the coupled thermal problems to be treated is described and benchmark results are proposed. Firstly, in the case of a flash method configuration, a purely conductive medium is considered. Then, both conductive and radiative transfers are solved on a one-dimensional problem and an effective thermal diffusivity is identified. For the sake of brevity, both validation cases are provided in Appendix. Finally, a verification of the solution by comparing with a commercial software like Comsol ® is provided on a 3D heterogeneous geometry.

To validate the MC algorithm to compute the transient solution in the case of a 3D complex geometry, solutions obtained by the MC method are compared with more common approaches like the deterministic methods used in commercial resolution software. As represented in Fig. 6(a), a heterogeneous 3D honeycomb geometry was considered. This geometry was chosen to be simple enough to be easily reproduced and used as a benchmark for further studies, but also complex enough to be comparable to a porous structure made of solid struts which is commonly studied in engineering science. It is demonstrated in both Appendices A.1 and A.2 that taking a transient heat flux or periodic boundary conditions into account does not cause any particular difficulties. This kind of study can be easily reproduced in future works on the development of the Monte-Carlo method for combined transfers.

The thermal problem is represented in Fig. 6 Dimensionless temperature differences, radiative Biot number and dimensionless temperature are defined as:

𝛥𝑇 * = 𝑇 hot -𝑇 cold 𝑇 ref ; 𝐵𝑖 = ℎ rad 𝑧 max 𝜆 s 𝑇 * = 𝑇 (𝑥 obs , 𝑡) -𝑇 cold 𝑇 hot -𝑇 cold ( 25
)
The dimensionless lengths and time are defined such as 𝑧 * = 𝑧∕𝑧 max , 𝑦 * = 𝑦∕𝑦 max , 𝑥 * = 𝑥∕𝑥 max and 𝑡 * = (𝛼 s 𝑡)∕𝑧 2 max . This problem was solved both by the present Monte Carlo method and by COMSOL multiphysics ® . In the case of the latter, conduction is solved by a finite element method and radiation by the radiosity method, which is well suited to this problem involving surface-to-surface radiative exchanges. Comsol simulation was carried out using 1900 mesh points forming 4563 tetrahedrons and 4008 triangles and automatic time step was used. In the case of the Monte-Carlo method, the associated algorithm to solve such a configuration is easily deduced from the methodology described in Section 2. All paths started at the time and position of interest and then a combined conductive and radiative random walk was performed. A path stopped when a known temperature such as 𝑇 cold , 𝑇 hot or 𝑇 ∞ was reached and this temperature was retained as a realization of the MC algorithm. Then, a new path was initialized. A low dimensionless temperature gradient was deliberately chosen (𝛥𝑇 * = 0.1) to ensure the assumption of linearity of the radiative heat exchanges to be valid.

Transient results obtained with the stochastic and deterministic methods for different positions 𝑧 * = [0.039, 0.27, 0.5, 0.75, 0.96] for 𝑥 * = 𝑦 * = 0.5 (centerline position cf. The dimensionless temperature as a function of the dimensionless time is represented for each radiative Biot number and position in Fig. 8. This figure shows a comparison between the Monte-Carlo method and the deterministic method at the steady state. The dimensionless temperature as a function of the dimensionless position along the centerline is drawn. In all cases, a good level of agreement between the two methods is observed, which validates the use of the stochastic method described to provide a precise solution for a 3D complex geometry and for the thermal problem targeted. The transient variation of the temperature (thermogram) is correctly calculated with the MC algorithm. We observe an increase in the rate of heat transfer (it should be noted that the scale of time axis is significantly different for each sub-graph) with the increase in radiation transfer, i.e. an increase in the total equivalent diffusivity. Moreover, a decrease in the temperature gradients in the volume due to the increase of radiative exchanges is also noticed. These results are consistent with the literature on radiative transfer in semi-transparent media. Fig. 9 shows one of the main interests of the Monte-Carlo method in our present case. In Section 2 we showed that the geometric complexity can be easily handled with tools provided by the computer graphics community (computation of the intersection between a path and the surface mesh) while the complexity of the problem can be easily managed thanks to the formulation of this kind of algorithm (multiple heat transfer modes with a single path). Here, a reduction of the computational time is observed when the role played by radiation increases. For simple problems, the effort to perform a MC simulation is more important than a deterministic one. Nevertheless, MC becomes advantageous with the increase of radiative propagation through the void phase.

This can be explained at the scale of one path. The calculation time is only linked to the capacity of a given path to reach a known temperature or heat sources and thus to end the random walk. Hence, it is linked to the speed of propagation for a given time 𝑡 of the heat transfer between the source and the position of interest (probe computation). Assuming there is pure conduction, a path can only walk through the solid following a Brownian motion until it reaches 𝑇 hot or 𝑇 cold . This could be quite tricky, especially with thin struts, and could require a great deal of computation time. In the case of coupled conductive and radiative heat transfer, the paths can follow a radiative path through the void phase, which allows them to travel along a straight line between two interfaces. The distance traveled in one step of the algorithm is very large, especially through large void channels, compared to the distance traveled during one step of the conductive random walks. Hence, a path can reach 𝑇 hot , 𝑇 cold or even 𝑇 ∞ in a smaller number of iterations of the MC algorithm. Moreover, this phenomenon increases with the increase of the probability to follow a radiative path.

Fig. 9 represents the computing time and the dimensionless computing time according to the position 𝑧 * to perform a given number of realizations of the Monte-Carlo algorithm and for three different radiative Biot numbers on a single CPU Intel Xeon X5650. We note that:

• The computation time for a position at the center 𝑧 * = 0.5 logically corresponds to a maximum. This position is indeed the furthest from the known temperatures.

• The computation times decrease when the radiative transfer increases. This can be explained by the increase of the number of sampled radiative paths as described previously.

• The dimensionless computation times become insensitive to the starting position with the increase of radiative transfers. This can also be explained by the increase of the number of sampled radiative paths. Paths would evolve preferentially in the void phase and reach a known temperature with a low number of steps, whatever the starting point. Fig. 9(b) shows the computing time required to compute the temperature for all the 11 positions for the same radiative Biot numbers in Fig. 9(a) with the MC algorithm. It is compared to the computing time required with the deterministic method (Comsol ® ) to obtain the whole field of temperature. The number of realizations of the MC algorithm was deliberately set to obtain a similar computing time between both methods for 𝐵𝑖 = 7.4 × 10 -3 . As we have explained, it can be noted that the computing time needed with the MC method strongly decreases with the increase of the radiative transfers thus demonstrating better performances in this case. The access of the whole temperature field is not limiting as long as only a probe computation is required in the case of a characterization procedure. This study thus provides advice on choosing a method for a possible user. The deterministic method is preferable for high porosity media and purely conductive transfer while the Monte-Carlo method described was found to give better performances with coupled transfers.

Application to the characterization of porous media

The purpose of this section is to apply the proposed procedure to a numerical thermophysical characterization of heterogeneous structures. Firstly, the construction of the 3D surface mesh and the inversion procedure are detailed. Then, a parametric investigation of the different data of the problem is performed on Kelvin cells, which are commonly used as the geometry of reference in the frame of foams.

Studied geometries and construction

3D foam geometries are constructed using 3D images defined by structured grids of 𝑉 𝑥 ×𝑉 𝑦 ×𝑉 𝑧 voxels, which can be obtained numerically (with a specific software) or experimentally with a scanning method like X-ray tomography. This article deals with the case of purely numerically generated foams with a structured matrix. They are based on a repetition of a unitary tetrakaidecahedron (Kelvin cell) or of a cubicshaped cell as used in the validation section (cf. Section 3. Nevertheless, as highlighted, the proposed method is easily expandable to a wide range of geometries such as randomly distributed matrix foams, fibers, packed spheres or even real X-ray tomographied foams. The initial raw data were obtained from a foam generator software named GenMat ® , which was developed by LTEN (Heat Transfer and Energy Laboratory at Nantes -UMR CNRS 6607) and IUSTI (Institut Universitaire des Systèmes Thermiques et Industriels -UMR CNRS 7343). Based on a watershed marker-based method, foams numerically generated by this software were validated by comparing with a reference SiC foam [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].

GenMat ® provides a stack of grayscale images corresponding to the desired foams. Two stacks of plain material corresponding to the two soleplates (cf. Section 1) were added. Then, the surface mesh was generated with a robust marching-cube method available in 3DSlicer ® (open-source software). Smoothing algorithms (surface preserved) or decimation algorithms (reducing the mesh refinement) can be also used and are available in Meshlab ® and Blender ® (open-source softwares). Finally, a good quality closed mesh was obtained as shown in Fig. 10.

Inversion of the thermogram -equivalent properties of the medium

The transient evolution of the temperature on the rear face is treated exactly as if it were an experimental thermogram. As the thermogram was obtained using a stochastic approach, the temperature calculated at any time is affected by a random error. It should also be noted that the standard deviation associated to each point is easy to obtain, which is one of the strengths of the Monte Carlo methods [START_REF]Monte Carlo Methods[END_REF]. This standard deviation is proportional to the inverse of the square root of the number of Monte Carlo realizations (or of the computation time, this one being proportional to the number of calculated random walks in the medium). From the inversion point of view, this error is similar to the noise observed in an experimental thermogram obtained with a flash method (cf. Appendix A.1).

The inversion of the thermogram requires a direct model of the equivalent homogeneous media like any kind of inverse problem. Here the direct model aims to calculate the thermogram at the rear face of a homogeneous sample sandwiched between two soleplates. This thermogram is indeed the quantity provided by the Monte-Carlo algorithm. The thickness and thermal properties of the two soleplates are known and correspond to the values used in the Monte-Carlo code. The homogeneous medium representing the foam is characterized by its total equivalent thermal conductivity and equivalent specific heat. The direct model used for the inversion procedure is a 1D threelayer purely conductive model. The thermal problem is solved with a quadrupole method [START_REF] Maillet | Thermal Quadrupoles: Solving the Heat Equation Through Integral Transforms[END_REF], which provides a simple semi-analytical solution. The equivalent total thermal diffusivity of the foam can be identified using this inverse method. The total equivalent conductivity can be deduced from the knowledge of the specific heat of the foam and its porosity. A Levenberg-Marquardt algorithm [START_REF] Moré | The Levenberg-Marquardt algorithm: implementation and theory[END_REF] was used to obtain the best match between the thermogram calculated on the actual sample structure (foam struts) with the Monte Carlo algorithm presented above and the thermogram calculated with the equivalent homogeneous medium. More details about inversion procedure can be found in the authors' paper devoted to characterization of foams at high temperature by using the Flash method [START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams at high temperature[END_REF]. The algorithm described in Section 2 and validated herein above was applied to a porous medium composed of a stack of five Kelvin cells between two soleplates as depicted in Fig. 10. Periodic boundary conditions were applied along the 𝑥 and 𝑦 directions such that an infinite medium was actually considered along these two directions (cf. Fig. 4).

The conductivity of the solid phase was fixed at 𝜆 s = 50 W m -1 K -1 (which is close to the thermal conductivity of SiC at 1000 K), the diameter of the cell was 𝐷 cell = 10 mm, the porosity of the foam was 𝜀 = 80%, the investigated range of the reference temperature was between 300 K and 2500 K. The specific heat of the solid was set at 10 6 J m -3 K -1 and the duration of the heat flux flash was set so that 𝑡 𝑑 = 𝑡 max ∕50 (cf. Eq. ( 27)). Let us note that taking constant thermal properties of the strut material along the whole temperature range from ambient to 2500 K is nonphysical. However, the goal of this study was not to describe the variation of the total equivalent conductivity of a real foam but to conduct a parametric study aiming at separating the influence of the different parameters (temperature, porosity, cell size, bulk conductivity and emissivity), that involves rightly to keep all parameters constant except the parameter of interest. Such study is easily feasible numerically whereas it is difficult to be carried out experimentally due to sample availability or due to correlations between parameters and reference temperature.

Fig. 11 shows the different numerical thermograms obtained with the MC algorithm and corresponding best fits after the estimation procedure based on the three-layer purely conductive and homogeneous model. A good level of agreement is observed. This figure shows clearly the acceleration of the thermograms when the reference temperature increases because of the increasing role of radiation. It should be noted that the temperature reached at the steady state is the same for all reference temperatures because the deposited energy was the same. Each thermogram was computed using 2000 Monte-Carlo realizations per time interval and 50 time intervals (thus 10 5 total random walks were computed in total).

The right part and left parts of Fig. 11 show an example of random walk in the structure at a low reference temperature (for which heat transfer is almost purely conductive) and at a high reference temperature (𝑇 ref = 2500 K) (for which radiation plays a significant role), respectively. As previously stated, when the reference temperature increases, the probability of following a radiative path increases. Thus an increasing number of paths (blue lines in insets in Fig. 11) cross the void phase. This travel is quasi instantaneous (due to the high velocity of light), so the mean celerity (i.e. the total equivalent thermal diffusivity) of the heat transfer increases. Hence, the increasing role played by the radiative transfer means the steady state is reached more quickly. This means the starting time of the paths can be decreased which causes a strong decrease in computation time as explained in Section 3. With a single CPU Intel Xeon X5650, a gain of a factor of 50 is observed between a pure conduction calculation and a calculation performed with 𝑇 ref = 2500 K. Hence, the MC algorithm is of interest when radiative heat transfer plays an important role, as was already explained in the previous section. An opposite trend is generally observed following deterministic methods for conduction-radiation algorithms which slow down when the radiation contribution increases due to error control.

For each reference temperature value, fitting the computed thermogram makes it possible to obtain a total equivalent diffusivity of the porous medium. In the following subsections, the influence of the foam's morphological (number of cells, porosity, cell size) and thermophysical (bulk thermal conductivity, emissivity) properties on the equivalent diffusivity are presented. For the sake of brevity, the thermograms will no longer be shown.

Influence of the number of pores in the thickness

The variation of the pore number in the thickness allows us to address the question of the validity of the diffusion approximation to describe the radiative transfer. This assumption is based on the idea that radiative transfer is comparable to a diffusive transfer such as conduction. This means the radiative transfer can be described by the simple definition of a ''radiative'' conductivity 𝜆 rad . In the case of an optically thick participating medium, Rosseland's law [START_REF] Coquard | Conductive and radiative heat transfer in ceramic and metal foams at fire temperatures[END_REF] is commonly used:

𝜆 rad = 16𝜎𝑇 3 ref 3𝛽 eff (26)
The effect of increasing the number of pores in the thickness is that the equivalent optical thickness of the porous medium also increases and thus, the diffusion approximation according to Rosseland becomes valid.

Five foams were generated which were exactly identical in terms of porosity and pore size. These foams differ simply by the number of patterns present in the thickness such as represented in Fig. 12. The conductivity of the solid was decreased to 1 W m -1 K -1 to increase the role of the radiative transfer. For these five foams, the total equivalent estimated conductivity was drawn according to the reference temperature. First, when the reference temperature is low, the equivalent total conductivity is identical between the different structures. This is due to the fact that in pure conduction the equivalent conductivity is independent of the number of patterns. A single cell is sufficient to be representative of the transfers. Secondly, with the increase of the radiative transfer, the total equivalent conductivity increases with the number of pores.

The Rosseland's law given in Eq. ( 26) predicts that the equivalent ''radiative'' conductivity, for a given reference temperature, does not vary with the thickness. Thus, if the medium is optically thick, the transfer by radiation only depends on the equivalent extinction coefficient and the reference temperature. We can observe graphically that the equivalent total conductivity for a foam with 3 cells and 4 cells are very close to that of a foam with 5 cells. This suggests, at least graphically, the validity of the diffusion approximation and a radiative behavior in agreement with the Rosseland prediction for the chosen set of parameters and investigated range of reference temperature. The proposed tool enables the quantification of an equivalent radiative conductivity, an exploration of the validity of the diffusion approximation and a simplification in the modeling. In the following, only foams with at least 4 cells along the thickness are considered.

The influence of the porosity

In this section the influence of the porosity and shape of the solid matrix on the coupled heat transfer is presented. Two kinds of cell geometries were considered: Kelvin cell and honeycomb cell (cubicshaped cell) as represented in Fig. 13. Cell size is set to 2 mm. A stack of 4 cells along the heat transfer axis was considered. For each structure, the porosity value of both structures varied from 0.2 to 0.8. The evolution of the total equivalent conductivity as a function of the reference temperature is drawn.

For a given porosity, a logical linear increase in the total equivalent conductivity with the third power of the reference temperature is observed (cf. Rosseland relation Eq. [START_REF] Ibarrart | Combined conductive-convective-radiative heat transfer in complex geometry using the Monte Carlo method : Application to solar receivers[END_REF] or according to other laws dedicated to radiative conductivity [START_REF] Poltz | The thermal conductivity of liquids-IV: Temperature dependence of thermal conductivity[END_REF]). At a low temperature, heat transfer arises from conduction only and a very well-known result is observed [START_REF] Ranut | On the effective thermal conductivity of aluminum metal foams: Review and improvement of the available empirical and analytical models[END_REF], i.e. the higher the porosity, the smaller the conductivity. Conversely at high temperature, the higher the porosity, the higher the conductivity is because the emptier the structure, the easier radiation can take place and transfer heat. We notice that the conductivity depends little on the porosity at a given temperature near 1200 K (4𝜎𝑇 3 = 400 W m -2 K -1 ) in this example, except for the lowest porosities. Moreover, the influence of the shape of the matrix is highlighted. For the same set of parameters, the total effective conductivity of the cubicshaped cell is greater than the Kelvin cell. This leads to the idea that a cubic shaped foam is more transparent than a Kelvin cell foam, i.e. it has a lower effective extinction coefficient. The slope of the curve presented in Fig. 13 means a radiative coefficient equivalent to the porous medium could be extracted following a judicious choice of an equivalent radiative law.

These results seem logical and expected. The proposed procedure involves retrieving the classic results from the literature to explore the influence of parameters like porosity and carry out a better quantification of the major heat transfer mechanism for a given structure. In the case of pure conduction, the estimated conductivity is compared to the one proposed by Pabst et al. [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] for a Kelvin cell with the same morphological parameters. As presented in Fig. 14, a good level of agreement is found which meant the result could be validated with a high level of confidence in the proposed method.

The influence of cell size

This section focuses on the influence of cell size. The geometry studied was composed of four Kelvin cells with a fixed porosity equal to 80%. We varied the cell size between 1 mm and 10 mm. The equivalent conductivity as a function of the third power of reference temperature is represented in Fig. 15. In pure conduction (low temperature), the equivalent conductivity does not depend on cell size. The solid matrix is organized in exactly the same way and the amount of material (porosity) is the same which means there is no reason to observe an influence of cell size. This is a well-known result in foam literature [START_REF] Ranut | On the effective thermal conductivity of aluminum metal foams: Review and improvement of the available empirical and analytical models[END_REF]. When the temperature increases, radiation transfers increase too and we observe that the larger the cell size, the higher the total equivalent conductivity of the homogeneous medium. As noted in the previous section, this result is logical. The mean free radiative path is indeed greater for larger cells, which corresponds to a lower extinction coefficient. It should be noted that, in theory, the equivalent optical thickness is the same for all these different structures. However, as the sample thickness increases when the size of the cell increases, this corresponds to a decrease in the extinction coefficient 𝛽 (the optical thickness being defined by 𝜏 = 𝛽𝑒) and thus an increase in radiative conductivity as predicted, for example, by the Rosseland model.

Influence of the bulk conductivity of struts

In Fig. 16, the results obtained for a stack of 4 Kelvin cells with a porosity equal to 0.8 and bulk conductivity of the struts varying from 0.1 W m -1 K -1 to 100 W m -1 K -1 are described. This study illustrates how conduction and radiation through the porous medium are in a competition to become the major heat transfer mechanism. The ratio 𝜆 eq ∕𝜆 𝑠 is presented in Fig. 16(a). If the temperature is low enough for heat transfer to be due to conduction alone, then the equivalent conductivity is proportional to the bulk conductivity and consequently the ratio 𝜆 eq ∕𝜆 𝑠 is logically independent of the thermal conductivity of the solid. This ratio only depends on structural properties such as porosity and the shape of the solid matrix. When the temperature increases, the ratio is no longer independent on the bulk conductivity because of the increasing role of radiation transfer, which becomes quickly the main transport mode when the thermal conductivity is low.

Fig. 16(b) presents the ratio of the radiative thermal conductivity by the total equivalent conductivity (percentage) as a function of the reference temperature for different conductivities of the bulk is represented. For a low bulk conductivity (e.g. 𝜆 𝑠 = 0.1 W m -1 K -1 ), radiation transfer very quickly becomes the main transfer method, even for low reference temperatures. For example, for 𝑇 ref = 500 K, radiation accounts for already almost 60% of the total transfer. On the other hand, for 𝜆 𝑠 = 100 W m -1 K -1 , conduction transfer remains predominant over radiation transfer even at high temperatures.

These observations can be understood through the scope of the stochastic method developed in this paper. Indeed, the competition between the two heat transfer modes appears because of the probability of following a radiative path 𝑃 rad . If the conduction term 𝜆 𝑠 ∕𝛿 𝑏 is very large compared to the radiative term ℎ rad , then the problem is purely conductive and it is weakly sensitive to radiative transfer and thus to temperature. On the other hand, considering an extreme case in which the radiation is very large in comparison to the conduction, the total heat transfer does not depend on the conductivity of the solid.

The influence of emissivity

In this section, the influence of the emissivity of struts on the total equivalent heat transfer is studied. In all previous cases, the emission of radiation was hemispherical (lambertian or isotropic on the hemisphere) and the emissivity was chosen equal to 1. On a structure composed of a stack of five Kelvin cells along the direction of heat transfer, simulations were performed for an emissivity and a reference temperature varying between 0 to 1 and 300 K to 2500 K, respectively. Fig. 17 shows the ratio 𝜆 eq ∕𝜆 𝑠 as a function of the reference temperature. In this work, perfectly diffuse reflection was implemented. Nevertheless, it is possible to specular reflection or combined specular/diffuse reflection into account using the proposed stochastic When transfers by radiation are weak, the emissivity logically has no influence on the equivalent total conductivity because conduction is the main mode of transfer. Then, with the increase of the reference temperature (i.e. radiation), an increase of the equivalent conductivity with the emissivity and, logically, with the temperature considered, is observed. For the asymptotic case 𝜀 = 0, no radiative transfer can occur. Then, a rapid increase in total equivalent conductivity between 𝜀 = 0 and 𝜀 = 0.25 is observed and reaches its maximum for 𝜀 = 1. In the case of an optically thick medium, this kind of phenomenon cannot be explained by the influence of the emission of the intern part of the two soleplates alone. An effective optical thickness of the porous medium, which can be seen as the statistical mean free paths traveled within the fluid phase and thus depends only on the morphological properties of the foam, appears here to depend on the emissivity. Evidence on this kind of result was also given by Vignoles [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], Patel et al. [START_REF] Patel | Development of correlations for effective thermal conductivity of a tetrakaidecahedra structure in presence of combined conduction and radiation heat transfer[END_REF] or Li and Wang [START_REF] Li | Equivalent thermal conductivity of open-cell ceramic foams at high temperatures[END_REF]. The latter also showed that the specular/diffuse mode of reflection has a major role on the effective extinction coefficient. This aspect could usefully be explored in more detail in further studies.

Conclusion

In this paper, a new tool for modeling and characterizing coupled transfers in complex geometric media such as porous media was proposed. The recent developments of the Monte-Carlo method make it possible to imagine, in the framework of linear thermal transfers, to solve coupled and multiphysics transfers with a single algorithm. This method is very interesting for the characterization of complex media for two reasons. Firstly its capacity to handle the increasing complexity of the solved physics (directional dependency of emission, semi-transparent medium, frequency dependence, anisotropy of the solid phase, etc.) without an increase in computation time. Secondly the calculation performances allowed by using tools from the graphic community (insensitivity of the refinement of the mesh, parallelization, etc.) make it an effective option as well. Moreover, the probe computation aspect, the statistical interpretation of the space of paths and the approach's capacity to simulate the combined heat propagation are very appropriate in the context of study and understanding of coupled mechanisms. For the first time, a numerical reproduction of the flash method in the case of coupled conductive and radiative transfers (opaque surfaces) was performed.

Firstly, the methodology was fully described. This highlighted the proposed formulation's various positive aspects such as its practical implementation and the understanding it provides of the intrinsic thermophysical properties through probabilities. However, the method requires the introduction of arbitrary lengths to generate the walk. Future studies will focus precisely on determining a time/accuracy trade-off for these parameters.

Next the proposed MC algorithm and its use as a numerical measuring instrument were validated by comparison with deterministic solvers such as FlexPde ® , Comsol ® or a semi-analytical solution provided by the quadrupole method. A validation of the transient solution was provided on a 3D heterogeneous benchmark case. The results showed a very good level of agreement with deterministic solvers and highlighted the interest of using the present method in the case of an increasing influence of the radiative transfer. New transport modes such as advection, conduction in the fluid phase and participating solid phase can be added.

Finally, this code was used to reproduce the usual flash method applied to experimentally determine the diffusivity of foams, as was recently performed by the authors. The thermal response of a foam sample sandwiched between two soleplates and submitted to a heat radiative flux flash was computed. The thermogram obtained numerically were fitted against the thermogram obtained with a direct model involving a homogeneous medium between the same soleplates. A total equivalent diffusivity (conduction+radiation) for the foam was thus obtained. The proposed procedure and its coupling with recently developed tools such as GenMat ® (foam generator) enabled us to propose a parametric study of how morphological and thermophysical parameters influence the equivalent total conductivity of the foam. Our results showed a good level of agreement with the associated literature. The proposed method can be adapted to new geometries and to new characterization configuration. Further studies could focus on new heterogeneous structures such as X-ray tomographied foams, numerical foams, packed sphere or fibers and new characterization methods such as the hot-guarded plate or the hot-wire method used to study combined heat transfer. Hence, correlations of effective equivalent properties of such media could be provided in parallel of experimental developments. can be taken into account during the inversion procedure (like during an experimental measurement) and implies an error in the estimated parameters. It is therefore possible to reduce computation time with the present approach but with reduced precision in the estimated parameters.

A.2. Coupled conductive and radiative transfers: 1D flash configuration, opaque parallel plates separated by void phase

The previous example only dealt with the simple case of pure conduction in a homogeneous medium. The more complex case of coupling between conduction and radiation is considered here. This academic case was studied by Vignoles [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] to validate his coupled conductive and radiative stochastic method and to identify the equivalent thermal diffusivity of foams [START_REF] Vignoles | Numerical study of effective heat conductivities of foams by coupled conduction and radiation[END_REF]. The problem is defined by solid plates of thickness 𝐿 1 which are homogeneous, parallel and opaque with a given emissivity 𝜖. Each plate is separated by a void phase of thickness 𝐿 tot -𝐿 1 and are assumed infinite in the other two directions. Thus, the considered thermal problem is unidirectional.

This thermal problem has an analytical solution. By writing the balance of the conductive and radiative fluxes (cf. Eq. ( 9)) at the interface, the equivalent diffusivity between two separated plates can be written as:

𝛼 eff = 𝛼 s (1 -𝛱) 2 ( 1 + 2 -𝜖 𝐵𝑖 )-1 (28)
with 𝛼 eff the equivalent diffusivity, 𝛱 = 1-𝐿 1 ∕𝐿 tot the volume fraction and 𝐵𝑖 the radiative Biot number defined such as:

𝐵𝑖 = ℎ rad 𝐿 1 𝜆 s = 4𝜖𝜎𝑇 3 ref 𝐿 1 𝜆 s ( 29
)
The simulation with the MC algorithm was performed on the geometry represented in Fig. 19(a). The simulation of the coupled heat transfer of the equivalent homogeneous medium was based on the quadrupole method [START_REF] Maillet | Thermal Quadrupoles: Solving the Heat Equation Through Integral Transforms[END_REF], which is very convenient for simulation of multi-layer and one dimensional heat transfer. A bi-layer medium was thus modeled. The first layer was defined with a thermal diffusivity 𝛼 eff , which accounts for the first soleplate and the void phase while the second, defined with 𝛼 s , accounts only for the second soleplate. The inverse procedure aims to identify the value 𝛼 eff of the homogeneous model that allows the best match with the thermogram obtained with the 3D simulation computed with MC. In Fig. 19(b), the estimated values of the equivalent diffusivity are represented as a function of the radiative Biot number for a given value of volume fraction 𝛱. These values and their associated uncertainties are compared with the analytical solution given in Eq. [START_REF] Tregan | Transient conducto-radiative heat transfer in a single Monte-Carlo algorithm: Handling the nonlinearity[END_REF]. For 𝐵𝑖 > 0.1, a good level of agreement was found between our estimated equivalent diffusivity and the analytical solution. Given the associated uncertainty, the results are statistically compatible. This means we can have great confidence in the described procedure involving direct modeling with combined MC algorithm and inversion based on the use of a homogeneous equivalent model.

Computation with the MC algorithm becomes difficult in the case of a low radiative Biot number (little body assumption 𝐵𝑖 < 0.1) which constitutes a limit of the present MC algorithm. For low value of ℎ rad , the probability of a path located at the boundary to jump through the void phase is very low and therefore many attempts are required to allow a given path to reach the front face, location of thermal excitation. This means the computation time needed to build a usable thermogram is increased too much. Convergence issues are highlighted even in this simple 1D case. Solutions exist like the introduction of arbitrary probabilities but have not been considered in the case of the described MC algorithm and could be the subject of further studies.

𝐻

  Cumulative density function of the time (-)𝑝 𝑋probability density function of a random variable 𝑋 (-)

Fig. 1 .

 1 Fig. 1. Principle of the flash method and walking technique for conduction.

Fig. 2 .

 2 Fig. 2. Conductive-radiative coupling, case of opaque surfaces: ray-tracing within the fluid phase (vacuum).

Fig. 3 .

 3 Fig. 3. Random walks techniques : (a) Fixed random walk (b) Floating random walk.

Fig. 4 .

 4 Fig. 4. Detailed description of the computation of one realization of the transient conductive-radiative MC algorithm : application to the flash thermogram.

  (b). 𝑧 max and 𝑥 max = 𝑦 max are the dimensions of the structure. Square channels with side 2 × M.Sans et al.

Fig. 5 .

 5 Fig. 5. Detailed description of the computation of one realization of the transient conductive-radiative MC algorithm : application to the flash thermogram.

Fig. 6 .

 6 Fig. 6. Studied benchmark case.𝐷 s define the void phase while square channels with side 𝐷 s represent the struts (solid phase). The solid is assumed to be opaque with a unit emissivity. As described in Section 2, conduction only occurs in the solid phase while radiative transfer only occurs in the void phase.As represented in Fig.6(b), the thermal boundary conditions are 𝑇 hot , the temperature imposed on the upper boundary, 𝑇 cold , the temperature imposed on the lower boundary and 𝑇 ∞ the ambient radiative temperature (infinitely, far from the geometry). The structure is under vacuum so that no natural convection can occur. Surfaces can exchange with the environment through radiation. The reference temperature and the ambient temperature are set so that 𝑇 ref = 𝑇 ∞ = (𝑇 hot +𝑇 cold )∕2, the average between the two imposed temperatures. At the initial time, the entire structure is in equilibrium at 𝑇 cold .

Fig. 7 .

 7 Fig. 7. Transient temperature along the 𝑧 centerline for different values of Biot number and for low temperature gradient 𝛥𝑇 * = 0.1.

  Fig. 6(b)) are drawn in Fig. 7. Each sub-figure corresponds to: • Fig. 7(a), 𝐵𝑖 = 7.4 10 -3 conduction is dominant • Fig. 7(b), 𝐵𝑖 = 1.6 conduction and radiation have a similar influence • Fig. 7(c), 𝐵𝑖 = 102 radiation is dominant

M. Sans et al.Fig. 8 .

 8 Fig. 8. Dimensionless temperature as a function of the dimensionless position for different radiative Biot numbers at steady-state.

Fig. 9 .

 9 Fig. 9. Computing time.

Fig. 10 .

 10 Fig. 10. Example of the studied geometry: four Kelvin cell sandwiched between two soleplates.

4. 3 .

 3 Parametric study and discussions 4.3.1. Influence of the reference temperature

Fig. 11 .

 11 Fig. 11. Numerical thermograms (symbols) and corresponding thermogram fits (plain lines) obtained with the three-layer purely conductive and homogeneous model (cp). Right and left insets: representation of random walk paths (blue lines) inside the structure. Right inset: 300 K temperature, radiation is negligible most of the random walk stays inside the struts. Left inset: 2500 K temperature, the random walk can jump straightly from one strut to another strut.

M

  .Sans et al.

Fig. 12 .

 12 Fig. 12. Influence of the number of pores in the thickness -𝜆 s = 1 W m -1 K -1 , 𝜀 = 80%, 𝐷 cell = 10 mm.

Fig. 13 .

 13 Fig. 13. Influence of the porosity on the total heat transfer -𝜆 s = 1 W m -1 K -1 , 𝐷 cell = 2 mm.

Fig. 14 .

 14 Fig. 14. Equivalent thermal conductivity of Kelvin cell according to porosity: pure conduction.

Fig. 15 .

 15 Fig. 15. Influence of cell size on the total heat transfer: Kelvin cell -𝜆 s = 1 W m -1 K -1 , 𝜀 = 80%.

Fig. 16 .

 16 Fig. 16. Influence of the bulk conductivity for a structure composed of 4 Kelvin cells -𝜀 = 70%, 𝐷 cell = 2 mm.

Fig. 17 .

 17 Fig. 17. Influence of the emissivity of struts -𝜆s = 50 W m -1 K -1 , 𝜀 = 70%, 𝐷 cell = 2 mm.

Fig. 18 .

 18 Fig. 18.Comparison of the rear face thermograms on a Flash configuration computed with the Monte Carlo method and the finite element method (commercial software).

Fig. 19 .

 19 Fig. 19. Comparison of the analytical solution and the solution estimated by inversion of the thermogram generated by MC.

𝑇 𝐼 (initial condition) else 𝑇 (𝐱, 𝑡) = 𝑇 (𝐱 + 𝛿𝐮, 𝑡 -𝜏) (recursivity) end if

	1 Sampling of a conductive path (infinite medium) Sample a direction of travel 𝐮 according to 𝑝 𝐔 Sample a time 𝜏 according to 𝑝 𝑇 if 𝑡 ≤ 𝜏 then 𝑇 (𝐱, 𝑡) =

  Sampling of a spatially averaged temperature: Initialization Sample location 𝑥 according to 𝑝 𝑋 (𝑥) Sample location 𝑦 according to 𝑝 𝑌 (𝑦) Set position 𝐱 obs 𝑇 rf (𝑡) = 𝑇 (𝐱 obs , 𝑡) (recursivity) 2.7.2. Computation of the Monte-Carlo weight All paths start at the time of interest 𝑡 and at the position 𝐱 obs sampled uniformly on the rear face (cf. Section 2.7.1). Let us define 𝛾, a given path that reaches the location 𝐱

M.

Sans et al.

densities. The initialization can then be easily performed following algorithm 4: Algorithm 4

Acknowledgments

The authors gratefully acknowledge Benoit Rousseau (LTEN -UMR CNRS 6607) for providing and support regarding the use of GenMat. We also wish to acknowledge Cyril Caliot (LMAP -UMR CNRS 5142), Mouna El Hafi (RAPSODEE -UMR CNRS 5302), Richard Fournier and Stephane Blanco (LAPLACE -UMR CNRS 5213) for their motivating advice and the fruitful discussions. We would like to thank the Meso-Star start-up and the EDStar platform for making their libraries Star-Engine available to us and for the support provided. High Performance Computing resources were partially provided by the EXPLOR centre hosted by the Université de Lorraine

Data availability

No data was used for the research described in the article.

Declaration of 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.

Appendix

A.1. Purely conductive case: Flash configuration

For this first case, a purely conductive homogeneous solid medium (no radiation) was considered. A transient heat flux was applied on the front face such that:

with 𝑡 𝑝 the time duration of the excitation and 𝑄 the constant heat flux deposited. The parameters were set as 𝜆 s = 1 W m -1 K -1 , 𝛼 = 1 mm 2 s -1 and 𝑡 𝑝 = 10 s. Studied geometries were simply a square and a cube whose sides had a length of 10 cm. Fig. 18 shows thermograms computed with the described MC method and with the finite elements method (commercial software Flexpde ® ). A dimensionless temperature such as 𝑇 * = (𝑇 rf (𝑡) -𝑇 𝐼 )∕(𝑇 max -𝑇 𝐼 ) was introduced with 𝑇 max the maximum temperature of the thermogram. In Fig. 18(a), on a square 2D geometry, the thermogram was discretized in 1000 time intervals and 10 5 realizations of the MC algorithm were performed for each time value. Convective losses on the front and rear face were considered (ℎ = 5 W m -2 K) with a fluid temperature set to the initial temperature 𝑇 𝐼 . The other boundaries were adiabatic. Thus, the heat transfer was unidirectional. A good level of agreement was observed between the two methods and the well-known shape of the rear face thermogram for the flash configuration was successfully obtained. In Fig. 18(b), on a cube 3D geometry, the thermogram was also discretized in 1000 time intervals but only 5000 realizations of the MC algorithm were performed. Convective losses were removed so the system reached a new temperature equilibrium at steady-state. In that case, a good level of agreement was also found between the two methods. Nevertheless, the solution computed with the Monte-Carlo method appears noisier than in the previous case due to a low number of realizations of the MC algorithm per time.

We would like to simply stress at this point the fact that the uncertainty associated with a MC computation is dependent on the number of realizations performed. The solution computed by the Monte Carlo algorithm and its associated error can be compared to the thermogram obtained by a measuring instrument. For the Monte Carlo algorithm, the associated uncertainty simply depends on the number of paths followed and, therefore, on the calculation time. The solution obtained by Monte Carlo tends towards the solution obtained by the deterministic solver with the increase of the number of realizations. This noise is independent and centered on the solution which means it