Raj Narayan Konduru

Olivier Farges

Gilles Parent

Vincent Schick

SEMI MESH-LESS MODEL FOR COMBINED CONDUCTION-CONVECTION-RADIATIVE HEAT TRANSFER : APPLICATION TO POROUS HEAT EXCHANGER

Keywords: Radiation, Convection, Heat exchanger, Porous media, High temperature, Monte Carlo Method

Existing deterministic methods can solve coupled heat transfer problems in complex 3D geometries, but the increasing complexity results in high computational time and memory usage due to grid-based approximations. This paper presents a semi-meshless method for solving coupled heat transfer problems. This method is based on a stochastic approximation, i.e., the Monte Carlo method. The resulting algorithm is less sensitive to the system's complexity. Thanks to advancements in computer acceleration technology, it is now easier to implement this method. This coupled algorithm was applied to hightemperature porous heat exchanger technology, which could be a potential solution for improved heat recovery. A Kelvin cell structure was used as the porous geometry and placed inside the heat exchanger. The heat recovered in terms of temperature rise was evaluated at different solid thermal conductivity and fluid inlet velocity parameters using the Monte Carlo method and compared to the deterministic method. The salient features of the Monte Carlo and deterministic methods were also discussed.

INTRODUCTION

Burning fossil fuels to meet energy demands releases vast amounts of greenhouse gases (GHGs) such as CO 2 , CH 4 , and NO 2 , causing egregious effects on the climate [START_REF] Kumar | Impact of covid-19 on greenhouse gases emissions: A critical review[END_REF]. GHG emissions are primarily generated by energy, industry, building, and transport sectors [START_REF] Lamb | A review of trends and drivers of greenhouse gas emissions by sector from 1990 to[END_REF]. As the global population continues to increase, energy demands also rise, increasing the likelihood of elevating GHG emissions. From 1990 to 2021, the warming effect of our planet's climate (Radiative forcing) due to long-lived GHGs rose by nearly 50 %, with carbon dioxide accounting for approximately 80 % of this increase. To reduce GHG emissions, with a primary focus on the energy sector, improving energy efficiency in hightemperature industries such as metallurgy, steel, cement, and glass is essential. The development of high-temperature heat recovery/storage/transport solutions is required to improve energy efficiency. In this work, we present a potential solution for heat recovery.

Shell-and-tube heat exchangers are widely used in process industries, as well as in many types of HVAC equipment, for high-temperature heat recovery at an industrial scale [START_REF] Trushad | Application of porous metal foam heat exchangers and the implications of particulate fouling for energy-intensive industries[END_REF]. Although stochastic foams are also used, they are less in demand. Deployment of new heat exchanger technologies could revolutionize the way industries use energy. Research on foam heat exchangers is steadily gaining momentum. To improve heat exchanger efficiency, the use of porous materials could be an attractive alternative to conventional heat exchanger technologies. To achieve better heat recovery, foams should be designed to have less pressure drop and a high temperature difference between the entry and exit of an element filled with such foams. However, optimizing these foams in terms of geometrical properties such as strut diameter, pore diameter, and foam design types like Kelvin-cell, structural, and stochastic foams, to obtain better heat recovery, is a challenging task.

Before performing optimization, it is essential to understand the physics and control the computational time. Deterministic methods, which are commonly used to solve coupled heat transfer and fluid flow * Corresponding Author: raj-narayan.konduru@univ-lorraine.fr through porous heat exchangers in 3D, use numerical approximations that consume significant computational time and memory. This computational effort is directly proportional to the complexity of the geometry. To address this problem, a semi-meshless method can be used for simulating coupled conduction-convection-radiation heat transfer. We use the finite volume method to solve the velocity fields and the Monte Carlo method [START_REF] John R Howell | Application of monte carlo to heat transfer problems[END_REF][START_REF] Mikhailovich | Monte Carlo Methods in Mechanics of Fluid and Gas[END_REF] to solve coupled heat transfer. The algorithm accounts for the coupling of conduction in the solid phase, convection, and radiation in the fluid/void phase. This approach is shown to use less memory and computation time for complex geometries.

PROBLEM DESCRIPTION

In this paper, we propose using porous structures, such as Kelvin-cell, structural, or stochastic foams, inside a heat exchanger to improve heat recovery. We consider these porous structures as complex geometries for numerical analysis, rather than solving the equations using porous properties. Specifically, we use the Kelvin-cell type of porous structure. The geometrical model of the Kelvin-cell foam was provided by CTIF (Centre Technique des Industries de la Fonderie), and samples made of this foam were manufactured in cast iron by CTIF. We conduct experiments in our lab using air as the fluid medium for convective analysis. For radiative analysis, we treat the fluid as a transparent medium and the solid as opaque. To solve for conduction, we use a solid medium with the properties listed in table 1. All properties mentioned in table 1 are constant, and temperature dependence has not been taken into account.

The solid medium corresponds to metallic materials such as cast iron, and its typical thermal properties are mentioned in table 1 except for λ s . The solid thermal conductivity λ s of cast iron is 40 W • m -1 • K -1 , but in our analysis, we also consider λ to be 1 W • m -1 • K -1 to decrease conduction and thus increase the role of radiation in heat transfer. In order to validate our proposed method, we present a case study as shown in fig. 1 only for better visual understanding. The foam is a structured lattice with a Kelvin-cell pattern, as shown in fig. 2. The cell diameter (d c ) is equal to 14 mm and the strut diameter (d strut ) is equal to 2.5 mm. The foam part has a length (l) of 200 mm, an inner diameter (d i ) of100 mm and its porosity (φ ) is 85 %. For the given operating conditions of maximum temperature (1073 K) and V in to be 2 m • s -1 , Knudsen number was calculated to be less than 0.01. This indicates that the fluid can be treated by the continuum model, which assumes that the fluid properties can be defined at every point in the fluid. Mach number was calculated to be 0.003, which indicates that the flow can be considered incompressible and the incompressible flow equations can be used.

