Loïc Seyer
email: seyer.loic@gmail.com

Vincent Gonneau

Franck Enguehard

Denis Rochais

Improvements of the Brownian walkers method towards the modeling of conductionradiation coupling

Keywords: transient heat transfer, heat conduction, conduction-radiation coupling, Brownian walkers, heterogeneous media, voxelized representation

This work directly follows the one presented in a previous article of V. Gonneau, D.

methods for the modeling of transient thermal conduction within a heterogeneous medium by the movement of Brownian walkers. The material structure is voxelized, and each walker transports an elementary enthalpy during its displacement within the structure. This enthalpy transport associated with the displacement of the walkers represents the conductive flow and makes it possible to simulate transient conduction with a quantitative stochastic approach. This article presents a new method for accounting for Dirichlet type boundary conditions that is quite efficient and that allows to relax quite substantially a constraint of maximum value of the time step of the Brownian walker simulations that had previously been pointed out in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. Other boundary condition treatments are also upgraded, in terms of speed (for Neumann type conditions) or generalization (for Robin type conditions). Then, a first non-linear conduction-radiation coupling model in an optically thick and grey semi-transparent medium is solved by Brownian walkers through two approaches: the first one based on a global conductivity function of temperature and allowing to confirm the validity of a stochastic transmission criterion at the interface of two constituents within a heterogeneous medium, and the second one based on the radiative volume power field. This second approach is preceded by the description and the validation of a procedure for the management of a volume power field using Brownian walkers. This procedure introduces the notion of negative Brownian walkers, which prove necessary for representing negative local values of the volume power field as is frequently the case with the radiative volume power.

Nomenclature

𝑎

Thermal diffusivity (m 2 .s -1 ) 𝛿𝑡 Time step (s)

𝐵𝑖

Biot number (-) 𝛿𝑥 1D space discretization (m)

Introduction

In the past few decades, heterogeneous media with complex 3D morphologies like composites, cellular materials, or ceramics have proved their excellent mechanical and thermal properties. They are used in many industrial applications like gas diffusion layers [START_REF] Ramousse | Estimation of the effective thermal conductivity of carbon felts used as PEMFC Gas Diffusion Layers[END_REF], solar energy volumetric receivers [START_REF] Kribus | The promise and challenge of solar volumetric absorbers[END_REF], heat exchangers [START_REF] Ortona | SiSiC heat exchangers for recuperative gas burners with highly structured surface elements[END_REF] or as thermal insulators at high temperatures (>1800 K) [START_REF] Ferrari | Sandwich structured ceramic matrix composites with periodic cellular ceramic cores: an active cooled thermal protection for space vehicles[END_REF]. To improve their properties or even extrapolate them to better performing new material structures, their thermal behavior has to be always better and more finely understood. To this end, dedicated experimental characterization setups are developed, and in parallel, numerical modeling approaches of varied inspirations taking account of the morphological descriptions of the materials are proposed to simulate the experiments.

Historically, the most used numerical methods are the deterministic ones, because they are usually very efficient and in most cases quite fast in terms of calculation times. However, these techniques reach their limits in the case of complex morphologies. In the context of our work, we aim to calculate quantitative transient temperature fields within heterogeneous semi-transparent media of complex morphologies submitted to high thermal loads, thus reaching high temperatures. In this situation, the radiative transfer can be predominant and must be accounted for. Deterministic approaches for the treatment of volume radiative transfer (finite differences, finite elements or finite volumes coupled to discrete ordinates or spherical harmonics) solve the Radiative Transfer Equation and hence quantify the radiative intensity field as a function of position and direction (and wavelength when the semi-transparent medium is not grey). These methods require a preliminary step of spatial discretization of the calculation domain, which is quite computationally intensive in the case of complex material morphologies [START_REF] Badri | High performance computation of radiative transfer equation using the finite element method[END_REF][START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF]. Consequently, following the tremendous increase in computing power of the last decades, new methods based on random events simulating the physical phenomena of interest and called Monte Carlo methods have emerged in the 1970s and are now considered as the reference methods in particular when the problem to solve presents a large number of degrees of freedom. Considering conduction and radiation heat transfer, these methods are based on enthalpy carriers (respectively Brownian walkers for conduction and photon bundles for radiation) the movement of which is ruled by specific statistical laws accounting for the physical phenomena acting on the carriers. Monte Carlo approaches present two important advantages: they are compatible with voxelized descriptions of heterogenous material structures (for example issued from X ray tomographies), hence avoiding the complex mesh generation procedure, and they are quite well adapted to radiative modeling.

Many studies have been performed to characterize conductive or radiative properties in the case of conduction-radiation coupling with deterministic approaches, e.g. with the Pn method [START_REF] Ferkl | Multiphase approach to coupled conduction-radiation heat transfer in reconstructed polymeric foams[END_REF], the two flux method [START_REF] Niezgoda | Modeling heat transfer within porous multiconstituent materials[END_REF], or the discrete ordinate method [START_REF] Hardy | 3D numerical modeling of the propagation of radiative intensity through a X-ray tomographied ligament[END_REF][START_REF] Boulet | Model of radiative transfer in fibrous media[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]. Recently, new developments based on finite elements [START_REF] Ouchtout | Finite element framework for modeling conductoradiative transfers within heterogeneous media at both discrete and continuous scales[END_REF] and finite volumes [START_REF] Kumar | A modified zonal method to solve coupled conduction-radiation physics within highly porous large scale digitized cellular porous materials[END_REF] have been proposed, stimulated by the increasing power of calculators. Some studies put forward hybrid methods, i.e. methods combining deterministic and stochastic approaches: a finite volume technique for conduction is coupled to a radiation Monte Carlo method in [START_REF] Haussener | Tomography-based heat and mass transfer characterization of reticulate porous ceramics for high-temperature processing[END_REF]. Recently, some stochastic models using a grid or a mesh for conduction have been proposed such as the Lattice Monte Carlo method [START_REF] Belova | The Lattice Monte Carlo method for solving phenomenological mass and heat transport problems[END_REF] or the Off-Lattice Monte Carlo method [START_REF] Gong | Review of recent developments on using an Off-Lattice Monte Carlo approach to predict the effective thermal conductivity of composite systems with complex structures[END_REF].

Even more recently, the so-called Stardis project [START_REF] Penazzi | Stardis: propagator evaluation for coupled heat transfer in large geometric models[END_REF] proposed an entirely stochastic conductionradiation coupling using Brownian walkers and the Monte Carlo method. Some of these authors proposed a review of fully Monte Carlo methods for linear situations [START_REF] Tregan | Coupling radiative, conductive and convective heat transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations[END_REF]. Yet this method is developed for mathematical 3D models (thus meshless and non-voxelized structures), and only evaluates the temperature at a selected probe point in space and time. Vignoles et al. [START_REF] Vignoles | A Brownian motion algorithm for twoscale modeling of chemical vapor infiltration[END_REF] developed a stochastic Brownian walker method for conduction with the advantage of being based on voxelized structures. Following this work, a fully stochastic approach adapted to voxelized structures has been proposed by Vignoles 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] using a Monte Carlo method for radiation and the Brownian walker method for conduction. In this work, the heterogeneous medium is supposed to be made of two constituents, one opaque and the other one transparent, and the enthalpy carrier is either a Brownian walker in the opaque phase or a photon bundle in the transparent phase. This method presents the great advantages to be entirely stochastic and usable with voxelized structures. On the other hand, it is limited to situations where the radiative transfer can be linearized, and the conduction-radiation coupling is only partial: each phase of the material is subject to only one heat transfer mode, and the coupling appears only at the interfaces between the two phases. In order to extend the field of application of all-stochastic conduction-radiation coupling approaches, Dauvois [START_REF] Dauvois | Modélisation du transfert thermique couplé conductif et radiatif au sein de milieux fibreux portés à haute température[END_REF] upgraded the radiation Monte Carlo method to non-Beerian media and made a first attempt to couple Brownian walkers and photon bundles within a binary structure made of two semi-transparent constituents. This preliminary work has revealed a lack of mastery of the Brownian walker kinetics when dealing with heterogeneous media (how do walkers behave at interfaces?) or when this kinetics must account for prescribed thermal boundary conditions (how must walkers behave at a boundary to account for a

Dirichlet condition for example?

). To answer these questions, Gonneau [START_REF] Gonneau | Modélisation du transfert thermique par marcheurs Browniens dans des milieux hétérogènes[END_REF] focused on the Brownian walker method alone and proposed stochastic behavior laws for the walkers at interfaces and boundaries in order the Brownian walker simulations to fulfill the classical temperature and heat flux conservation relations at interfaces on the one hand and to represent classical (Dirichlet, Neumann and Robin) boundary conditions on the other hand. The next step is now to combine the recent advances concerning Brownian walker kinetics and the more common radiative Monte Carlo method and to revisit the all stochastic method first proposed by Vignoles 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] but with a full conduction-radiation coupling and no need for any radiative transfer linearization.

This work revisits, improves and validates certain theoretical concepts concerning Brownian walker kinetics that have been introduced in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. Section 2 is devoted to the presentation of improvements in the treatment of Dirichlet, Neumann and Robin boundary conditions. It is shown in particular that our new treatment of Dirichlet boundary conditions allows to relax quite substantially a constraint of maximum value of the time step of the Brownian walker simulations pointed out in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. Section 3 presents our first conduction-radiation coupling activities within an optically thick and grey semi-transparent medium through two approaches, the first one based on the notion of temperature dependent radiative thermal conductivity, and the second one involving the radiative volume power field. This second approach requires to insert a volume power field in the transient energy conservation equation, the values of this field being potentially negative. This situation will lead us to introduce the notion of negative Brownian walkers and to specify their kinetics and how to account for them.

Improvement of methods accounting for classical boundary conditions with Brownian walkers

In [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], we have detailed how to model a Dirichlet condition (imposed temperature), a Neumann condition (imposed heat flux) and a Robin condition (convective heat transfer at a wall) with Brownian walkers. This section will describe how to accelerate calculations in the presence of Dirichlet and Neumann conditions, and how to extend the Robin condition to convective heating.

Moreover, it will be shown that our improved treatment of the Dirichlet condition substantially relaxes a strong constraint over the time step of the Brownian walker simulations.

As a reminder of the results of [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], when simulating transient conductive transfer with Brownian walkers, the accuracy of the results is strongly related to the value of the dimensionless ratio 𝛼 defined as:

𝛼 = √2𝑎𝛿𝑡 𝛿𝑥 ( 1
)
where 𝑎 is the thermal diffusivity of the material, 𝛿𝑡 the time step of the Brownian walker simulation and 𝛿𝑥 the spatial discretization step. The parameter 𝛼 represents the ratio of the characteristic diffusion length over the duration 𝛿𝑡 to the length 𝛿𝑥. It represents the order of magnitude of the number of voxels that a walker crosses during one time step. When working on the voxelized representation of a 3D material structure (for example issued from a X ray tomography), the spatial discretization step 𝛿𝑥 is imposed by the resolution of the tomography equipment, so that imposing a value of 𝛼 is equivalent to imposing a value of the time step 𝛿𝑡. Furthermore, when simulating a transient 1D Cartesian conduction problem using a finite difference algorithm and an explicit scheme, a well-known result is that for a given spatial discretization step 𝛿𝑥, the algorithm is stable if the time step of the calculation 𝛿𝑡 fulfills the so-called Courant-Friedrichs-Lewy (CFL) condition

𝑎𝛿𝑡 𝛿𝑥 2 ≤ 1 2 .
This condition imposes a maximum value to 𝛿𝑡 (equal to 𝛿𝑥 2 2𝑎 ⁄ ), and the CFL condition, when reexpressed as a function of the dimensionless parameter 𝛼 introduced in (1), becomes 𝛼 ≤ 1.

After examining the results obtained by Brownian walkers for several transient 1D Cartesian conduction problems involving different boundary conditions, we had come to the conclusion that 𝛼 must be close to 2 in order to obtain a good tradeoff between accuracy of the results (𝛼 sufficiently small, but not too small in order to avoid artificial acceleration of the heat transfer) and rapidity of the calculations (𝛼 large enough, but not too large in order to prevent temperature jumps at the boundaries) [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. Consequently, the maximum value of the time step of Brownian walker calculations is 4 times the one imposed by the CFL condition. On the other hand, this condition 𝛼 ≈ 2 remains an important constraint because it generally imposes to have small time steps, which in turn imposes numerous time iterations when the asymptotic steady-state temperature field is searched for. Given that the Brownian walker method is transient by definition, this constraint usually leads to long computation times. In the following, an improved procedure for managing Dirichlet boundary conditions, that allows to relax the constraint 𝛼 ≈ 2 quite substantially, is introduced.

Improvement of the Neumann boundary condition management

Before dealing with the Dirichlet boundary condition, the treatment of Neumann boundary conditions with Brownian walkers presented in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] is revisited. This section has two purposes: (i) first, to present a new method of management of the Neumann boundary condition which proves to be much faster in terms of computation time than the previous one, (ii) second, to show that in the case of a problem involving only Neumann boundary conditions, the constraint 𝛼 ≈ 2 is actually inexistent.

To represent the heating by an imposed heat flux, we start from the method described in Section 2.5.2. of [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], i.e. by injecting, between instants 𝑡 and (𝑡 + 𝛿𝑡), a number of Brownian walkers per unit surface 𝛿𝑀 𝑠 (𝑡) given by:

𝛿𝑀 𝑠 (𝑡) = 𝜑(𝑡)𝛿𝑡 𝛿ℎ 𝑟 ( 2
)
where 𝜑(𝑡) is the heat flux imposed at the boundary and 𝛿ℎ 𝑟 the elementary enthalpy carried by each walker. In [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], we mentioned that the injection length 𝑙 of the walkers should follow a specific statistical law the expression of which is given in Equations ( 10)-( 11) of [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF].This procedure gives quite accurate results, but the method of sampling injection lengths 𝑙, relying on the numerical inversion of the cumulative distribution function related to the random variable 𝑙, has a serious impact on the computation duration. A much simpler and faster injection method is now presented. The idea is inspired by the injection method consecutive to a volume power field which will be further detailed in Section 3. of this article. It should be reminded that the displacement of a Brownian walker representative of a 1D exclusively conductive transfer between two consecutive instants 𝑡 and

(𝑡 + 𝛿𝑡) is given by: 𝑥(𝑡 + 𝛿𝑡) -𝑥(𝑡) = √2𝑎𝛿𝑡 𝑟 𝑁 [START_REF] Kribus | The promise and challenge of solar volumetric absorbers[END_REF] where 𝑟 𝑁 is a dimensionless random number following the unit normal distribution. The idea of our new method is to inject 𝛿𝑀 𝑠 (𝑡) Brownian walkers at the boundary (as before), but the injection length 𝑙 mentioned above, the evaluation of which requires the numerical inversion of its cumulative distribution function, is not computed anymore. Instead, the walkers are injected at random instants (𝑡 + 𝑟 𝑈 𝛿𝑡) during the time step (𝑟 𝑈 being a random number following the uniform distribution over ]0 ; 1[), and they are submitted to an ordinary Brownian walk during the rest of the time step.

Consequently, their displacement, starting at the boundary, is obtained by:

𝑥(𝑡 + 𝛿𝑡) -𝑥(𝑡 + 𝑟 𝑈 𝛿𝑡) = √2𝑎(1 -𝑟 𝑈 )𝛿𝑡 𝑟 𝑁 ( 4
)
where 𝑥(𝑡 + 𝑟 𝑈 𝛿𝑡) is the position of the boundary, which is also the position of the injected walker at instant (𝑡 + 𝑟 𝑈 𝛿𝑡). We do not show a comparison of the two methods in this article because they give the same exact results. On the other hand, in our validation case described below, the new method was about forty times faster than the ancient one relying on use of the Newton-Raphson inversions of the 𝑙 cumulative distribution function.

In order to illustrate the validity of the new approach, we consider the problem of a 1D Cartesian where 𝜆 is the thermal conductivity of the wall. A dimensional analysis shows that for the particular problem treated here, the field 𝜃 * depends only on the dimensionless space 𝑥 * and time 𝑡 * variables.

Consequently, for the Brownian walker simulation of this problem, it can be shown that only four dimensionless quantities need to be provided: the number of voxels covering the thickness of the wall