V in = 2 m • s -1 Medium ρ C p λ ε µ (in kg • m -3 ) (in J • kg -1 • K -1 ) (in W • m -1 • K -1 ) (in Pa • s)

NUMERICAL MODELS

Numerical modelling is widely applied technique to tackle complex problems by computational simulation of such scenarios. It uses mathematical models to describe the physical conditions of the scenario using equations. Some partial differential equations are difficult to solve directly, these numerical models can approximate the solution. In this work, two types of numerical methods were used and are discussed in following sections.

3.1. Deterministic Method: These methods involve numerical approximations such as Finite Difference method (FDM), Finite Volume method (FVM), and Finite Element method (FEM). In this work, the Finite Volume method (FVM) was used to solve the coupled heat transfer problem using a Conjugate heat transfer model. ANSYS Fluent® software was employed for the simulation. It is important to note that only steady-state simulations were carried out in this study. This numerical simulation involves three steps.

Geometry. The Kelvin-cell structure, which serves as the solid medium is shown in fig.

2. An enclo- sure was developed around this structure to serve as the fluid medium, and both media were perfectly interconnected to form a perfect interface. The geometry was checked for errors and, once confirmed to be error-free, moved on to the next operation.

Meshing. Meshing is a process of dividing the volume of interest into smaller sub-volumes called elements. In this work, polyhedral elements were used as a mesh type for better stability when dealing with complex geometries [START_REF] Sosnowski | Polyhedral meshing in numerical analysis of conjugate heat transfer[END_REF], [START_REF] Wang | Comparison of hexahedral, tetrahedral and polyhedral cells for reproducing the wind field around an isolated building by les[END_REF]. Mesh errors were checked and a grid independence test was performed.

An element size of 0.5 mm and 1.0 mm were found to provide a stable solution, and the latter was chosen for further simulations due to its lower computational memory and time compared to the former.

Solver and Post-processing. In this study, the Fluent® solver was employed to simulate the heat transfer process. The Conjugate heat transfer model was utilized for both constant density and perfect gas density cases. To account for radiation heat transfer, the surface-to-surface radiation model was employed. The results obtained from the simulation were then post-processed and discussed in further sections.

It should be noted that the simulation was performed on the full 3D geometry. The geometry is not axially symmetric, but has to two planes of symmetry. The entire structure was modeled in our study to test its capability in handling complex geometries. Indeed, our aim is to develop a model that can work for any complex and even stochastic structures. For testing purposes, we considered a Kelvin cell structure. However, for non-structured foams, symmetry boundary conditions cannot be applied, and the entire domain must be computed. Hence, it is better to compare the computational resources required for whole domains and complex geometries to evaluate the model's capabilities.

Semi Mesh-less Method:

This method is a combination of deterministic methods using FVM (Mesh-based approximation) to obtain velocity fields and a stochastic method, i.e., Monte Carlo method (Mesh-less), to solve the coupled heat transfer problem, as proposed by Loris Ibarrart et al. [START_REF] Ibarrart | Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm : getting insensitive to geometrical refinement[END_REF], to obtain the temperature at any prescribed point in the domain. As per our application, we could calculate the temperature at the outlet directly using this method instead of calculating the whole field as performed The Monte Carlo method only requires the surface mesh of the interface separating the solid medium from the fluid medium, whereas deterministic methods need a volumetric mesh of both media. The random walk technique is the main principle behind this method. To compute the temperature at location x, a large number of paths crossing the complex geometry are calculated. The temperature is the expectation of a random variable T(x). This random variable is determined from a path inside the calculation domain. A path starts at the location and at the time where and when the temperature needs to be computed. Each path can follow a conductive path in the solid (considered here opaque) or a radiative path in the void phase (considered here perfectly transparent) or a convective path. The walk/path step size and direction are based on random number generators and probabilities derived from partial differential equations. Some famous random walk techniques are shown in fig. 3 for one realization only. Multiple realizations are required to initialize, and then the expectation is computed. It is suggested that the reader refers to [START_REF] Ibarrart | Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm : getting insensitive to geometrical refinement[END_REF] for a detailed mathematical formulation and analytical validation of the coupled algorithm.

For the sake of brevity, we will not discuss the derivation in this paper. For better understanding, the algorithm is described in fig. 4.

Fig. 3: Random walk techniques

Yes