𝑁 𝑥 = 𝐿 𝛿𝑥
, the number of time steps of the calculation 𝑁 𝑡 , the value of the parameter 𝛼 and the number of walkers 𝑀 𝑟 necessary to account for a temperature elevation of 𝜃 𝑟 within one voxel. Figure 1 shows the results of two Brownian walker simulations compared to the 𝜃 * field yielded at different instants 𝑡 * by a standard finite difference calculation. In the two graphs of Figure 1, the values retained for 𝑁 𝑥 and 𝑀 𝑟 are the same, namely 100 and 10 5 respectively. On the other hand, (𝛼, 𝑁 𝑡 ) = (10,100) on This figure should be printed in color.

2-column fitting image.

One can see on both graphs of Figure 1 that the two walker calculations coincide very well with the finite difference results, although the 𝛼 ratio is substantially higher than 2. This is an important All calculations presented in this paper were done with one standard 3.5 GHz processor. The RAM load is very low in all cases. All deterministic calculations were done with home-made Python algorithms using implicit schemes, and all stochastic calculations with home-made C++ algorithms, using only open source libraries. Parameters were kept the same for both methods. The deterministic calculation durations are not specified. In the 1D spatial configurations treated here, they were about a few seconds. So we must concede that our stochastic method is not competitive in terms of CPU time compared to other deterministic approaches, particularly when the asymptotic steady-state regime is wanted. Yet, our stochastic approach presents a strong advantage that makes it very interesting for our applications, namely its very low memory need. More classical deterministic approaches are strongly limited by this memory need, and thus do not allow the study of representative material structures, especially when treating transient conduction-radiation modeling involving complex 3D heterogeneous media.

Revisiting the treatment of a Dirichlet boundary condition with Brownian walkers

The strategy presented in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] (Section 2.5.1.) for applying a Dirichlet (imposed temperature) condition to a boundary voxel consisted in regulating the number of walkers within this voxel at the end of each time step. Doing so, it was shown that the parameter 𝛼 should be close to 2 in order the walker simulation to give correct results. But it has been shown in the preceding section that the constraint 𝛼 ≈ 2 is inexistent in the presence of Neumann or adiabatic boundary conditions. So the objective of the following studies is to revisit the way to account for Dirichlet boundary conditions with Brownian walkers in the hope to make the constraint 𝛼 ≈ 2 disappear.

Dirichlet condition, previous method: "regulation"

As explained in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], this method consists, at the end of each time step, after the displacement of the Brownian walkers, in simply regulating the numbers of Brownian walkers within the voxels submitted to the Dirichlet condition. Thus, if 𝑀 𝑖 Brownian walkers represent the imposed temperature elevation 𝜃 𝑖 , and if at the end of the time step a voxel 𝑛 submitted to the Dirichlet condition contains 𝑀 𝑛 Brownian walkers, then (𝑀 𝑖 -𝑀 𝑛 ) Brownian walkers are added to the voxel 𝑛.

In order to illustrate the results and limits of this approach, we consider the problem of a 1D Cartesian homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 0 and submitted at instants 𝑡 > 0 to a Dirichlet condition at 𝑥 = 0 and to an adiabaticity condition at 𝑥 = 𝐿. The dimensionless temperature elevation 𝜃 * is defined as

𝜃 * = 𝜃 𝜃 𝑟
where 𝜃 = 𝑇 -𝑇 0 is the temperature elevation and 𝜃 𝑟 is the temperature elevation imposed at 𝑥 = 0. Two Brownian walker calculations were performed with the following common data: number of voxels covering the thickness of the wall 𝑁 𝑥 = 100, number of walkers necessary to account for a temperature elevation of 𝜃 𝑟 within one voxel 𝑀 𝑟 = 10 5 . On the other hand, the values of the parameter 𝛼 and of the number of time steps 𝑁 𝑡 were (𝛼, 𝑁 𝑡 ) = (2,9000) for the first calculation and (𝛼, 𝑁 𝑡 ) = (6,1000) for the second calculation (note that the product 𝛼 2 𝑁 𝑡 is conserved between the two calculations, which ensures that both calculations cover the same temporal interval). The two graphs of Figure 2 compare the walker results to the finite difference calculation. This figure should be printed in color.

2-column fitting image.

In Figure 2(a) corresponding to 𝛼 = 2, the walker results are quite well superposed with the finite difference ones, which is not at all the case in Figure 2(b) where 𝛼 = 6. This observation has already been made in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] and has led to the necessity to choose the time step 𝛿𝑡 in such a way that the ratio 𝛼 is close to 2 with resulting computation durations taking long times as indicated in the legend of Figure 2.

Dirichlet condition, new method: "injection"

A new method of management of Dirichlet boundary conditions is now presented. The idea is inspired by the Neumann condition injection method described in Section 2.1.. Instead of an abrupt walker number adjustment as with the "regulation" method, Brownian walkers are injected at the frontier submitted to the Dirichlet condition at random instants during the time step, and then they move in a Brownian way via the standard procedure. Some of them will eventually land in the voxel submitted to the imposed temperature, the others will diffuse further. This procedure is repeated until the number of walkers within the voxel reaches the number 𝑀 𝑖 representing the imposed temperature elevation at the voxel. In the particular case where the temperature of the voxel is higher than the imposed temperature, negative Brownian walkers must be injected. The notion of negative walker will be introduced and validated in Section 3.2.1.. The strategy described here is strictly equivalent to injecting a (positive or negative) heat flux at the Dirichlet boundary until the temperature at this boundary attains the prescribed one.

In order to illustrate the validity of this approach, let us consider the same problem as the one described in the preceding section. In Figure 3, two walker calculation results are presented. For these two calculations, the parameters 𝑁 𝑥 and 𝑀 𝑟 are 100 and 10 5 respectively. On the other hand, (𝛼, 𝑁 𝑡 ) = (6,1000) in Figure 3 This figure should be printed in color.

2-column fitting image.

With the new method, the results obtained with 𝛼 = 6 (Figure 3(a)) are quite satisfactory, even at times close to 0 when the number of walkers to be injected is probably very high (let us recall here that imposing an instantaneous temperature jump at a boundary implies the application of an infinite heat flux at this boundary at the very first instants). The results remain satisfactory when 𝛼 = 24 (Figure 3(b)), but the first time steps are not so well reproduced. This is simply due to the fact that compared to the calculation of Figure 3(a), the time step has been multiplied by 16, hence weakening the performance of the walker calculation at low instants when temporal gradients are very high. This set apart, both walker calculations coincide quite satisfactorily with the finite difference results, at the difference of the walker result obtained with the "regulation" method and with 𝛼 = 6 (Figure 2(b)). This is an important advancement for at least two reasons. First because of the impact on the computation duration: the results of Figure 3(b) are obtained after 63 time steps whereas 9000 time steps were necessary in Figure 2(a) to cover the same time interval (𝛼 2 𝑁 𝑡 being conserved between the two simulations), thus the computation is 140 times faster. Second, the new methodology theoretically eliminates another difficulty highlighted in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] (Section 3.2.) and appearing when dealing with high diffusivity contrasts in heterogeneous media: a variation of the 𝛼 ratio in space or time should be much more tractable with the new procedure.

Robin condition (convective heat exchange)

In [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] (Section 2.5.3.), a procedure for simulating convective heat losses with Brownian walkers was proposed in the specific case where the low temperature of the problem was the initial temperature of the medium and was also the temperature of the fluid, therefore where this temperature could be associated with a number of Brownian walkers equal to zero. This was a simplifying assumption in the sense that there was no need to introduce any Brownian walker originating from the fluid; only those originating from the medium and susceptible to leave it because of the convective exchange had to be considered. In order to account for this convective loss, a stochastic criterion for the transmission of a walker at the medium/fluid interface was introduced and demonstrated [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. This criterion involves a dimensionless quantity 𝑃 𝑐𝑜𝑛𝑣 expressed as:

𝑃 𝑐𝑜𝑛𝑣 = ℎ√𝜋𝛿𝑡 𝑒 ( 5
)
where ℎ is the convective heat transfer coefficient and 𝑒 the thermal effusivity of the medium. The comparison of a random number 𝑟 𝑈 following the uniform distribution over ]0 ; 1[ with 𝑃 𝑐𝑜𝑛𝑣 determines the behavior of a Brownian walker encountering the medium/fluid interface:

 If 𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 , the walker is transmitted to the fluid, thus it is no more considered in the calculation,  Otherwise it is specularly reflected.

We now present how to deal with the case where the fluid temperature cannot be represented by zero Brownian walkers. This is the case for example when the fluid heats the medium. To this end, two substantially different methods will be described: the first one is non-stochastic and is based on the expression of the convective heat flux, the second one is entirely stochastic and is inspired from the "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " criterion mentioned above.

Robin condition, "flux" method

Here, the expression of the convective heat flux is used:

𝜑 𝑓𝑙 (𝑡) = ℎ(𝑇 𝑓𝑙 -𝑇(𝑡)) (6)
where 𝑇 𝑓𝑙 is the temperature of the fluid and 𝑇(𝑡) the temperature of the voxel submitted to the Robin boundary condition. Then the method described in Section 2.1. for imposing a prescribed heat flux is applied: more precisely, knowing the temperature profile at instant 𝑡, the heat flux at this instant is evaluated using expression (6), a number of walkers corresponding to the value of the heat flux are injected at the Robin boundary, and from the Brownian movement of all the walkers (those present before the injection plus the ones injected) the temperature profile at instant (𝑡 + 𝛿𝑡) is deduced. This procedure can therefore be qualified as "explicit" in the sense that the convective heat flux at instant 𝑡 is used to evaluate the temperature profile at instant (𝑡 + 𝛿𝑡).

In order to illustrate the validity of this approach, we consider the problem of a 1D Cartesian This figure should be printed in color.

2-column fitting image.

An overall good correspondence between Brownian walkers and finite differences is observed.

There is a visible mismatch between walker and finite difference results at the very first instants for 𝐵𝑖 = 2, but this mismatch disappears quickly in the following time steps. The explanation for this is quite similar to the one given in the section concerning Dirichlet boundary conditions: a high 𝐵𝑖 value means a fast evolution of the temperature of the concerned boundary at the first instants which can not correctly be accounted for when the value of 𝛼 is high. In order to circumvent this difficulty, a solution consists in decreasing the value of 𝛼, as illustrated in Figure 5. In these two calculations, 𝑁 𝑥 and 𝑀 𝑟 have the classical values of 100 and 10 5 respectively, and 𝐵𝑖 has the high value of 10. In This figure should be printed in color.

2-column fitting image.

Finally, let us note that such results can also be obtained with 𝑇 𝑓𝑙 < 𝑇 0 , i.e. when the fluid cools the medium. In this situation, negative Brownian walkers must be injected, as we will see in Section 3.2.1..

Robin condition, "cell" method

In this section, the stochastic criterion involving the quantity 𝑃 𝑐𝑜𝑛𝑣 expressed in ( 5) is revisited.

The fluid temperature being now represented by a non-zero number of walkers, new "cells" outside the medium representing the fluid are introduced. In 1D this cell is a single voxel, in 3D it is a surface of voxels. The voxel(s) representing the fluid contain(s) a number of Brownian walkers 𝑀 𝑓𝑙 representative of its (their) temperature 𝑇 𝑓𝑙 . These walkers are going to move at each time step according to the same rules as those within the medium. In the case of a 1D wall submitted to a convective heating at 𝑥 = 0, and using dimensionless quantities, an artificial voxel of size 𝛿𝑥 * covering the interval -𝛿𝑥 * ≤ 𝑥 * ≤ 0 is created. Then the following procedure is applied at each time step:

1. The Brownian walkers of the medium and of the fluid cell are moved using the standard procedure (Equation ( 3)).

2. If a walker encounters the boundary 𝑥 * = 0, then it is submitted to the "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " transmission test described above.

3. If a walker encounters the boundary 𝑥 * = -𝛿𝑥 * , then it is reflected.

4. Once the walkers have reached their final positions, the temperature of the fluid cell is adjusted using the injection method described in Section 2.2.2.: new positive or negative Brownian walkers are injected at random instants during the time step at the initial position 𝑥 * = -𝛿𝑥 * . During their displacements, they follow the rules described above at the boundaries 𝑥 * = 0 and 𝑥 * = -𝛿𝑥 * .

It is worth mentioning that, when simulating Brownian walker kinetics, the thermophysical properties of the medium are needed to calculate the displacements of the walkers and to determine the number of walkers representative of a voxel temperature. On the other hand, in a thermal problem involving a convective exchange characterized by a heat transfer coefficient ℎ, the thermophysical properties of the fluid do not appear explicitly in the problem (in fact they are hidden within the coefficient ℎ), which means that it must be possible to simulate a convective heat exchange with walkers without introducing any fluid thermophysical properties. Consequently, we decide to take the fluid properties equal to those of the medium. This way, the displacement of a Brownian walker within the fluid can be expressed with the same diffusivity as within the medium. This also allows to quantify the number of Brownian walkers to consider in the fluid cell 𝑀 𝑓𝑙 : the reference number 𝑀 𝑟 being defined as the number of walkers accounting for a temperature elevation of (𝑇 𝑓𝑙 -𝑇 0 ) within one voxel of the medium, the fluid voxel having the same size as a medium voxel, and the density and specific heat of the fluid being taken equal to those of the medium, it follows immediately that

𝑀 𝑓𝑙 = 𝑀 𝑟 .
In order to test the approach described above, we consider again the problem of a 1D Cartesian homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 0 and submitted at instants 𝑡 > 0 to a convective heat exchange at 𝑥 = 0 with a fluid of temperature 𝑇 𝑓𝑙 > 𝑇 0 and to an adiabaticity condition at 𝑥 = 𝐿. Three walker simulations are performed with the following common parameters: This figure should be printed in color.

𝑁 𝑥 =

2-column fitting image.

As can be seen in Figure 6, the transient phase is quite well reproduced for 𝐵𝑖 = 0.1, but not at all for 𝐵𝑖 = 2 and 𝐵𝑖 = 10, which seems to indicate that the "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " transmission test is adequate only for low values of 𝑃 𝑐𝑜𝑛𝑣 . Moreover, in the third simulation, 𝑃 𝑐𝑜𝑛𝑣 is substantially greater than 1, which means that the test "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " is always positive and that walkers encountering the interface are systematically transmitted. This raises the following point: two thermal problems differing only by the values of their Biot numbers and such that the 𝑃 𝑐𝑜𝑛𝑣 parameters associated with these two Biot numbers are both greater than 1 will yield exactly the same results when treated by Brownian walkers using the "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " test, which of course is not conceivable. Consequently, the test "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " is not adapted for high values of 𝐵𝑖.

In order to improve this test, let us imagine a solid (medium 1), at uniform temperature 𝑇 1 , occupying the half-space 𝑥 ≤ 0, in contact with a fluid (medium 2), at uniform temperature 𝑇 2 , occupying the half-space 𝑥 ≥ 0. 𝑇 2 is assumed to be lower than 𝑇 1 . It was shown in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] that the surface flux crossing the interface 𝑥 = 0 in the direction {medium 1 → medium 2} during the time step 𝛿𝑡 due to the movements of the walkers of medium 1 is

𝜑 1→2 = 𝑒 1 (𝑇 1 -𝑇 𝑟 ) √𝜋𝛿𝑡
where 𝑇 𝑟 is the reference temperature for enthalpy. Identically, the surface flux crossing the interface 𝑥 = 0 in the direction {medium 2 → medium 1} during the time step 𝛿𝑡 due to the movements of the walkers of medium 2

is

𝜑 2→1 = 𝑒 2 (𝑇 2 -𝑇 𝑟 ) √𝜋𝛿𝑡 .
As explained above, in our approach the thermophysical properties of the two media are considered equal, so that the indices "1" and "2" may be dropped in the effusivities.

Consequently, the net surface flux in the direction {medium 1 → medium 2} due to the kinetics of the walkers only is

𝜑 𝑛𝑒𝑡 = 𝜑 1→2 -𝜑 2→1 = 𝑒(𝑇 1 -𝑇 2 )

√𝜋𝛿𝑡

, and from this expression a thermal resistance 𝑅 𝑛𝑒𝑡 appears, expressed as

𝑅 𝑛𝑒𝑡 = 𝑇 1 -𝑇 2 𝜑 𝑛𝑒𝑡 = √𝜋𝛿𝑡 𝑒
. Now, if the fluid (medium 2) is moving, the convective heat transfer at the interface is characterized by a coefficient ℎ, and the convective surface