Fig. 4: Semi Mesh-less Algorithm Conduction is solved by using the floating random walk method in the solid medium as shown in fig. 3(a) [START_REF] Haji | The floating random walk and its application to monte carlo solutions of heat equations[END_REF]. As air is considered a transparent medium, the surface-to-surface interaction model is used to calculate radiation in the fluid medium [START_REF] Michael | Modest[END_REF]. Convection is solved by transforming the energy equation to a random walk on the grid in the fluid domain as shown in fig. 3(b). To solve the energy equation in eq. ( 1), prior knowledge of the velocity field is required, as the path moves randomly through the domain, as shown in fig. 3. Thus, an interpolation technique is needed to estimate velocities between grid points. To handle complex data, we employed scattered data interpolation [START_REF] Amidror | Scattered data interpolation methods for electronic imaging systems: a survey[END_REF]. We performed a validation test on a 3D complex geometry, specifically the flow through a Kelvin-cell foam placed at the center of a long tube with a diameter of 100 mm and a length of 800 mm. The inflow velocity V in was set to 0.1 m • s -1 , and the fluid medium was air with a Reynolds number of approximately 660. We validated our interpolation technique by comparing the velocity profile along the horizontal axis with the expected results, as shown in fig. 5b. For reference, the velocity contour for this test case is shown in fig. 5a. The results of the interpolation technique were found to be in good agreement with the results obtained from the deterministic method (Fluent®). However, it should be noted that an increase in the amount of data results in longer computational times. In comparison to the radial basis interpolation technique [START_REF] Amidror | Scattered data interpolation methods for electronic imaging systems: a survey[END_REF], this method is preferred as it requires less computational memory. The radial basis technique needs to calculate the Euclidean distance between each point for interpolation, which results in the development of a large matrix. The second step involves discretizing eq. ( 1) using FDM and determining direction probabilities to enable the random walk. This is a type of fixed random walk.

∇.(-λ f ∇T (x) + ρ f C p f VT (x)) = 0, x ∈ Ω f (1)
The rearranged equation is written as, T (x, y) = P x+ T x+ + P x-T x-+ P y+ T y+ + P y-T y-+ P z+ T z+ + P z-T z-

where,

P k+ = α f δ 6α f δ + u + v + w P k-= α f δ + j 6α f δ + u + v + w k ∈ [x, y, z] j ∈ [u, v, w] resp. α f = λ f ρ f C p f
where k is location coordinate vector and j is the velocity vector. The finite difference method formulation (Method-1)) and Peclet formulation (Method-2) were both applied to the energy equation in eq. ( 1), as suggested by [START_REF] Ibarrart | Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm : getting insensitive to geometrical refinement[END_REF]. Both methods were found to be in good agreement. Method-2, which uses Peclet formulation, was found to be more accurate than Method-1, which uses FDM, by a margin of 0.1 % when compared to the finite volume method (FVM). Nonetheless, Method-2 requires more computational time than Method-1. Therefore, we used Method-1 in this study based on a trade-off between computational time and accuracy.

VALIDATION

Using both deterministic and semi mesh-less methods, simulations were performed by considering inlet ambient temperature of air at 300 K, T imp on the solid outer surface at 1073 K, air V in at 0.1 m • s -1 , 0.2 m • s -1 and 0.5 m • s -1 and λ s at 1 W • m -1 • K -1 and 40 W • m -1 • K -1 to simulate the case. No gravitational effects were considered. Near central horizontal profile as shown in fig. 1 is used to draw the plots. This profile is only located in fluid medium. Monte Carlo algorithm results are compared to Fluent® results, the latter being considered the reference solution. As per the fundamentals of this algorithm, velocity fields is used to compute the temperature but in return, no effect of temperature on velocity fields was considered (one-way coupling). So, this algorithm is equivalent to the constant density deterministic approach. Temperature in the domain was computed in deterministic approach using both variable density (perfect gas law) and constant density.

The objective to is study the performance of the algorithm and compare results generated by deterministic approach (Fluent®). Primarily results obtained by Fluent® will be discussed followed by comparison with Monte carlo algorithm. As per Mach number, we have assumed the flow is incompressible. Using Fluent®, we have performed the simulations for velocity magnitude and temperature for λ s at 1 W • m -1 • K -1 and V in at 0.1 m • s -1 and 0. It was observed that a significant difference in velocity magnitude at lower V in but no difference in temperatures for both chosen inlet velocities by keeping density as constant (DC) and variable (DV). For the case considering variable density, increase in velocity is due to the thermal expansion of the gas. It is lower at high V in because the temperature rise is lower and it is due to the energy that is transferred by conduction from the outer part of the tube to the inner part and the amount of energy transported is limited by the thermal conductivity. If the gas flow rate increases, the temperature rise should decrease.

λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.1 m • s -1
Indicates for line plots Fluid flow 0.3 0.5

(b) Velocity contour (DV) :

λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.1 m • s -1
(c) λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.1 m • s -1
(d) λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.1 m • s -1
(e) λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.5 m • s -1
f) λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.5 m • s -1 (
Fig. 6: Variation in temperature and velocity magnitude keeping density as constant and variable using deterministic method.

Because constant density and variable density are equivalent from the thermal point of view, for further analysis, we will consider the simulations of constant density as reference solution for comparing Monte Carlo algorithm.

This algorithm was validated against Fluent® results. The algorithm used 3000 realizations to compute the temperature at each chosen probe point which belongs to fluid medium only as shown in fig. 1. All results were at 1 % standard deviation. But in fig. 7, we consider ±2 σ interval in plot representations to have a better confidence (i.e. 95% instead of 67%). As per our application, the main interest is all about the temperature rise between inlet and outlet of the tube. Based on this, L1 error was computed between DM results and MC results using eq. ( 3). (a) (b) As from figs. 7a to 7d, the trend in maintained through out foam length. Since our goal is to compute the outlet temperature of the foam, eq. ( 3) was used to calculate the error at the outlet section only. For λ s at 1 W • m -1 • K -1 , L 1 error was found to be 2 %, 3 % and 5 % for V in at 0.1 m • s -1 , 0.2 m • s -1 and 0.5 m • s -1 respectively. Similarly at λ s at 40 W • m -1 • K -1 , L 1 error was found to be 2 % for V in at 0.1 m • s -1 . It is also observed that changing λ s does not lead to a significant change in the error %, but changing V in increases the error. It was also observed that increases the λ s also increase in temperature at the outlet. In general, this algorithm has good agreement with low input speeds and the percentage error increases with increasing input speed. The reason could be that increasing the input velocity in complex geometries leads to large fluctuations in temperature and velocity, which makes the algorithm less efficient to predict accurately. Secondly, It was ensured that the Deterministic Method (DM) solution has converged and well approximated but it is not the exact solution.

L 1 = |T DM | -|T MC | |T DM | (3)
λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.1 m • s -
λ s = 1 W • m -1 • K -1 , T imp = 1073 K, V in = 0.2 m • s -1
Table 2 summarizes the features of Monte Carlo algorithm compared with deterministic method. In this application, the objective is to quantify the foam efficiency which depends both on the temperature rise along the foam and the flow rate. In the DM method, the entire temperature field needs to be calculated for every change in parameter, which can be computationally expensive and time-consuming. However, with the proposed method, we only need to obtain the outlet temperature at a few selected points, which significantly reduces the computational resources required. In contrast, the DM method calculates the whole field, which is not necessary for our objective of optimizing the porous structures. Furthermore,

• m -1 • K -1 and 40 W • m -1 • K -1
, we can use the same 0.1 m • s -1 velocity field data from the DM method as the input for Monte Carlo simulation, and then we can obtain the temperature for both solid thermal conductivity's without recalculating the velocity field, which saves a considerable amount of time and resources.

CONCLUSIONS

Semi mesh-less method was applied to porous heat exchanger at high temperatures and compared with deterministic method at stationary state. This algorithm was then compared with results obtained by FVM method. Results were discussed along with the reasons of agreements and disagreements. Linearization of radiative term is one of the main assumption of this algorithm. In occurrence of high non-linearity situations, this assumption would not be valid. The mesh-less method could be very interesting to use in terms of reducing computational time and memory while maintaining the accuracy given by FVM methods. Future is development of full mesh-less method because for now the calculation is still limited by the CFD calculation for velocity fields. Porous heat exchanger experiments will be performed as an upcoming work. The experimental results will be used to validate the model discussed in this work.

The transient behaviour of these system will be recorded experimentally to analyse the storage capacity.

ACKNOWLEDGEMENT

The Monte Carlo algorithm was implemented using Star-Engine libraries developed by Meso-star Company https://gitlab.com/meso-star/star-engine.git. Région Grand Est is acknowledged for funding Raj Konduru Phd Thesis. This work is part of the THERMOUSSE project funded by ICEEL. This work was performed using HPC resources from the EXPLOR mesocentre of University of Lorraine. https://explor.univlorraine.fr/

NOMENCLATURE

Fig. 1 :

 1 Fig. 1: Porous Heat exchanger -Sample simulation Temperature contour : 300 K(Blue color) ; 1073 K (Red color) for V in = 2 m • s -1

  Fig. 2: Geometry Model

1 Fig. 5 :

 15 Fig. 5: Scattered data interpolation test.

  5 m • s -1 as shown in figs. 6c to 6f. The referenced velocity contours for fig. 6c are shown in fig. 6a and fig. 6b.

  Velocity contour (DC) :

1 Fig. 7 :

 17 Fig. 7: profile computed using Monte Carlo algorithm compared deterministic method.

αρu

  Thermal diffusivity [m 2 • s -1 ] Density [kg • m -3 ] V Velocity [m • s -1 ] C p Heat Capacity [J • kg -1 • K -1 ] Velocity component in x-axis direction [m • s -1 ] v Velocity component in y-axis direction [m • s -1 ] w Velocity component in z-axis direction [m • s -1 ]x Coordinate x-axis [m] y Coordinate y-axis [m] z Coordinate z-axis [m]

Table 1 :

 1 Solid and Fluid media properties

	Solid	8000	450	1 & 40	0.9	-
	Fluid	1.225	1005	0.025	-	1.81 * 10 -5

  1

		1000	DM -DC MC
			MC ±2𝜎
		900
	Temperature (in K)	600 700 800
		500
		400
		0.300 300	0.325	0.350	0.375	0.400	0.425	0.450	0.475	0.500
					Position (in m)

Table 2 :

 2 Features of Monte Carlo Method over Deterministic method for every change in parameter, the DM method requires calculating the velocity field from the beginning. For instance, if we need to obtain the temperature for two different solid thermal conductivity's, 1 W

	Type	Monte Carlo Algorithm	Deterministic method
	Geometry input file	Surface geometry (.obj)	Solid geometry (.stp)
	Mesh type	Surface mesh	Volume mesh
	Geometry complexity effect on memory	Consumes less memory	Consumes more memory
	Results obtained	Probe locations	Total field
	Computational time	45 min for each plot in fig. 7	6 hours (Steady State)
		150 probe points for each plot
	Parallelization	Yes	Yes
	No. of processors used	32	32
	Input data	Velocity fields required	Not required ; V & T are calculated


L Penazzi

O Farges

Y Jannot

J Meulemans

V Schick

RESOLUTION OF CONDUCTIVE-RADIATIVE HEAT TRANSFER IN A SEMI-TRANSPARENT MEDIUM: FUNCTIONAL ESTIMATION OF THE RADIATIVE SOURCE TERM WITH MONTE CARLO METHOD

Keywords: Conduction, Radiation, Semi-transparent medium, Monte Carlo method, Transfer function

To estimate materials' thermo-physical properties at high temperatures, there is a need of robust direct models that can be solved quickly and account for the strong coupling between conduction and radiation. This study investigates using the Monte Carlo method to solve a conduction-radiation model in a semitransparent material in a reduced computational time. Transient conduction in an absorbing/emitting and non-scattering grey medium in an academic configuration is investigated. The radiative source term of the heat equation is estimated as a function of the temperature field to the power 4, using a Finite Differences algorithm and a single Monte Carlo simulation to solve radiation. The numerical study shows that the results are in good agreement with the literature and that the functional estimation of the radiative source term allows coupled model resolution in less computing time that the standard approaches, such as Finite Volumes or Elements Methods. Indeed, the coefficients of the function were obtained in 4h with the Monte Carlo computation on an academic calculation server. 15 different configurations were then solved in less than 1min in transient state configurations and less than 6min at steady state configurations with the Finite Differences Method on a conventional laptop. This function is obtained without any compromise on the complexity of the coupled physics or the geometry. Because of the latter features, using this method in inversion processes becomes easible, as the direct model must be evaluated a large number of times.

INTRODUCTION

To optimize the use of insulating materials in high temperature processes, developing measurements methods to characterize the thermal properties of these materials at their operating temperature (beyond 1000 • C) is crucial [START_REF] Machin | The European project on high temperature measurement solutions in industry (HiTeMS) -A summary of achievements[END_REF]. One of the main challenges lies in the development of direct models taking into account the strong conduction-radiation coupling in the sample at high temperature [START_REF] Penazzi | Influence of radiation heat transfer on parallel hot-wire thermal conductivity measurements of semi-transparent materials at high temperature[END_REF][START_REF] Zhang | Numerical study of the influence of thermal radiation on measuring semi-transparent thermal insulation material with hot wire method[END_REF]. The heat transfer by radiation in a participating medium (emitting/absorbing and scattering) needs to be considered as well as the potentially complex geometry of the experimental set-up and the sample. Moreover, in order to be used in an estimation methodology, this coupled model needs to offer a decent compromise between accuracy and speed of execution during model resolution.

So far, in most of the inverse processes associated to the thermal properties estimations (thermal conductivity or diffusivity, for instance), either purely conductive model are used or coupled models using assumptions on the optical properties (no scattering events, grey medium etc.) [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF][START_REF] Jannot | Thermal diffusivity measurement of insulating materials at high temperature with a four-layer (4L) method[END_REF]. Hence, a robust coupled direct model in which these assumptions may be lifted and without limitation to asymptotic approach (optically thick or *Corresponding Author: lea.penazzi@univ-lorraine.fr thin media) is needed. Previous studies demonstrated the feasibility and relevance of using the Monte Carlo method to solve coupled and transient thermal problems in the context of complex media characterization [START_REF] Sans | Experimental characterization of the coupled conductive and radiative heat transfer in ceramic foams with a flash method at high temperature[END_REF].

In the present work, the development and resolution of a coupled conduction-radiation model in a semitransparent media is addressed. The studied configuration is a standard academic case corresponding to a two-dimensional square enclosure containing an absorbing/emitting and non-scattering grey medium bounded by black walls. Heat conduction is considered at transient state and is coupled to volumetric radiation. When Mishra et al. [START_REF] Mishra | Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem[END_REF] investigated this transient problem, the heat transfer equation was solved using alternatively a Finite-Volume Method (FVM) and a Lattice Boltzmann Method (LBM), while the radiative source term was calculated using the the Collapsed-Dimension Method (CDM). This reference case was as well studied by Asllanaj et al. [START_REF] Asllanaj | Transient Radiation and Conduction Heat Transfer in a Gray Absorbing-Emitting Medium Applied on Two-Dimensional Complex-Shaped Domains[END_REF]. They used the FVM to solve the radiative source term based on a cell vertex scheme and associated to a modified exponential scheme. The heat transfer equation was solved using low or high order finite elements. Their results were in agreement with Mishra et al. [START_REF] Mishra | Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem[END_REF] results.

In this study, the transient conduction is solved by Finite Differences Method (FDM) and the volumetric radiation is solved by the Monte Carlo Method (MCM). The originality and interest of this approach lies in the functional estimation of the radiative source term with the MCM allowing to evaluate the temperature field in a highly reduced computational time compared to conventional methods. This functional estimation is carried out without any compromise neither on the complexity of the coupled physics nor the geometry. The purpose of this study is to present a proof-of-concept of the method on an academic case that can be expanded to more complex cases in the future. First, the model of the academic configuration is presented, then, the algorithm of resolution is introduced including the functional estimation of the radiative source term with the MCM coupled with the FDM resolution. Finally, the results are presented, compared to the numerical solutions of Mishra et al. [START_REF] Mishra | Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem[END_REF] used as a benchmark and conclusions are drawn.

THE FUNCTIONAL ESTIMATION MODEL

Studied configuration

The studied configuration is illustrated on Fig. 1, which corresponds to a square enclosure of unit length ( =1 m) containing a non scattering gray medium ( s = 0). The optical thickness = = ( a + s ) is equal to 1, as the absorption coefficient a is equal to 1m -1 , and the medium boundaries are black surfaces ( = 1) with imposed temperatures. The reference temperature hot is the temperature of the bottom wall  hot . The other surface temperatures are cold and equal to cold = hot ∕2 (on the surface  cold ). Heat conduction is at transient state and at initial time, the temperature of the medium is set equal to cold . The heat transfer model equations are the following :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ p - = -⃗ ∇ ⋅ ⃗ R ∀⃗ ∈ , ∀ > 0 (⃗ , 0) = cold ∀⃗ ∈ , = 0 (⃗ , ) = hot ∀⃗ ≡ ⃗ hot ∈  hot , ∀ > 0 (⃗ , ) = cold ∀⃗ ≡ ⃗ cold ∈  cold , ∀ > 0 (1) (2) (3) (4)
The radiative source term from Eq. (1) may be obtained from a radiative balance of a homogeneous and isothermal volume  expressed as follows:

 = ∫  ( ⃗ ∇ ⋅ ⃗ R )  = ∫   ∫ 4 ⃗ a  eq (⃗ ) -(⃗ , ⃗ ) (5)
The medium being grey, the optical properties and the radiative intensity do not depend on frequency. In this equation, ⃗ corresponds to an elementary volume element and ⃗ , to the solid angle. The black body intensity  eq at the location ⃗ is expressed as:

 eq (⃗ ) = SB 4 (⃗ ) (6)
With SB , the Stefan-Boltzmann constant and (⃗ ) the temperature at the ⃗ location. (⃗ , ⃗ ) is the radiative intensity at the location ⃗ in the direction ⃗ and is solution of the following Radiative Transfer Equation (RTE):

⃗ . ⃗ ∇(⃗ , ⃗ ) = -( a + s )(⃗ , ⃗ ) + a  eq + s ∫ 4 P(⃗ , ⃗ | ⃗ ′ )(⃗ , ⃗ ′ ) ⃗ ′ (7)
Where ⃗ is the unit direction, s , the scattering coefficient and P, the scattering phase function. The refractive index r is equal to 1.

Heat conduction numerical resolution

The temporal term , is discretized with a explicit first- order Euler method and the laplacian term is approximated with a second-order difference approximation. The heat equation corresponding to Eq. ( 1) becomes :

n i,j -n-1 i,j t - n-1 i-1,j + n-1 i+1,j -4 n-1 i,j + n-1 i,j-1 + n-1 i,j+1 2 x = - 1 p  n-1 ray,i,j (8)
With the spatial step x (identical for both dimensions, i.e. x = y ), the time step, t and the thermal diffusivity, = p . As represented in Fig. 1, the 2D medium is discretized into a number × of conductive meshes. The conduction mesh is cartesian regular with = . The conduction meshes are indexed with ranging from 0 to -1 and ranging from 0 to -1. Each node [ , ] (represented using purple circles on Fig. 1) is located at the center of a conductive mesh and each mesh has a surface equal to 2 x . In the finite differences method, the assumption made is that the temperature of each conductive mesh is homogeneous. Reformulating Eq. ( 8), leads to the following expression for the temperature at a position ( , ) and time :

n i,j = 1 -4 t 2 x n-1 i,j + t 2 x n-1 i-1,j + n-1 i+1,j + n-1 i,j-1 + n-1 i,j+1 -t p  n-1 ray,i,j (9)
In this case, the Courand-Friedrichs-Lewis (CFL) condition states that ≤ 1∕4 so that this explicit scheme remains stable. The previous Eq. ( 9) remains valid for ∈ { 1 ∶ -2} and ∈ { 1 ∶ -2}. For the boundary nodes located either on the cold surfaces (when ⃗ cold ∈  cold , see Fig. 1) on the hot surface (when ⃗ hot ∈  hot , see Fig. 1), the temperature values are set respectively to cold and hot .

The volumetric source term  ray (in W m -3 ) is computed at each iteration of the finite differences algorithm for each conductive node [ , ]. To this end, the radiative source term  is estimated for each radiative node and then equally allocated between the four conductive nodes nearby, such as :

 n-1 ray,i,j = n-1 ,i-1,j-1 + n-1 ,i-1,j + n-1 ,i,j + n-1 ,i,j-1 4 2 x ( 10
)
The radiative mesh is chosen slightly shifted from the conductive mesh and in each radiative mesh, the radiative source term  is homogeneous. In Fig. 1, the ( -1) × ( -1) radiative meshes (represented in green squares) are indexed with from 0 to -2 and from 0 to -2.

Estimation of the radiative source term with the MCM For each of the ( -1)( -1)

radiative node, a Monte Carlo computation is required to estimate the radiative source term  . The expression of the source term in Eq. ( 5) needs to be translated in statistical terms to estimate it with a Monte Carlo algorithm. Given the expression of the solid angle ⃗ = sin depending on the polar and azimuthal angles denoted respectively and , the radiative source term may be expressed as :

 = ∫  X (⃗ ) ⃗ ∫ 2 ( ) ∫ 0 ( ) 4 2 x a  eq -(⃗ , ⃗ ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ Monte Carlo weight k (11)
The previous equation can be given a statistical interpretation by introducing the probability density functions. The probability density X (⃗ ) = 1 = 1 2 x corresponds to the sampling of a starting position in the 2D radiative mesh and ( ) = sin 2 and ( ) = 1 2 are the probability density corresponding respectively to the sampling of the polar and azimutal angles necessary to the unit direction ⃗ sampling. The previous Eq. ( 11) may be interpreted as the expectation of a random variable , the Monte Carlo weight, such as  = ( ). A Monte Carlo computation consists in sampling MC optical paths, providing MC realizations 1 , 2 ... MC of , and estimate the radiative source term  as the mean value of these Monte Carlo weights realizations, such as  ≈ = 1 MC ∑ MC =1 k as well as the associated uncertainty

= 1 √ MC 1 MC ∑ MC =1 2 k -2 1 2
. The essence of the path sampling, represented on the Fig. 2 is as follows.

The heat flux  in W estimated in Eq. ( 5) becomes surfacic and expressed as follows

 = ∫  ( ⃗ ∇ ⋅ ⃗ R )  = ∫   ∫ 4 ⃗ a  eq (⃗ ) -(⃗ , ⃗ ) in W m -1 .
A start position ⃗ start is sampled uniformly within the medium  according to X (⃗ ). The current position ⃗ c value is initialized at the value ⃗ start . Then, a direction ⃗ is sampled with the successive sampling of the angles and according to the probabilities densities ( ) and ( ). An extinction length is sampled according to the exponential Beer-Lambert law (depending on the extinction coefficient = a + s ) and the current position ⃗ c value is actualized to ⃗ c + ⃗ . A test is performed to identify if the current location ⃗ c is stil in the volume  or at one of the medium boundary, either  cold or  hot . If ⃗ c ∈ , an absorption event occurs and the path stops at the location ⃗ c (see 1 in Figs. 2 et 3). The path's end location ⃗ end is set equal to ⃗ c and the Monte Carlo weight is stored. If ⃗ c is not in the volume , a test is performed to know at which boundary is the current location. An absorption event occurs either at the  hot boundary and the path stops at the location ⃗ c ≡ ⃗ hot (see 2 in Figs. 2 et 3) or at the  cold boundary and the path stops at the location ⃗ c ≡ ⃗ cold (see 3 in Figs. 2 et 3), then the corresponding Monte Carlo weight is stored. As indicated in Eq. ( 11), the Monte Carlo weight k stored at the end of each path contains the difference between the emitted radiant intensity  eq and the absorbed radiant intensity (⃗ , ⃗ ) for a given location ⃗ and direction ⃗ . Translated in algorithmic terms, k contains the difference between the black body intensity emitted in a mesh [i start , j start ] and the black body intensity absorbed in a mesh [i end , j end ], i.e. k = 4 2

x a  eq [i start , j start ] - eq [i end , j end ] .

The corresponding functional estimation of the radiative source term

The radiative source term  needs to be evaluated for each of the ( -1)( -1) radiative mesh at each iteration on the Finite Differences (FD) algorithm. At each calculation of the source term, MC optical paths are sampled. These operations may be time consuming and can be questioned when one notices that at each iteration of the FD algorithm, only one set of values change in input of the MC algorithm : the list of temperature at each node to the power 4. Given the fact that each time the ( -1)( -1) radiative source terms

 [i,j] ≡  [i start , j start ]
the only values changing in input are the equilibrium radiative intensities  eq computed, i.e. the list of temperatures [i,j] 4 ≡ [i end , j end ] 4 . This can be noticed by injecting Eq. ( 6) in Eq. ( 11) and express the estimated radiative source term at the [i start , j start ] as follows :

 [i start , j start ] ≈ 1 MC MC ∑ =1 4 SB 2 x a [i start , j start ] 4 -[ end, , end, ] 4 ⇔  [i start , j start ] ≈ 4 SB 2 x a ( [i start , j start ] 4 - 1 MC MC ∑ =1 [ end, , end, ] 4 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ A ) (12)
The second part A of the previous equation can be reformulated as :

A = 1 MC MC ∑ =1 ( end, = 0; end, = 0) [0, 0] 4 + ... + ( end, = -1; end, = -1) [ -1, -1] 4 (13)
In which  is the Heaviside function equal to 1 if the condition is satisfied and 0 otherwise. For the first term, at the end of the th path, if end, = 0 and end, = 0,  is equal to 1.

1

MC MC ∑ =1 [ end, , end, ] 4 = 1 MC 0,0 [0, 0] 4 + ... + -1, -1 [ -1, -1] 4 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟ ∑ -1 =0 ∑ -1 =0 , MC [ , ] 4 (14)
In this way, the latter expression can be transformed so that at the end of the MC calculation, counters stored the number of times a path ended in this or that mesh [ , ]. For instance, 0,0 stores the number of path that ended in mesh indexed [0, 0] out of a total of MC paths.

Hence, replacing the latter expression in A from Eq. (12) leads to the following expression :

 [i start , j start ] ≈ 4 SB 2 x a ⎛ ⎜ ⎜ ⎝ [i start , j start ] 4 - -1 ∑ =0 -1 ∑ =0 , MC [ , ] 4 ⎞ ⎟ ⎟ ⎠ (15)
In order to estimate the  [i start , j start ] as a function of the whole temperature field, a single Monte Carlo computation may be performed in which, instead of storing the Monte Carlo weight as expressed previously, the only coefficients needed to be stored are each , corresponding to the number of times the radiative paths have been absorbed in the mesh [ , ] during the Monte Carlo calculation (including the mesh boundaries).

Concretely, before starting the finite differences algorithm, ( -1)( -1) Monte Carlo computations are perfomed only once to estimate each  [i,j] as a function of the temperature list at each conductive node to the power 4, such as  [i,j] = f ( n-1 i,j ) 4 . For each of the ( -1)( -1) radiative source terms

 [i,j], at each Monte Carlo computation, × coefficients ,

MC

are computed. These , MC coefficients may be assimilated to volumetric view factor informing on how each radiative mesh element sees another. Similarly as previously, the associated uncertainty of each Monte Carlo computation may be estimated. In this way, the temperature field to the power 4 corresponding to the matrix ( n-1 i,j ) 4 is injected in the function at each iteration of the FD algorithm to estimate the  [i,j]. indefinitely. The hot temperature is considered as the reference temperature. The other faces of the square enclosure are kept constant at cold during this time.

To quantify the relative contribution of conduction versus radiation, a dimensionless parameter is introduced, the Stark number (also named as the conduction-radiation parameter) expressed as  = . For large values of  , conduction predominates over radiation and, on the opposite, for small values of  , radiation predominates. Three cases are presented in which  is set equal to 0.01, 0.1 and 1. In order to compare these cases a dimensionless time parameter is introduced : = 2 a . The results are compared at different instants = { 0.001, 0.005, 0.015, 0.04} and three differents values corresponding for each to the steady state, = 0.06, 0.17, 0.2 for respectively  = 0.01, 0.1, 1. For these given values of the Stark numbers and dimensionless time parameters, the temperature profiles on the centerline ( = ∕2 and ∈ [0, ]) are computed with the FEMC/FD algorithm and presented on Fig. 4 . Each of these temperature profiles (denoted with red-filled bullets) are compared to results from two different numerical methods from Mishra et al. [START_REF] Mishra | Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem[END_REF], i.e. the FVM method (represented with grey-filled triangles) and the LBM method (represented with blue continuous lines). The results from the FEMC/MC method are in good agreement with results from the reference case.

As presented in Tab. 1, one of the major benefit of the FEMC/FD algorithm is to enable a significant computational time reduction compared to the original MC/FD algorithm. Indeed, to produce this whole set of results for differents values of Stark number and dimensionless time step, only a single FEMC computation was necessary to estimate the function's coefficient and the calculation lasted 4 h (see Tab. 1) for MC = 10 000 000 and a total of 21 × 21 conductive nodes. After this, each of the 15 graphs (3 values of  , for each 5 values of ) were computed with the FD algorithm in less than approximately 1min for transient state and a maximum value of 340s at steady state (see = 0.2 and  = 1 on Tab. 1).

CONCLUSIONS

The present study was a proof of concept in which the radiative source term of the heat equation was evaluated as a function of the temperature field to the power 4 without linearization of the radiative heat transfer. This feature proved to be interesting in order to reduce the computation time required for the resolution of the coupled conduction-radiation heat transfer in a semi-transparent medium at transient state. Indeed, this function was obtained without any assumption on the radiative behaviour of the physics and is not limited to grey media, the dependency of the radiative properties on the frequency may be taken into account. The scattering character of the medium will be studied in further configurations and does not pose any particular difficulties as these may be compared to literature results from [START_REF] Mishra | Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem[END_REF]. Furthermore, the development of this method is not limited to academic configuration and may be extended to more complex geometries (2D/3D). For all of these reasons, to consider the development and implementation of the FEMC/FD method in inversion processes becomes practicable as the direct model needs to be 1) as realistic as possible and 2) a compromise between precision and being solved within a reasonable time. Further developments may be continued to implement this direct coupled model in a practical configuration in order to carry out a sensitivity study to the various influencing parameters before the inverse process. The method will be developed specifically for use in inversing parallel hot wire measurements [START_REF] Penazzi | Influence of radiation heat transfer on parallel hot-wire thermal conductivity measurements of semi-transparent materials at high temperature[END_REF].