flux in the direction {medium 1 → medium 2} is 𝜑 𝑐𝑜𝑛𝑣 = ℎ(𝑇 1 -𝑇 2 ) = 𝑇 1 -𝑇 2 𝑅 𝑐𝑜𝑛𝑣
where 𝑅 𝑐𝑜𝑛𝑣 is the classical thermal resistance related to convection and is expressed as

𝑅 𝑐𝑜𝑛𝑣 = 1 ℎ
. In order to merge the effects of the kinetics of the walkers and of the convective heat transfer, an equivalent electrical circuit composed of the two resistances 𝑅 𝑛𝑒𝑡 and 𝑅 𝑐𝑜𝑛𝑣 connected in parallel between the temperatures 𝑇 1 and 𝑇 2 is proposed. This topology makes sense in the following two extremes situations: (i) if ℎ = 0, 𝑅 𝑐𝑜𝑛𝑣 → ∞ and the heat transfer between media 1 and 2 is only through the resistance 𝑅 𝑛𝑒𝑡 , which is consistent with the fact that there is no convective heat transfer; (ii) if ℎ → ∞, 𝑅 𝑐𝑜𝑛𝑣 ≪ 𝑅 𝑛𝑒𝑡 and the heat transfer between media 1 and 2 is essentially through the resistance 𝑅 𝑐𝑜𝑛𝑣 , meaning that convection dominates the other heat transfer phenomena at the interface. With this parallel model, the total surface flux from medium 1 to medium 2 is

𝜑 𝑡𝑜𝑡 = 𝑇 1 -𝑇 2 𝑅 𝑛𝑒𝑡 ∥ 𝑅 𝑐𝑜𝑛𝑣 = (ℎ + 𝑒 √𝜋𝛿𝑡 ) (𝑇 1 -𝑇 2 )
and the convective surface flux is 𝜑 𝑐𝑜𝑛𝑣 = ℎ(𝑇 1 -𝑇 2 ), so that the ratio However, the walker results still do not exactly coincide with those from finite differences. Progress has been made, but the method still needs improvement. More generally, one of the main challenges in the future is a better management of interfaces. In Section 3.2. of our previous article [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], a solution based on a transmission probability was presented for the determination of the behavior of Brownian walkers at interfaces. Yet, Section 3.3. of this same paper highlighted the weaknesses of this approach when considering a high contrast between the thermophysical properties of the two constituents on either side of the interface. Our future research aims to solve this difficulty, and we believe that progress in this direction could improve the results of Figure 6 using our "cell" method.

Finally, it is also worth mentioning that, unlike the "flux" method presented in Section 2.3.1., the "cell" method introduced here does not raise any particular difficulty at the first time steps, which allows us to obtain better results more rapidly at the first time steps for high 𝐵𝑖 values. This figure should be printed in color.

2-column fitting image.

Conduction-radiation coupling in an optically thick and grey semi-transparent medium

In all of the studies carried out above, conduction was considered as being the only heat transfer mode present within the medium. However, the media targeted in our studies are semi-transparent and devoted to applications at high temperatures where the transfer by radiation cannot be neglected and can even become the dominant transfer mode depending on the microstructure and the radiative properties of the constituents. This section presents our first modeling studies, by Brownian walkers, of a non-linear conduction-radiation coupling in the transient state within a voxelized structure, in an optically thick and grey medium.

Conduction-radiation coupling using the radiative conductivity

In this study, it is supposed that the semi-transparent sample is optically thick. In these conditions, the radiative flux within the medium can be represented (except in the vicinity of its walls) by the Rosseland law (also called the radiative Fourier law):

𝜑 ⃗ 𝑅 = -𝜆 𝑅 (𝑇) ∇ ⃗ ⃗ 𝑇 (9)
where 𝜑 ⃗ 𝑅 stands for the radiative flux and 𝜆 𝑅 (𝑇) is the radiative conductivity, which depends strongly on the temperature 𝑇. If it is supposed in addition that the material constituting the sample is grey, then 𝜆 𝑅 (𝑇) has an expression of the type 𝜆 𝑅 (𝑇) = 𝑘𝑇 3 , 𝑘 being a constant of expression 𝑘 =

16𝑛 𝑜 2 𝜎 𝑆𝐵 3𝛽(1-𝑔𝜔)
where 𝑛 𝑜 is the refractive index of the medium, 𝛽 its extinction coefficient, 𝜔 its scattering albedo, 𝑔 the asymmetry factor of its scattering phase function and 𝜎 𝑆𝐵 the Stefan-Boltzmann constant [START_REF] Howell | Thermal Radiation Heat Transfer, Seventh Edition[END_REF].

By applying this model, the total conductivity 𝜆 𝑇 (𝑇), defined by 𝜆 𝑇 (𝑇) = 𝜆 + 𝜆 𝑅 (𝑇) and giving account of the conduction-radiation coupling, depends on the temperature, and the same goes for the total diffusivity 𝑎 𝑇 (𝑇) and the total effusivity 𝑒 𝑇 (𝑇). Thus the calculation of the displacement of a walker between instants 𝑡 and (𝑡 + 𝛿𝑡) will depend on the temperature of its departure voxel at instant 𝑡 but also on the temperatures of all the other voxels that it is going to cross. The situation treated here is actually analogous to the one studied in Section 3.1. of [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]: it must be considered that each voxel is made of a different constituent of properties 𝑎 𝑇 (𝑇(𝑡)) and 𝑒 𝑇 (𝑇(𝑡)) where 𝑇(𝑡) is the temperature of the voxel at instant 𝑡. Consequently, any change in voxel during the displacement of a walker must be regarded as a problem of interface between two different constituents with the associated "𝑟 𝑈 < 𝑒 2 𝑒 1

" transmission test where 𝑟 𝑈 is a random number following the uniform distribution over ]0 ; 1[ and 𝑒 1 and 𝑒 2 are the effusivities of the voxels before and after transmission respectively [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF]. Moreover, the diffusivities 𝑎 1 and 𝑎 2 of these two voxels being a priori different, the same goes for the 𝛼 values of the two constituents, and this discrepancy was a source of difficulty in [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] because of the necessity to have all 𝛼 values close to 2. Thanks to our recent progress concerning the treatment of Dirichlet boundary conditions presented in Section 2.2., the constraint 𝛼 ≈ 2 is now obsolete.

In order to validate this approach, we consider the problem of a 1D Cartesian homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 𝐿 and submitted at instants 𝑡 > 0 to an imposed temperature 𝑇 0 > 𝑇 𝐿 at 𝑥 = 0, the temperature 𝑇 𝐿 being maintained at 𝑥 = 𝐿. The intensity of the conduction-radiation coupling is chosen such that 𝜆 𝑇 (𝑇 0 ) = 5 𝜆 𝑇 (𝑇 𝐿 ), condition that makes it possible to evaluate the multiplicative constant 𝑘 in the expression of the radiative conductivity. In this situation, the ratio each elementary volume of the semi-transparent medium emits and absorbs radiation, which results in an overall (positive or negative) radiative power injected within the volume. This volume radiative power, resulting from the radiative transfer, deforms the asymptotic temperature profile in a convex way (let us recall that without radiation, this asymptotic temperature profile would be a straight line).

In Figure 8(a), the time step 𝛿𝑡 is constant and equal to 𝑡 𝑓 𝑁 𝑡

. A non-negligible delay is observed between the Brownian walker temperature profiles and their counterparts obtained by finite differences during the transient phase. Note that this shift disappears when the asymptotic steadystate regime is reached (black points and curve). In Figure 8(b), the same simulation was done but with a time step evolving with the advancement of the simulation. More precisely, the time step was chosen to follow a geometric evolution over the calculation interval [0 ; 𝑡 𝑓 ], i.e. of the form 𝛿𝑡 𝑛 𝑡 = 𝛿𝑡 0 × 𝑞 𝑛 𝑡 , the values of 𝑡 𝑓 , 𝑁 𝑡 and 𝑞 allowing the evaluation of 𝛿𝑡 0 (note in particular that 𝑞 = 1 corresponds to a constant time step). As seen in Figure 8(b), with this geometric evolution and with the parameter 𝑞 equal to 1.05, the delay between the walker and finite difference results is much more acceptable than in Figure 8(a), and the computation duration is shorter. Our interpretation of the slightly imperfect superposition of the results is that a high temperature gradient appears at the boundary 𝑥 = 0 at the first instants, which results in a high total effusivity contrast between neighboring voxels at this boundary. The "𝑟 𝑈 < 𝑒 2 𝑒 1

" transmission criterion between two neighboring voxels involving the ratio of their total effusivities, this criterion is always fulfilled in one direction (the one from the medium to the boundarybecause 𝑒 𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 ≫ 𝑒 𝑚𝑒𝑑𝑖𝑢𝑚 ), but almost never in the opposite. Let us also recall that an instantaneous temperature jump prescribed at a boundary requires to impose an infinite flux at this boundary at the first instants, so it is not so surprising to encounter difficulties at the first time steps.

The very good overall agreement between the walker and finite difference results validates the procedure proposed in our previous paper [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] for managing the reflection or the transmission of a walker at the interface between two constituents within a heterogeneous medium. In addition, this procedure has permitted to solve a non-linear transient thermal problem. That said, it has several weaknesses. First, it is limited to the very particular situation of an optically thick and grey medium for which the Rosseland radiative conductivity is expressed as a simple function of the temperature.

In addition, it requires treating any change of voxel as an interface encounter problem and to apply the associated procedure. This aspect generates a substantial increase in the calculation time.

The conduction-radiation coupling within semi-transparent media is commonly treated by the introduction of a radiative power density source term 𝑝 𝑅 (𝑥, 𝑡), expressed in W.m -3 , in the energy balance equation [START_REF] Howell | Thermal Radiation Heat Transfer, Seventh Edition[END_REF]. In order to extend the capacities of the stochastic approach by Brownian walkers to the resolution of more general conduction-radiation coupling problems, it is preferable to establish a procedure that models a supply of enthalpy by a radiative power density source within the medium.

Conduction-radiation coupling using the radiative power density field

As explained above, the conduction-radiation coupling is commonly treated by the introduction of a radiative power density source term 𝑝 𝑅 (𝑥, 𝑡). To achieve such a model with Brownian walkers, a method allowing to take account of a power density field in an energy balance equation must first be defined and validated. Then the problem of Section 3.1. will be treated again, this time using the radiative power density field.

Treatment of a thermal power density field with Brownian walkers -Negative walkers

A positive thermal power density provides enthalpy to the medium, so it can be simulated, using Brownian walkers, simply by the ad hoc creation of walkers in the concerned voxels. On the contrary, a negative thermal power density extracts enthalpy from the medium and results in the disappearance of walkers in the concerned zones. A supply of thermal power density of intensity 𝑝 in a voxel of volume 𝑉 𝑣 deposits an enthalpy 𝐻 𝑝 equal to 𝑝𝑉 𝑣 𝛿𝑡 during the time step 𝛿𝑡. This added enthalpy results in the injection of 𝑀 𝑝 = 𝐻 𝑝 𝛿ℎ 𝑟 new walkers at each time step. In one voxel, each walker are injected at a random position and a random instant during the time step, consequently their displacement is given by equation ( 4). This procedure has been proposed and validated in [START_REF] Gonneau | Modélisation du transfert thermique par marcheurs Browniens dans des milieux hétérogènes[END_REF].