Fig. 1

 1 Fig. 1 Schematic view of the 2D configuration defining the academic thermal problem from [7, 8]. Conductive nodes are represented with purple circles and radiative nodes, with green squares.

X 2 3 Fig. 2

 232 Fig. 2 Path sampling algorithm corresponding to a th realization of the Monte Carlo algorithm to estimate  .

Fig. 3

 3 Fig. 3 Examples of random paths inside the square enclosure according to the algorithm from Fig. 2.

  hot =100 K, as well as a =1 m -1 , if the value of the Stark number  is set, the thermal conductivity value may be calculated with the following expression =

Fig. 4

 4 Fig. 4 Dimensionless temperature distribution along the centerline position at different dimensionless times and for three different values of the Stark number  .

Table 1

 1 Computation times (in seconds) for the pure absorbing case with black walls.

		Stark number  (-)
	Dimensionless time (-) 0.01 0.1	1
	0.001	1.7
	0.005	8.2
	0.015	24.4
	0.04	65.9
	0.06	99.7
	0.17	284
	0.2		340
	1 The CPU time to obtain the functional Monte
	Carlo coefficients is 4 h on an Intel®Xeon™E5
	2.10 GHz
	2 Computing times are given for an
	Intel®Core™i5-10210U CPU 1.60GHz