Situations with a high thermal emission within a semi-transparent medium can result in locally negative values of the radiative power density, hence in the suppression of a substantial quantity of walkers, so that the number of walkers can become locally negative. This particular configuration however does not result in any physical inconsistency. Indeed, a voxel that does not contain any walker corresponds to the reference temperature 𝑇 𝑟 , not the absolute zero. However, the walker counting procedure must be adapted to this situation. Consider the vector 𝑀(𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗ = {𝑀 𝑛 (𝑡)} where 𝑀 𝑛 (𝑡) represents the number of walkers in voxel 𝑛 at instant 𝑡. To obtain its counterpart at instant (𝑡 + 𝛿𝑡), the following procedure is applied:

1. 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is initiated at 0 ⃗ .

For each voxel 𝑛:

a. If 𝑀 𝑛 (𝑡) > 0, 𝑀 𝑛 (𝑡) walkers are displaced, and for each of these walkers 1 is added to the population of its arrival voxel in the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ .

b. Otherwise, (-𝑀 𝑛 (𝑡)) walkers are displaced, and for each of these walkers 1 is subtracted from the population of its arrival voxel in the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ . These walkers behave as "negative walkers" in the sense that, when managed that way, they behave as transporting the enthalpy (-𝛿ℎ 𝑟 ).

In order to illustrate the validity of this approach, we consider the problem of a 1D Cartesian homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 0 and submitted at instants 𝑡 > 0 to an imposed time-independent thermal power density profile 𝑝(𝑥) having the following expression:

𝑝(𝑥) = 𝑝 0 𝑠𝑖𝑛 ( 2𝜋𝑥 𝐿 ) (10)
where 𝑝 0 is a constant. Moreover, the two boundaries 𝑥 = 0 and 𝑥 = 𝐿 are assumed to be adiabatic.

The expression retained for 𝑝(𝑥) is positive for 0 ≤ 𝑥 ≤ and 𝛼 = 10. The results of the walker simulations are compared to the ones obtained by finite differences at different instants in Figure 9. The very good overall agreement between the two approaches validates the procedure proposed above. We now have the tools to model a conductionradiation coupling using a radiative power density field with Brownian walkers. Single column fitting image.

Conduction-radiation coupling using the radiative power density field and Brownian walkers

In this section, as in Section 3.1., the semi-transparent sample is assumed to be grey and optically thick. The radiative flux is then expressed by the radiative Fourier law 𝜑 ⃗ 𝑅 = -𝜆 𝑅 (𝑇) ∇ ⃗ ⃗ 𝑇, and the radiative conductivity as 𝜆 𝑅 (𝑇) = 𝑘𝑇 3 . Yet, unlike in Section 3.1., we are not going to use the notion of total conductivity 𝜆 𝑇 (𝑇); instead, we are going to write the energy balance equation including the radiative power density source term 𝑝 𝑅 , this term being equal to (-∇ ⃗ ⃗ • 𝜑 ⃗ 𝑅 ). In a 1D

Cartesian problem, this equation writes:

𝜌𝐶 𝑝 𝜕𝑇 𝜕𝑡 = 𝜆 𝜕 2 𝑇 𝜕𝑥 2 + 𝑝 𝑅 (11)
with:

𝑝 𝑅 = - 𝜕𝜑 𝑥 𝑅 𝜕𝑥 = - 𝜕 𝜕𝑥 (-𝑘𝑇 3 𝜕𝑇 𝜕𝑥 ) = 𝑘 4 𝜕 2 (𝑇 4 ) 𝜕𝑥 2 (12)
Consequently, the problem can in principle be solved using the method presented in Section 3.2.1. by transforming the values of the radiative power density field into numbers of (positive or negative)

Brownian walkers to inject in the voxels. However, this approach "𝑝 𝑅 + walkers" brings specific difficulties: using the Brownian walker method produces noisy temperature profiles, making the numerical evaluations of the double derivatives of 𝑇 4 very hazardous. This statistical noise can be attenuated with a high quantity of Brownian walkers, at the expense of an increased computation duration. An alternative solution is to apply polynomial adjustments to the noisy temperature profiles and then to make the derivations on the polynomials. Using 8 degree polynomial adjustments was generally adapted.

In order to illustrate the validity of this approach, we consider the same problem as in Section 3.1.. A second difficulty specific to this particular configuration and concerning the calculation of the radiative power density field 𝑝 𝑅 appeared. When dealing with a thermal power density field, the standard procedure consists in evaluating the temperature profile at instant (𝑡 + 𝛿𝑡) based on the power density at instant 𝑡 (explicit scheme). But here, because of the instantaneous temperature jump imposed at the boundary 𝑥 = 0, the 𝑝 𝑅 term takes enormous values close to this boundary at instant 𝑡 = 0, which results in a divergence of the calculation. In order to solve this problem, the following method was applied. The vector 𝑀(𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗ = {𝑀 𝑛 (𝑡)} is again introduced, where 𝑀 𝑛 (𝑡) represents the number of walkers in voxel 𝑛 at instant 𝑡. To deduce 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ from 𝑀(𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗ , the following procedure is applied:

1. 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is initiated at 0 ⃗ .

2. For each voxel 𝑛, |𝑀 𝑛 (𝑡)| walkers are displaced, and for each of these walkers 1 or (-1) is added to the population of its arrival voxel in the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ according to the sign of 𝑀 𝑛 (𝑡).

3. The procedure described in Section 2.2.2. is used to impose the Dirichlet boundary conditions to the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ . This step is crucial to compute the correct radiative power density field because it avoids the divergence of the calculation.

4. The radiative power density field is evaluated based on the temperature profile contained in the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ .

5. The (positive or negative) walkers created by the radiative power density field are injected and displaced, and the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is updated accordingly.

6. Finally, the procedure described in Section 2.2.2. is used once again to make sure that the Dirichlet boundary conditions are imposed to the vector 𝑀(𝑡 + 𝛿𝑡) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ at the end of the time step.

Figure 10 shows the results of two walker simulations done with the parameters listed in Section 3.1.. The value retained for 𝑀 𝑟 (10 5 ) proved to be sufficiently high, so that the double derivatives of the 𝑇 4 profiles could be evaluated without resorting to polynomial adjustments.

The results of Figure 10(a) were obtained with a constant time step. One notes a discrepancy between the walker and finite difference profiles at the first instants of the transient state, particularly at the very first time step. Again, this discrepancy is attributed to the difficulty to calculate the second derivative of the 𝑇 4 profile at the very first instants when the temperature gradient is very high. This figure should be printed in color.

2-column fitting image.

Conclusion and perspectives

In this article, we have taken interest in the improvement of methods initially presented in a previous contribution [START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF] regarding the simulation by Brownian walkers of transient thermal conduction within heterogeneous media represented by voxelized structures. In particular, a new procedure has been described for the management of Dirichlet boundary conditions, allowing to get rid of a constraint on the choice of the time step initially imposed by the criterion 𝛼 = √2𝑎𝛿𝑡 𝛿𝑥 ≈ 2. This important progress makes it possible to perform much faster simulations and to solve difficulties related to strong effusivity contrasts in heterogeneous media. A new procedure, much simpler and faster, for the management of Neumann boundary conditions has also been proposed. We have also extended the field of application of our stochastic approach for the management of Robin boundary conditions, both in terms of the temperature of the fluid and of the intensity of the heat exchange. In the last section of this paper, a method allowing to account for a thermal power density field using Brownian walkers has been introduced, and most importantly, our first simulations of a non-linear transient conduction-radiation coupling in an optically thick and grey semi-transparent medium using two different methods based only on the use of Brownian walkers have been successfully achieved.

In future research, progress has to be made in the modeling of convective heat exchanges with Brownian walkers for high Biot numbers. But the long-term objectives of our work go beyond a stochastic conductive model and aim at coupling stochastic approaches simulating conduction and radiation in semi-transparent 3D voxelized structures. Now that our conduction model by Brownian walkers is quite complete, we are going to examine its performance over 3D voxelised representations of heterogeneous media, and we are also going to investigate the ways to exchange energy between photon bundles for radiation and Brownian walkers for conduction for a precise and entirely stochastic conduction-radiation coupling.

  homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 0 and submitted at instants 𝑡 > 0 to an adiabaticity condition at 𝑥 = 𝐿 and to a constant heat flux 𝜑 at 𝑥 = 0. The following dimensionless parameters 𝑥 * = 𝑥 𝐿 and 𝑡 * = 𝑎𝑡 𝐿 2 are introduced. The ratio 𝛼 defined in (1) can be regarded as a dimensionless representation of the time step 𝛿𝑡 for a given space step 𝛿𝑥. Lets us introduce the dimensionless temperature elevation field 𝜃 * defined by 𝜃 * = 𝑇-𝑇 0 𝜃 𝑟 where 𝑇 is the temperature field within the wall and 𝜃 𝑟 is a reference temperature elevation defined by 𝜃 𝑟 = 𝜑𝐿 𝜆

Figure 1 (

 1 Figure 1(a), whereas (𝛼, 𝑁 𝑡 ) = (20,25) on Figure 1(b). Note that multiplying 𝛼 by 2 (10 → 20) means multiplying 𝛿𝑡 by 4, so that the time interval covered by the calculation is conserved when 𝑁 𝑡 is divided by 4 (100 → 25).

Figure 1 :

 1 Figure 1: Study of the influence of the value of the 𝛼 ratio on the results of a Brownian walker simulation with Neumann and adiabatic boundary conditions. Figure 1(a) (left): 𝛼 = 10, computation duration = 59 s. Figure 1(b) (right): 𝛼 = 20, computation duration = 15 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  result because the possibility to multiply 𝛼 by 10 in the calculations (2 → 20) means the possibility to multiply the time step by 100 and consequently, in most cases, to divide the computation duration by the same amount. The two walker calculations presented in Figure 1 took 59 s (calculation (a) -𝛼 = 10) and 15 s (calculation (b) -𝛼 = 20). The division of the computation duration by a factor of 4 is perfectly consistent with the doubling of the value of 𝛼.

Figure 2 :

 2 Figure 2: Study of the influence of the value of the 𝛼 ratio on the results of a Brownian walker simulation with Dirichlet and adiabatic boundary conditions. The Dirichlet boundary condition is treated by the "regulation" method. Figure 2(a) (left): 𝛼 = 2, computation duration = 19220 s.

Figure 2 (

 2 Figure 2(b) (right): 𝛼 = 6, computation duration = 2020 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  (a) and (𝛼, 𝑁 𝑡 ) = (24,63) in Figure 3(b).

Figure 3 :

 3 Figure 3: Study of the influence of the value of the 𝛼 ratio on the results of a Brownian walker simulation with Dirichlet and adiabatic boundary conditions. The Dirichlet boundary condition is treated by the "injection" method. Figure 3(a) (left): 𝛼 = 6, computation duration = 2173 s.

Figure 3 (

 3 Figure 3(b) (right): 𝛼 = 24, computation duration = 134 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  homogeneous wall of thickness 𝐿, initially isothermal at temperature 𝑇 0 and submitted at instants 𝑡 > 0 to a convective heat exchange at 𝑥 = 0 with a fluid of temperature 𝑇 𝑓𝑙 > 𝑇 0 and to an adiabaticity condition at 𝑥 = 𝐿. The dimensionless temperature elevation 𝜃 * is now defined as 𝜃 * = 𝜃 𝜃 𝑟 where 𝜃 = 𝑇 -𝑇 0 is the temperature elevation and 𝜃 𝑟 is equal to (𝑇 𝑓𝑙 -𝑇 0 ). A dimensional analysis shows that in this problem, a new dimensionless parameter appears: the Biot number 𝐵𝑖 = ℎ𝐿 𝜆 where ℎ is the heat transfer coefficient and 𝜆 the thermal conductivity of the medium. Figure 4 shows the results of two walker simulations. The parameters common to these two calculations are 𝑁 𝑥 = 100, 𝑀 𝑟 = 10 5 and 𝛼 = 20. On the other hand, (𝑁 𝑡 , 𝐵𝑖) = (2000,0.1) in Figure 4(a) and (𝑁 𝑡 , 𝐵𝑖) = (200,2) in Figure 4(b). The fact that 𝑁 𝑡 is taken much larger in the first calculation than in the second is justified by the values retained for the Biot number: in the first calculation, 𝐵𝑖 is low, hence the duration to reach the asymptotic steady state is long. It is of course the contrary in the case of the second calculation.

Figure 4 :

 4 Figure 4: Study of the influence of the value of the Biot number on the results of a Brownian walker simulation with Robin and adiabatic boundary conditions. The Robin boundary condition is treated by the "flux" method. Figure 4(a) (left): 𝐵𝑖 = 0.1, computation duration = 4069 s. Figure 4(b) (right): 𝐵𝑖 = 2.0, computation duration = 438 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

Figure 5 (

 5 Figure 5(a), 𝛼 and 𝑁 𝑡 are 15 and 200 respectively, and in Figure 5(b), these two parameters are 5 and1800 respectively (note that the product 𝛼 2 𝑁 𝑡 is conserved between the two calculations, which ensures that both calculations cover the same temporal interval). Decreasing the value of 𝛼 is clearly efficient, but at the expense of an increased computation duration.

Figure 5 :

 5 Figure 5: Study of the influence of the value of the parameter 𝛼 on the results of a Brownian walker simulation with Robin and adiabatic boundary conditions. The Robin boundary condition is treated by the "flux" method. Figure 5(a) (left): 𝛼 = 15, computation duration = 443 s. Figure 5(b) (right): 𝛼 = 5, computation duration = 3879 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  100, 𝑀 𝑟 = 10 5 and 𝛼 = 20. The values of 𝑁 𝑡 and 𝐵𝑖 are (𝑁 𝑡 , 𝐵𝑖) = (2000,0.1) for the first simulation, (200,2) for the second simulation, and (100,10) for the third simulation. From these data, the parameter 𝑃 𝑐𝑜𝑛𝑣 , which has for expression 𝑃 𝑐𝑜𝑛𝑣 = √𝜋 2 ⁄ 𝛼 𝐵𝑖 𝑁 𝑥 ⁄ , can be quantified: 𝑃 𝑐𝑜𝑛𝑣 = 0.025 for the first simulation, 𝑃 𝑐𝑜𝑛𝑣 = 0.496 for the second simulation, and 𝑃 𝑐𝑜𝑛𝑣 = 2.482for the third simulation.

Figure 6 :

 6 Figure 6: Study of the influence of the value of the Biot number (and its associated 𝑃 𝑐𝑜𝑛𝑣 value) on the results of a Brownian walker simulation with Robin and adiabatic boundary conditions. The Robin boundary condition is treated by the "cell" method. Figure 6(a) (top left): 𝐵𝑖 = 0.1, computation duration = 2147 s. Figure 6(b) (top right): 𝐵𝑖 = 2.0, computation duration = 254 s.

Figure 6 (

 6 Figure 6(c) (bottom): 𝐵𝑖 = 10, computation duration = 128 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

Figure 7

 7 Figure7shows the results of the Brownian walker calculations using the new "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 ′ " transmission test in the same configurations as the ones examined in Figure6. The numerical values of 𝑃 𝑐𝑜𝑛𝑣 ′

Figure 7 :

 7 Figure 7: Study of the impact of the new "𝑟 < 𝑃 𝑐𝑜𝑛𝑣 ′ " transmission test on the results of Brownian walker calculations with Robin and adiabatic boundary conditions and for different values of the Biot number. Figure 7(a) (top left): 𝐵𝑖 = 0.1, computation duration = 2117 s. Figure 7(b) (top right): 𝐵𝑖 = 2.0, computation duration = 239 s. Figure 7(c) (bottom): 𝐵𝑖 = 10, computation duration = 118 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  07 for 𝑇 𝐿 and 4.34 for 𝑇 0 . In other words, with the data retained, the heat transfer is essentially conductive in the vicinity of the cold boundary and dominated by radiation in the vicinity of the hot boundary. The numerical values of the input data are: 𝐿 = 0.1 m, 𝑁 𝑥 = 100, 𝑡 𝑓 = 2000 s (𝑡 𝑓 is the final instant of the simulation), 𝑁 𝑡 = 100, 𝑀 𝑟 = 10 5 , 𝑇 0 = 2000K, 𝑇 𝐿 = 500K, 𝜆 = 1 W.m -1 .K -1 , and 𝜌𝐶 𝑝 = 10 6 J.m -3 .K -1 .

Figure 8 :

 8 Figure 8: Temperature profiles in an optically thick and grey medium calculated by Brownian walkers using the radiative conductivity approach. Figure 8(a) (left): with a constant time step; 𝛼 = 6.6 to 14.8, computation duration = 818 s. Figure 8(b) (right): with a progressive time step (tn = n-th instant of the temporal discretization), 𝛼 = 1.2 to 32.7, computation duration = 530 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

Figure 8

 8 Figure8shows the results of two Brownian walker simulations of the problem. In both calculations, a change of curvature (from concave to convex) of the 𝑇(𝑥 * ) profile with time is observed. The convex asymptotic steady-state profile is a classical result due to the radiative transfer:

2 ≤𝑇

 2 𝑥 ≤ 𝐿, so negative walkers must be introduced at each time step in the right half of the wall. A dimensional analysis concludes that the dimensionless temperature elevation field 𝜃 * defined by 𝜃 * = of the dimensionless space and time variables 𝑥 * = 𝑥 𝐿 and 𝑡 * = 𝑎𝑡 𝐿 2 only. The following numerical data were retained for the Brownian walker simulations: 𝑁 𝑥 = 100, 𝑁 𝑡 = 30, 𝑀 𝑟 = 10 5

Figure 9 :

 9 Figure 9: Temperature profiles in a wall submitted to a thermal power density field taking positive and negative values. Computation duration = 2 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

Figure 10 (

 10 b) shows the results of the same simulation but with a time step evolving with the advancement of the simulation (of the form 𝛿𝑡 𝑛 𝑡 = 𝛿𝑡 0 × 𝑞 𝑛 𝑡 as explained in Section 3.1.). The walker results are now significantly better, even at the very first time step (one point set apart). Moreover, the computation duration is significantly faster. Compared to the method based on the total conductivity presented in Section 3.1., this approach appears to be more precise in the transient regime and quite faster, with a gain of around 6 depending on the choice of the time steps. Thus, a second procedure for solving a non-linear transient conduction-radiation problem in an optically thick and grey semitransparent medium with Brownian walkers is validated.

Figure 10 :

 10 Figure 10: Temperature profiles in an optically thick and grey medium calculated by Brownian walkers using the radiative power density approach. Figure 10(a) (left): with a constant time step; 𝛼 = 6.35, computation duration = 131 s. Figure 10(b) (right): with a progressive time step (tn = nth instant of the temporal discretization), 𝛼 = 1.2 to 14, computation duration = 78 s. Discrete markers = Brownian walkers; continuous curves = finite differences.

  Following a reasoning identical to the one developed in[START_REF] Gonneau | Modeling of heat transfer within heterogeneous media by Brownian walkers[END_REF], one comes to the conclusion that the former "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " transmission criterion should be replaced by a new "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣

		𝜑 𝑐𝑜𝑛𝑣 𝜑 𝑡𝑜𝑡	=	ℎ +	ℎ √𝜋𝛿𝑡 𝑒	=	𝑃 𝑐𝑜𝑛𝑣 𝑃 𝑐𝑜𝑛𝑣 + 1	with 𝑃 𝑐𝑜𝑛𝑣 =	ℎ√𝜋𝛿𝑡 𝑒	(7)
										′	" criterion with
	𝑃 𝑐𝑜𝑛𝑣 ′	expressed as:
						𝑃 𝑐𝑜𝑛𝑣 ′	=	𝑃 𝑐𝑜𝑛𝑣 𝑃 𝑐𝑜𝑛𝑣 + 1	(8)
		Let us note several interesting features of the new quantity 𝑃 𝑐𝑜𝑛𝑣 ′	which has emerged from our
	analysis: (i) when 𝑃 𝑐𝑜𝑛𝑣 ≪ 1, 𝑃 𝑐𝑜𝑛𝑣 ′	≈ 𝑃 𝑐𝑜𝑛𝑣 , which means that when 𝐵𝑖 is small, Brownian walker
	calculations with either of the two tests "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 " and "𝑟 𝑈 < 𝑃 𝑐𝑜𝑛𝑣 ′	" will give approximately the
	same results; (ii) unlike 𝑃 𝑐𝑜𝑛𝑣 which can take any positive value, 𝑃 𝑐𝑜𝑛𝑣 ′	is always lower than 1, which
	ensures that whatever the value of 𝐵𝑖, a fraction of the walkers encountering the interface will be
	reflected; (iii) 𝑃 𝑐𝑜𝑛𝑣
										𝜑 𝑐𝑜𝑛𝑣 𝜑 𝑡𝑜𝑡	is
	equal to:

′ is a strictly increasing function of 𝑃 𝑐𝑜𝑛𝑣 and hence a strictly increasing function of 𝐵𝑖, so that two different Biot numbers will produce two different 𝑃 𝑐𝑜𝑛𝑣 ′ quantities, both lower than 1, which in turn will produce two different results using Brownian walkers; (iv) finally, when 𝐵𝑖 becomes very large, 𝑃 𝑐𝑜𝑛𝑣 ′ → 1, which means that all walkers encountering the interface are transmitted: this is perfectly consistent with our approach described above for simulating a Dirichlet boundary condition with walkers.

Acknowledgements

This work was financially supported by CEA (the French national atomic energy commission) and CNRS (the French national scientific research center).