ACKNOWLEDGMENT

High Performance Computing resources were partially provided by the EXPLOR centre hosted by the University of Lorraine (Project: 2020EMPPX2210).

Finite differences and Monte Carlo algorithm

Calculation of each ⃗ new [i,j] (see Eq. ( 9)) end end ← + t end As described in the Algorithm 2.1, first, for each of the ( -1)( -1) element of the radiative mesh, all of the ( -1)( -1) , MC coefficients are estimated with the MC algorithm. This way, the radiative source term may be estimated afterwards as a function of the temperature field to the power 4. This method will be referred hereafter as the Function Estimation Monte Carlo (FEMC) method. Then in the FD algorithm, the ⃗ old vector is firstly initialized with its values set equal to ⃗ new , corresponding now to the (n -1) nth iteration (i.e. n-1 i,j ). Then, for each radiative mesh element [i,j], each source term  [i,j], is estimated thanks to the functional relation presented previously and the coefficients estimated with the FEMC method. Afterward, for each conductive node [i,j] the radiative source term  ray [i,j] is computed. Eventually, the corresponding ⃗ new [i,j] profile at the (n) nth time step of the finite differences algorithm can be computed (i.e. n i,j ). The algorithm proceeds to the next time step by adding t to the current time until the maximum time max is reached. If max corresponds to the steady state time a tolerance criterion is set, the ratio

should be lower than 1 × 10 -5 .

RESULTS

The pure grey absorbing medium  with black walls presented in Section 2.1 is studied (see Fig. 1). The medium is initially at a temperature cold =50 K. At a time > 0, the temperature on the bottom wall ( = 0) rises from cold =50 K to hot =100 K within one dimensionless time step and then remains constant