École Sigle

Doctorale

Renée El

Melhem Bât

Pascal Blaise

R Gourdon

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

M Fabrice

Mme Emmanuelle
email: emmanuelle.canet@univ-lyon1.fr

Étage

M Luca

Jean-Yves Buffière
email: jean-yves.buffiere@insa-lyon.fr

M Philippe Boisse

Viviane Polsinelli
email: viviane.polsinelli@univ-lyon2.fr

Brigitte Dubois

J Y Toussaint

M Christian Montes

Prof K S Wilson

Prof V Chiu

Naso

Rémi Coquard

Cellular media such as plastic, ceramic and metal foams present specific characteristics that make them interesting for a number of applications related to thermal engineering. Their ability to minimize natural convection makes them ideal candidates for insulation applications, and while plastic foams are already widely used in low temperature applications, low conductivity ceramic foams are increasingly used in high temperature applications. On the opposite, the high specific surface and permeability to fluid of open cell foams makes them interesting heat transfer enhancers, using metal or high conductivity ceramics as material. In addition, their permeability to light makes them an ideal candidate for thermal radiation based applications, such as porous burners or solar energy collectors.

In many of these application, thermal radiation heat transfer can have a significant influence on the heat transfer process, ranging from a 20-30% effect at room temperature to becoming the predominant mode of heat transfer at higher temperatures. To model radiation heat transfer in these media, both accurate radiation models and accurate morphological models of the structure of the foam are required. This work provides an original contribution on both these accounts.

An extensive discussion of the literature on numerical methods for radiation heat transfer in cellular media is presented, with a special focus on Monte Carlo methods. Homogeneous Phase (HPA) and Multi Phase (MPA) methods are discussed. Further efforts are required to accurately model and digitally replicate of foam morphologies. Our goal is to propose a unified framework that allows, within a single chain of tools, to digitally generate three commonly occurring types of foam structures, covering a large range of real materials: highporosity open cell foams, high-porosity closed cell foams, low-porosity open-cell structures. For highporosity open cell foams, a novel automated parametric digital generation technique was developed based on phenomenological investigation of the foam's morphological parameters, and validated against a dataset consisting of raw morphological data obtained by tomographic analysis. The generation capabilities were then applied to parametrically investigate the influence of morphological

Résumé Etendu

Les matériaux cellulaires macroporeux tels que les mousses polymère (Fig. 1), carbone (Fig. 2), céramique (Fig. 3) ou métallique (Fig. 4) représentent une classe large et variée de matériaux architecturés. Leurs caractéristiques (forte porosité, grande surface spécifique, …), les rendent très intéressantes pour de nombreuses applications dans le domaine de l'ingénierie thermique. des bâtiments (Fig. 5). Pour les applications à très haute température, les mousses céramiques réfractaires sont de plus en plus utilisées, notamment pour les fours. Les mousses à pores ouverts, présentant une grande surface spécifique, une forte porosité ouverte et une grande tortuosité, permettent d'obtenir de très hautes valeurs du coefficient de convection thermique. Produites avec une phase solide à forte conductivité, telle que le cuivre, l'aluminium ou la céramique carbure de silicium, ces mousses sont des alternatives prometteuses aux techniques conventionnelles (ailettes) (Fig. 6) [START_REF] Lu | Heat transfer in open-cell metal foams[END_REF]. Elles peuvent également être utilisées pour les bruleurs poreux (Fig. 7) [9][10] et les récepteurs solaires volumiques (Fig. 8) [11]. Fig . 6 -Echangeur de chaleur en mousse Fig. 7 -Bruleur poreux Fig. 8 -Tour solaire Jülich [11] Dans beaucoup des applications considérées, l'échange de chaleur par rayonnement peut être important: il représente une contribution de 20-30% du transfert thermique à température ambiante, il devient un mode de transfert prépondérant à haute température.

Pour modéliser correctement l'échange de chaleur dans ces milieux, il faut utiliser un approche multi-échelle [12][13] [14][15] [16], qui intègre les informations sur la structure microscopique du matériau pour modéliser le comportement macroscopique à une échelle plus grande. Pour cela, des modèles précis du rayonnement, ainsi que de la morphologie de la structure de la mousse, sont nécessaires. Des avancées significatives dans la modélisation numérique de ces deux aspects ont été obtenues ces dernières années.

Les premiers modèles morphologiques considéraient la répétition d'une seule cellule élémentaire périodique [17][2] [18][19] [20][21]. Par la suite, l'introduction de la tomographie (CT -Computed Tomography) a permis obtenir des reconstructions numériques 3D d'échantillon réelles de matériaux poreux [22][23] [24][25] [26][27] [28]. Grace aux ressources de calcul croissantes, des modelés numériques de plus en plus réalistes ont été développés, [29][30] [31], tels que ceux basés sur les partitions de Voronoi [32][33][34] [35], prenant en compte notamment le caractère aléatoire de la vraie structure des mousses.

L'étude du rayonnement a avancé à la fois numériquement et théoriquement. Les premières approches étaient basées sur la subdivision de la structure en "diffuseurs" et la sommation de la contribution des "diffuseurs" individuels [17][2] [18][19] [20][21], pour obtenir une conductivité thermique radiative équivalente basée sur le modelé de Rosseland. Les approches plus récentes, font pour la plupart appels à des simulations de Monte Carlo (MCRT -Monte Carlo Ray Tracing) de la propagation de l'énergie lumineuse dans des structures réalistes Ces dernières permettent de calculer les propriétés radiative équivalentes (coefficient d'absorption, de diffusion, fonction de phase) [36][22][23] [24][25] [26][27] [28], intervenant dans l'Equation du Transfert Radiatif (RTE -Radiative Transfer Equation). En outre, plusieurs modifications de la RTE, telles que l'Approche Multi Phase (MPA -Multi Phase Approach) [37][38] [39][40] et l'Equation du Transfert Radiatif Généralisé (GRTE -Generalized Radiative Transfer Equation) [41][42], ont été proposées pour améliorer la précision de la modélisation radiative dans ces milieux complexes et hétérogènes.

Ce travail de thèse est consacré à la fois à la modélisation de la morphologie et à la modélisation du rayonnement.

Dans la Section 1, un état de l'art sur les méthodes de détermination des propriétés radiatives est effectué, avec une attention particulière portée aux méthodes de Monte Carlo.

(A) (B) (C) Trois types de méthodes directes, inverses ou hybrides peuvent être utilisés pour calculer les propriétés radiatives équivalentes. Les méthodes directes déduisent les propriétés équivalentes directement à partir de la morphologie du milieu, les méthodes inverses sont basées sur la formulation d'un problème radiatif inverse, alors que les méthodes hybrides font appel aux deux techniques.

Concernant les méthodes de Monte Carlo, un état de l'art permet de mettre en évidence des différences notamment dans le choix de l'origine de rayons, de la technique d'intégration statistique ou encore de la condition limite aux bords. Aussi des géométries de test constituées de Sphères Opaque Identiques Superposées. Effet, une solution analytique exacte étant disponible pour cette particulière géométrie, elle nous permet d'étudier la convergence de plusieurs techniques et de déduire les choix de modélisation les plus appropriés (Fig. 10). Dans la section 2, une méthode de génération numerique de la morphologie de mousse est présentée. Trois différents types d'architecture de mousses de prosité différentes sont générés (Fig. 12),. Dans cette méthode, suite au choix d'une distribution initiale de la taille des cellules (obtenue à partir de l'analyse tomographique d'échantillons réels), on génère pour cette distribution un empilement compact de sphères [START_REF] Modest | Radiative heat transfer[END_REF], qui constitue la base du procédé de génération.

Avec cet empilement de sphères en entrée, on peut, à travers une succession de voxelisations, de filtrages d'image et finalement de remaillages utilisant iso2mesh [START_REF] Brewster | Thermal radiative transfer and properties[END_REF], obtenir des maillages de structures à pores ouverts à basse porosité.

Alternativement, l'empilement de sphères est utilisé comme donnée d'entrée pour le logiciel voro++ [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF] pour créer des diagrammes de Voronoi-Laguerre, qui sont à leur tour raffinés et stabilisés avec le logiciel Surface Evolver [START_REF] Brakke | The surface evolver[END_REF], pour obtenir enfin un maillage triangulaire des cellules. Une dernière étape est dédiée à l'ajout de l'épaisseur aux parois des cellules permettant ainsi la génération d'architecture 3D de mousses à pores fermés à haute porosité.

Alternativement, les faces des cellules sont éliminées et seul le squelette formé par les bâtons 1D subsiste. Alors, en suivant ce squelette, on ajoute des bâtons 3D, dont on contrôle paramétriquement la forme, et enfin on les soude aux intersections à travers une méthode « shrink-wrapping » [START_REF] Randrianalisoa | Radiative characterization of random fibrous media with long cylindrical fibers: Comparison of singleand multi-RTE approaches[END_REF], pour obtenir des modèles réalistes de mousses à pores ouverts à haute porosité. On compare les morphologies résultantes avec une base de données issue de l'analyse tomographique de 4 échantillons réels, en confrontant la distribution de la taille des cellules (Fig. 14) et la distribution de la connectivité des cellules (Fig. 15). Représentatif périodique (VER / REV -Representative Elementary Volume) (Fig. 17). Elle permet de calculer les grandeurs radiatives macroscopiques (tels que la transmittance, la réflectance, les facteurs de configuration, etc.).

Les résultats issus des méthodes HPA et MPA existantes présentant des écarts avec ceux issus de la méthode de référence. Des améliorations des méthodes HPA et MPA sont proposées dans le but d'accroitre la précision de calcul toute en gardant le maximum possible de simplicité. Notre approche repose sur deux idées principales:

-L'utilisation de techniques hybrides directes/inverses pour déterminer les coefficients radiatifs (Fig. 18) afin de capturer les effets de diffusion multiples.

-Le traitement des effets liés à l'histoire de la propagation de la radiation, et notamment de son origine (phase solide ou fluide). -Transmittance/Réflectance d'un échantillon (Fig. 20).

-Calcul de facteurs de forme dans une cavité cubique (Fig. 21).

-Echange de chaleur par rayonnement entre deux parois à températures différentes (Fig. 22).

Les modèles améliores se montrent régulièrement plus précis des modèles homogènes existants. (1992). The surface evolver. Experimental mathematics, 1 [START_REF] Placido | Thermal properties predictive model for insulating foams[END_REF], 141-165.

[48] Kobbelt, L., Vorsatz, J., Labsik, U. and Seidel, H.P., 1999, September. A shrink wrapping approach to remeshing polygonal surfaces. In Computer Graphics Forum (Vol. 18, No. 3, pp. 119-130).

parameters on the radiative properties (namely, the extinction coefficient) of an opaque open-cell foam. Highly accurate analytical relationships were subsequently deduced and validated by comparison with results obtained from tomography samples.

Modeling radiation in foams with a semi-transparent solid phase is substantially more complex, due to radiation propagating inside the solid. Efficient and flexible reference methods are required to validate the results obtained by homogenization models. A new, purely numerical, Direct Monte-Carlo Homogenization reference technique is proposed, that allows to simulate radiation within arbitrary cavities and calculate macroscopic radiative quantities (such as transmittance, reflectance, configuration factors etc.) based on a Representative Elementary Volume (REV) of cellular material.

The technique is validated against full scale Monte Carlo simulations and compared to the existing HPA model. The comparison reveals significant inaccuracies of conventional HPA.

On the basis of numerical and phenomenological observation of the sources of discrepancy between homogenized models and underlying discrete physical phenomena, improvements of the existing Homogeneous Phase and Multi Phase approach are proposed, their objective being to maximize accuracy while introducing as little additional complexity as possible. This is obtained through extensive use of inverse methods and the addition of one equation to take into account specific phenomena taking place in the semi-transparent solid phase. Variance of the distribution of x

ERR(x)

Error to be minimized as a function of independent variable(s) x

[x]

Iverson brackets [x] = { 1 0 𝑖𝑓 𝑥 𝑖𝑠 𝑡𝑟𝑢𝑒; 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.

General introduction

Macroporous cellular materials such as plastic (Fig. In heat insulation applications, these materials present a natural ability to simultaneously minimize natural convection, thanks to the small diameter of their pores, and thermal conduction in the solid, thanks to their low relative density to the spatial arrangement of the solid matter. At the same time they guarantee good mechanical characteristics, thanks to the interconnectedness of the solid matrix.

Foams made from plastics such as polyurethane are already used in a large number of heat insulation

applications [1][2][3][4][5][6][7]
, ranging from relatively conventional building (Fig. 0.5) and tank insulation to thermal protection systems for space vehicles. In very high temperature insulation applications, such as furnace linings, foams made from refractory ceramics such as alumina are being increasingly adopted. In many of the applications considered, thermal radiation heat transfer can have a significant influence on the heat transfer process, contributing up to 20-30% of the total heat transfer in insulation applications at room temperature and becoming the predominant mode of heat transfer at very high temperatures or in specific setups such as volumetric solar receivers.

To correctly model radiation heat transfer in these media, it is necessary to use a multi-scale approach [12][13][14] [15][16], using information about the microscopic structure of the medium to obtain useful properties to model the macroscopic behavior of the medium at a larger scale. For this to be possible, both accurate radiation models and accurate morphological models of the structure of the foam are necessary. In recent years, significant progress has been made in the numerical modeling of these media on both these accounts.

Understanding of the morphology was for a long time tied to simplistic models, usually based on regular lattices constituted by the repetition of a single periodic cell [START_REF] Placido | Thermal properties predictive model for insulating foams[END_REF] Thanks to increasing computational resources, it has been possible to create more realistic numerical models of the morphology [29][30] [31], such as those based on Voronoi cells [32][33][34] [35], which can take into account the randomness of the real foam.

Understanding of the radiation modeling has advanced both computationally and theoretically. The first approaches were based on the subdivision of the structure in multiple "scatterers" and the summing of the contribution of these "scatterers" [START_REF] Placido | Thermal properties predictive model for insulating foams[END_REF] Transfer Equation (GRTE) [41][42], have been developed to specifically cope with challenges of modeling radiation in these complex, heterogeneous media.

This work provides original contributions on both the morphological modeling front and the radiation modeling front.

In section 1, an extensive review of numerical methods for radiation heat transfer in cellular media is presented, with a special focus on Monte Carlo methods. The RTE [43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF] and general lines of multi-scale approach (Fig. 0.9) with respect to radiation heat transfer in cellular media are presented. Then, a review of existing direct numerical methods for the calculation of radiative properties of opaque cellular media is conducted, revealing significant variation in the existing methodologies, which have to be clarified to identify methods appropriate to the materials that will be studied in the following. Implications of semi-transparent behavior of the solid phase of the foam are briefly discussed, relevant models from recent literature are presented, and some limits of current knowledge are pointed out, suggesting the opportunity of devoting some effort to the development of more accurate radiative models for these materials.

The accuracy in the representation of cellular morphologies has increased considerably in recent years, evolving from early models based on simple periodic cells to realistic structures captured through tomography. Lately, methods for the direct digital generation of cellular structures have been gathering considerable interest. In section 2, we propose a framework of digital replication of foam morphologies. A single framework is proposed that allows, within a single chain of tools comprising voro++ [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF] and Surface Evolver [START_REF] Brakke | The surface evolver[END_REF], to parametrically control the digital generation of three In section 3, improvements to existing homogenized models for foams containing a semi-transparent solid phase are discussed. In these foams, modeling is more complex, due to radiation propagating inside the solid. Due to these additional difficulties, comparison of homogenized methods with reference radiative simulations have been sparse, and the few results indicate that the accuracy of conventional homogeneous models can be unsatisfactory in these materials [40][48].

Introduction

The current baseline state-of-the-art approach to model radiation in porous media such as foams is based on the multi-scale approach (Fig. 1.1) and utilization of the Radiative Transfer Equation (RTE).

General lines on the usage of RTE can be found for example in textbooks [12] An overview of radiative properties determination for porous media can be found in the monograph by Dombrovsky and Baillis [12]. Reviews dealing with radiative properties of highly porous foams can be found in [13][14][15] [16].

In the past, significant effort has been devoted to the development of inverse techniques using experimental data for hemispherical or directional transmissivity/reflectivity experiments ). The photons are simulated according to the Geometric Optics Approximation, which ensures that for wavelengths much smaller than the objects' size the photon propagate in a straight line (i.e. diffraction and other near field effects can be ignored).

Numerical methods used to predict radiative properties initially were closely related with inverse methods, making use of Monte Carlo results rather than experimental results as input data for inverse fitting [54][55]. Later, hybrid methods appeared, based on a combination of direct mean free path calculation and inverse techniques to determine the radiative properties [START_REF] Coquard | Radiative characteristics of opaque spherical particles beds: a new method of prediction[END_REF][57] [START_REF] Coquard | Radiative Characteristics of Beds of Spheres Containing an Absorbing and Scattering Medium[END_REF]. Finally, in recent years, full direct methods, which only directly use free path distribution (in the form of cumulative distribution [22][36] [37][42] [41] or mean free path [15][32] [40]) and scattering angle distribution, have been gaining favor.

In the following, we firstly present an extract of our published works, featuring an extensive review of literature on direct numerical approaches for the calculation of equivalent properties in opaque foams and an initial attempt at improvement of existing analytical correlations. Then, we discuss in a more general way direct and inverse numerical approaches to calculate radiative properties, and, finally, some challenges specifically tied to the case of cellular media with a semi-transparent solid phase.

Chapter 1 -Radiative properties modeling of open cell opaque foam: review and new analytical law

This chapter is an extract of an article published in the International Journal of Thermal Sciences [15].

It deals with direct methods to determine radiative properties of opaque foams. A review of different numerical and analytical methods is proposed, explaining the methodologies and evidencing the common points, limits and assumptions. Numerical methodologies are firstly applied to sets of spherical particles and compared with benchmark exact analytical solution. In a second step 3D

Voronoi open cell foams are generated, the various methods analyzed are tested and compared. Some attention is dedicated to the evaluation of effects due to varying degrees of irregularity in the structure and ligament. Finally, a new analytical law is proposed to determine radiative extinction coefficient

of 3D Voronoi open cell foam without significant additional computational effort. This relation is expected to be useful for preliminary optimization/design purposes.

Introduction

Cellular foams are a key material for many energy-engineering applications. Their high porosity (or low relative density) and large specific area play an important role from the thermal point of view. For example, high porosity closed cell polymer foams are used as efficient insulating materials

[1][2][3][4][5][6][7]
. Open-cell solid foams can be designed to have very low up to high values of thermal conductivity, depending on the conductivity of the solid [START_REF] Ashby | Metal foams: a design guide[END_REF][60][61] [START_REF] Wulf | Experimental and numerical determination of effective thermal conductivity of open cell FeCrAl-alloy metal foams[END_REF]. Thus, they are employed in a variety of energy related applications, such as volumetric solar energy receivers for CSP plants [11], compact heat exchangers [START_REF] Lu | Heat transfer in open-cell metal foams[END_REF], porous radiant burners [9][10] and fire barriers [63][64]. Accurate modeling of thermal properties is obviously highly desirable for the optimization of the performance in these applications. Considering the high porosity (typically in a range from 85% up to 98%), radiative heat transfer contribution can be significant, and in some cases even prevalent over other heat transfer modes. For this reason, a large number of analytical and numerical models have been dedicated to characterization of radiative heat transfer in open-cell solid foams.

Most studies focus on determining appropriate equivalent continuous medium properties.

General lines on the use of the Radiative Transfer Equation (RTE) can be found for example in textbooks [12][43] [44][45]. Various authors have developed specific adaptations of RTE for dispersed media, such as the Multi Phase Approach (MPA) [38][39] [40] and the Dependence Included Discrete Ordinates Method (DIDOM) [START_REF] Singh | Modelling radiative heat transfer in packed beds[END_REF]. However, RTE is usually considered sufficiently accurate for most practical cases, if the relevant coefficients (radiative properties) are correctly determined [39][66] [START_REF] Kamiuto | Correlated radiative transfer in packed-sphere systems[END_REF].

As such, most of the literature has been focused on finding efficient and reliable ways to determine radiative properties. An overview of radiative properties determination for porous media can be found in the monograph by Dombrovsky and Baillis [12]. Reviews dealing with radiative properties of highly porous foams can be found in [13][14].

The radiative properties may be theoretically predicted and/or identified from directional/hemispherical transmittance/reflectance measurements [17][23] [START_REF] Sacadura | Experimental characterization of thermal radiation properties of dispersed media[END_REF][68] [START_REF] Hendricks | Inverse radiative analysis to determine spectral radiative properties using Discrete Ordinates techniques[END_REF], often employing the Fourier Transform Infrared (FTIR) method.

The large number of analytical models available in literature mostly refers to independent scattering in randomly dispersed media [START_REF] Baillis | Thermal radiation properties of dispersed media: theoretical prediction and experimental characterization[END_REF]. High porosity foams are modeled as a random dispersion of particles, whose contributions are summed up to obtain the effective radiative properties. This approach, originally proposed by Glicksmann et al. [17], who modeled the foam as a set of dodecahedral cells, was also followed by Placido et al. [START_REF] Placido | Thermal properties predictive model for insulating foams[END_REF] for polymer closed foam, and Baillis et al.

[18] [19] for open cell carbon foams. Coquard et al. [20] and Loretz et al. [21] extended these results by considering models with different cells and strut shapes.

Independent scattering approaches remain largely prevalent in literature, probably thanks to their comparative simplicity. A typical limitation of these studies is the difficulty to account for shadowing effects. Independent scattering approaches also typically require the knowledge of a number of geometrical parameters of the foam, such as the strut diameter or cell size, which are difficult to determine univocally, because of the intrinsically random quality of the real foam structures.

To overcome these limitations, alternative numerical methods have been developed and, thanks to increasing computational power, have lately gained large popularity. In particular, numerical methods based on Monte-Carlo techniques for the determination of radiative properties are becoming very popular in order to study either real structures obtained from tomographic imaging or computer generated structures that closely mimic the micro structure of the real foams.

Tancrez and Taine [36] proposed to use the Radiative Distribution Function Identification (RDFI) model and determined radiative properties of spherical packed beds. Zeghondy et al. [22][23] and Petrasch et al. [24] applied the RDFI approach to tomographic data. Coquard et al. [25][26] [27] proposed to use an alternative Monte Carlo approach based on mean free path calculation.

In addition to Monte Carlo methods, some alternative numerical methods have been presented in literature. Most notably, Loretz et al. [28] presented a geometric approach to rapidly calculate extinction coefficient from open cell foam tomographic data.

Techniques based on tomographic data provide satisfactory agreement with experimental data, but their dependence on high quality scans of existing foam samples makes them of limited utility for design purposes. Some recent studies have sought to overcome these limitations by digitally reproducing the foam structures using different approaches, including mathematical morphology operations applied on existing tomography data [29][30], simulation of the bubbling process [31],

regular [32][33] and irregular [34][35] Voronoi partitions. By computer generating a number of structures and running numerical simulations [29][30] [35] it is possible to obtain useful results for the optimization of energy transfer. Irregular 3D Voronoi structures seem to be particularly promising for this purpose as they replicate the structures of real foam pretty well [START_REF] Randrianalisoa | Microscale direct calculation of solid phase conductivity of voronoi s foams[END_REF].

In the current paper, a special emphasis is put on the predictive models of radiative properties of cellular foams with open cells and high porosities (ε > 85%). In the light of the state of the art, a representative selection of numerical and analytical methods is presented in a comparative fashion, focusing on their similarities and differences, their strong points, limits and assumptions.

Subsequently, the methods are numerically compared. First, sets of spherical geometries are generated. For such geometries, an exact analytical solution exists [10] that can be used as a benchmark to evaluate the numerical methodologies. 3D Voronoi open cell structures are also considered, to provide a more realistic representation of the foam. Structures with two different degrees of irregularity and two different ligament shapes (circular and triangular) are considered. The various methods analyzed are tested and compared. Based on the results obtained, guidelines are proposed to allow optimal choice of the numerical method. Additionally, some corrections to commonly used analytical relations are proposed, that should improve their accuracy without significant additional computational effort. This new relation should provide useful guidance for preliminary optimization/design tasks.

Radiative properties modeling

Porous media such as foams are usually considered as equivalent continuum media and the radiative transfer equation (RTE) can be used [43][45]: is much greater than one, then radiation can be approximately treated as straight rays propagating through the medium, i.e. Geometric Optics Approximation prevails and diffraction effects can be neglected. This is the case for strut diameter in the 100μm range, as the peak radiation intensity at room temperature is at a wavelength of around 10μm, and the length of a ligament is typically much larger than its diameter. In addition, the struts are organized in quasi-polyhedral cells with quasi-polygonal windows comprising, so that the typical spacing between them is comparable to their length. Thus near field effects can be neglected. On the other hand, shadowing effects cannot be neglected.

Ω • ∇I = 𝛼𝐵 0 -𝛽𝐼 + 𝜎 4𝜋 ∫ I(Ω′)Φ(Ω ′ , Ω)𝑑Ω′ 4𝜋 (
In the present paper, solid phase is considered opaque (as in the case of metal and carbon foam). The reflection is assumed to be diffuse and, additionally, reflectivity is considered independent of incidence angle, as it is assumed in several previous works [22][24][25] [26][36] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF]. Under this assumption, the following relation holds for the determination of the scattering albedo:

𝜔 = 𝜌 𝑠 (1.
2) Additionally, while Monte Carlo methods allow to numerically calculate the scattering phase function, it has been shown [24] that under the assumption of diffuse reflection, the numerically calculated scattering phase function for open cell foam structure closely matches the scattering phase function for opaque large spheres or randomly oriented convex opaque particles [43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF]:

𝛷(𝜃) = 8 3𝜋 (sin 𝜃 -𝜃 cos 𝜃) (1.
3)

The extinction coefficient, β, is more difficult to determine: different methodologies to calculate it will be the focus of the following discussion.

Analytical method

As recalled in the introduction, a large number of analytical relationships have been presented to calculate the extinction coefficient directly from the knowledge of the geometrical characteristics of the structure. The foam is assumed as a random dispersion of particles under the hypothesis of independent scattering. These relationships are based on a common framework, which can be summarized as follows:

-Create a representation of the structure as a collection of randomly dispersed simple scatterers (usually simple particles such as cylinders to model the struts and cubes or spheres to model the junctions).

-Evaluate the number of particles/scatterers per unit volume Nv on the basis of a representation of the overall structure in terms of regular polyhedral cells composed by struts and, sometimes, junctions.

-Evaluate the average projected area G of a single scatterer.

-At this point, in the hypothesis of independent scattering one can sum up the contributions of the single scatterers.

Assuming geometric optics prevails and diffraction effects can be neglected, the extinction coefficient is obtained as 𝛽 = G • Nv [START_REF] Brewster | Thermal radiative transfer and properties[END_REF].

In fact, for any convex (no negative curvature) randomly oriented scatterer one can obtain the average projected area as 𝐺 = 𝑆 4 ⁄ [START_REF] Hottel | Radiative transfer[END_REF], where S is the particle surface area. It follows:

𝛽 = 𝐺 • 𝑁 𝑣 = 𝑆 4 ⁄ • 𝑁 𝑣 = 𝑆 𝑣 4 ⁄ (1.4)
Where Sv is the specific surface area i.e. total surface area per unit volume.

It is worth noting that some works have studied the limits of independent scattering and proposed scaling correlations in order to determine the extinction coefficient of spherical beds [START_REF] Singh | Modelling radiative heat transfer in packed beds[END_REF][66] [START_REF] Kamiuto | Correlated radiative transfer in packed-sphere systems[END_REF] or fibrous media [START_REF] Coquard | Radiative properties of dense fibrous medium containing fibers in the geometric limit[END_REF] for relatively large solid volume fractions. In particular, Brewster [START_REF] Brewster | Volume scattering of radiation in packed beds of large, opaque spheres[END_REF] presented a review of pre-existing correlations together with a simple correlation that allows correcting for "volume scattering effects" in spherical packed beds:

𝛽 = 𝑆 𝑣 4 ⁄ • 𝑆 𝑟 = 𝑆 𝑣 4 ⁄ • 1 𝜀 ⁄ = 𝑆 𝑣,𝑓 4 ⁄ (1.5)
Where Sv,f is the specific surface area per unit volume fluid. This relation will be used as analytical reference in the following as it was shown in [36] to be rigorously valid for geometries constituted of randomly distributed Overlapping Opaque Spheres. A similar corrected relation does not exist in the case of cellular materials such as foams: it is among the objectives of this study to propose such a relation.

Scattering albedo and phase function are generally considered to be left unchanged by volume scattering and shadowing effects, so that relations (1.2) and (1.3) can be used to determine them.

Projection method

Loretz et al. [28] proposed a simple and fast approach to calculate the extinction coefficient.

Exploiting the typical digital form of open cell foam tomographic data (collections of 2D slices), this method is based on incremental "projection" of subsequent slices and direct evaluation of optical obstruction by pixel counting. Fig. 1.1 shows a sketch of the projection process taken from the original work [28].

Fig. 1.1 -Projection of tomography slices in [28] As the projection proceeds, the surface fraction of white pixels 𝑊 𝑛 remaining after Ns slices of thickness ∆𝑡 can be equated to the fraction of radiation transmitted, so that according to the Beer Law:

𝛽(𝑠) = - Log(𝑇𝑟𝑎𝑛𝑠𝑚𝑖𝑡𝑡𝑒𝑑 𝑅𝑎𝑑𝑖𝑎𝑡𝑖𝑜𝑛) 𝑇ℎ𝑖𝑐𝑘𝑛𝑒𝑠𝑠 = - Log(𝑊 𝑛 ) N 𝑠 × ∆𝑡 (1.6)
Low computational effort and direct use of the data structure are the main advantages of this method.

It also makes it possible to readily evaluate anisotropy in the radiative behavior.

The method has been originally conceived for volumetric data structures obtained from X-Ray tomography. To treat the polygonal mesh geometries used in the current study, a preparatory step involving the conversion of the meshes to a voxel structure has been required. Sensitivity to the resolution of the voxelization step has been reported in the numerical results section. Additionally, isotropic geometries have been considered, so the final value of 𝛽 presented for the algorithm is taken in all cases as the average of the x, y and z directions. Dispersion is quite small in any case: typical values of average absolute deviation and standard deviation are within 2.5% of the average.

Monte Carlo method

Ray casting methods to determine radiative coefficients use a Representative Elementary Volume (REV) to extract relevant properties, and the rays casted are considered as representative of typical emission/absorption/scattering events happening in the structure.

Different mathematical methods can be used to treat the resulting data. Other methodological variations observed in literature include the distribution of ray origins and the handling of the intersection of the rays with the bounding box of the geometries (REVs) undergoing numerical analysis. These issues will be discussed in detail in the following.

Variations on ray casting method -Ray origin

The choice of origin for the "representative" ray distribution is not obvious, and in fact, different assumptions have been used in literature. Two intuitive choices can be presented (Fig. 1.2):

(A) Casting rays starting from random points inside the fluid/void phase [22][24] [36] (B) Casting rays starting from the solid/fluid interface [26][27] These two choices can be linked to two somewhat distinct radiation regimes:

(A) To a regime of external radiation flux imposed over a substantially passive (i.e. absorbing and scattering, but not emitting) foam.

(B) To a regime of radiatively active (high temperature) foam, i.e. behaving at the same time as emitting and absorbing/scattering fluid interface.

Additionally, while for (A) the obvious choice of angular distribution of rays is the uniform distribution, for case (B) various different distributions can be proposed (uniform, lambertian, etc.).

In the following, a uniform distribution has been utilized, as proposed in [26]. The effects of such modeling choices will be analyzed hereafter.

Variations on ray casting method -Boundary handling

Another potentially relevant issue that has been rarely explicitly addressed in papers is the handling of rays that exit the bounding box of the volume under consideration (RVE) without intersecting the solid phase.

Some authors [24] purposefully choose the location of origin of the rays to minimize this fraction of rays and do not consider these rays in subsequent calculations. In some cases, when the total optical thickness of the data set under consideration is low, such an outcome could be difficult to obtain leading to noticeable errors in the estimation of radiative properties. It must also be considered that the size of the ray origination zone should be large enough to be itself a RVE with respect to the entire structure. For these reasons, in a number of cases it is highly desirable to be able to originate the rays from any point in the structure, and thus a methodology to handle the rays exiting the bounding box without intersecting the solid phase is required.

In the current study we have chosen to compare three examples of possible methodologies. The three methodologies are schematically presented in Fig. 1.3 and can be summarized as follows:

(A) Symmetrically reflecting the ray as it hits the boundaries. [32] (B) Wrapping around the ray in a periodic fashion, i.e. restarting the ray with the same direction from the boundary opposed to the one that has been intersected. To ensure preservation of solid/fluid continuity and improve convergence, this method also requires the starting location on the opposite side to be randomized, so that the ray starts on the opposite face on a random point of the void domain. [25] (C) Randomly casting a new ray with a random direction from a random point inside the void domain, treating this completely newly casted ray as a continuation of the previously casted ray.

[32] While method (A) preserves both local and direction continuity, method (B) preserves direction and some location information (the ray that exits a boundary newly appears at a boundary), while (C) sacrifices all information regarding both direction and position, preserving only local phase coherency, thus essentially eliminating ordering phenomena due to the boundary handling.

Calculation of radiative properties

The calculation of the extinction coefficient does involve some mathematical complexity.

Two distinct approaches based on Monte Carlo method can be distinguished in literature for the calculation of the extinction coefficient. In the RDFI method [36] and other similar methodologies (discussed in the introduction), an extinction distribution function is defined as:

𝐺 𝑒 (𝑠) = 1 4𝜋𝑉 𝑓 ∫ ∫ ∫ 𝛿[𝑠 ′ -𝑠 0 (𝑟, 𝑢)]𝑑𝛺(𝑢)𝑑𝑟𝑑𝑠′ 4𝜋 𝑉 𝑓 𝑠 0 ≅ 1 -𝑒 -𝛽 𝐼 𝑠 (1.7)
Which represents the fraction of rays intercepted at lengths lesser or equal then s. 𝛽 𝐼 is found by fitting an exponential distribution to the actual ray distribution using a least squares method.

On the other hand, in the method presented [25] and other similar methodologies, the value of the extinction coefficient is found as the inverse of the average beam length.

𝛽 𝐼𝐼 = 1 ∑ 𝑠 𝑛 𝑁 𝑛 𝑁 = 1 𝑠 𝑎𝑣𝑔 (1.8)
It must be noted that, while these two formulations are mathematically different, each one corresponds to a typical physical description of the scattering/absorption phenomenon in semitransparent media. Namely, the formulation given in [36] corresponds to a description of 𝛽 as ray interception probability, while the formulation given in [25] corresponds to a description of 𝛽 as the inverse of photon mean free path.

In fact, it can be shown that, for media that strictly obey a Beer-Lambert distribution of ray population, the two formulations are equivalent. For any given extinction distribution Ge(s), we can write:

𝛽 𝐼 𝑖𝑠 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡: ∫ [𝐺 𝑒 (𝑠) -(1 -𝑒 -𝛽 𝐼 𝑠 ) ] 2 𝑑𝑠 ∞ 0 = min 𝛽 𝐼 ′ {∫ [𝐺 𝑒 (𝑠) -(1 -𝑒 -𝛽 𝐼 ′ 𝑠 ) ] 2 𝑑𝑠 ∞ 0 } (1.9) 𝛽 𝐼𝐼 = 1 ∫ 𝑠 𝑑𝐺 𝑒 (𝑠) 𝑑𝑠 𝑑𝑠 ∞ 0 (1.10)
Assuming a Beer -Lambert distribution function, i.e. 𝐺 𝑒 (𝑠) = 1 -e -𝛽𝑠 , we obtain:

𝛽 𝐼 𝑖𝑠 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡: ∫ [𝑒 -𝛽𝑠 -𝑒 -𝛽 𝐼 𝑠 ] 2 𝑑𝑠 ∞ 0 = min 𝛽 𝐼 ′ {∫ [𝑒 -𝛽𝑠 -𝑒 -𝛽 𝐼 ′ 𝑠 ] 2 𝑑𝑠 ∞ 0 } ; 𝛽 𝐼 = 𝛽 (1.11) 𝛽 𝐼𝐼 = 1 ∫ 𝑠𝛽e -𝛽𝑠 ∞ 0 𝑑𝑠 = 1 1 𝛽 = 𝛽 (1.12)
Therefore, in general, the more closely the real ray distribution approximates a real exponential distribution (i.e. the more accurately the radiative behavior of the medium can be predicted using an HPA approach), the better the agreement between the two methods.

In the light of this finding, we chose to use the inverse mean free path calculation as it affords some advantages. Notably, not requiring a fitting step, it allows to directly and univocally obtain the value of beta, thus resulting simpler to implement and more robust. This is particularly true since the least squares minimization to be performed is non-linear.

Numerical results and comments

As announced in the introduction, a set of geometries constituted of randomly distributed identical overlapping opaque spheres (IOOS) has been generated, because of its simplicity and because the analytical solution is known and could be used as benchmark.

In a second step, to replicate as closely as possible real foam structures, cellular structures based on Voronoi partitions of space have been considered as in recent works [34][35].

Spheres

Geometry description

Six subsets of spherical geometries have been generated. Each set is constituted of 2000 spheres of identical radius, whose centers are randomly distributed in a cube of 1000x1000x1000 units 3 . An example of the resulting geometries is presented in Fig. 1.4. The spheres can overlap (Fig. 1.4a).

Additionally, all the spheres intersecting the 1000x1000x1000 bounding box have been "cut" to the bounding box (Fig. 1.4b). Each set has a different level of porosity and specific surface, determined by the radius of the spheres. Table 1.1 summarizes the relevant geometric data. 1.1 -Summary of the geometrical characteristics of the random spheres geometries.

Results and comments

Computed values of the extinction coefficient are presented in Fig. 1.5 as a function of porosity. The analytical results are calculated making use of relation (1.5). The Monte Carlo calculations that led to the results shown in this figure were carried out considering 10 6 rays and periodic wrapping at the boundaries of the box, with rays starting from the solid surface. Influence of variations in Monte Carlo methods will be discussed in the following. The expected trend of the extinction coefficient increasing with decreasing porosity is found. It can be seen that, for this set of benchmark geometries, all the proposed methodologies converge to very 0.0E+00 similar values. In fact, the values obtained using numerical methodologies all lie within a ± 1.5% interval centered on the analytical benchmark.

As can be readily noted, the results diverge (though slightly) as porosity decreases. This can be most probably attributed to a mix of meshing effects and boundary cutting effects that positively correlate with increasing sphere diameter.

Thanks to the availability of an analytical solution, the spheres were utilized for a series of convergence tests of the methodologies proposed. For the tests, an "average error" quantity was evaluated as follows:

𝐴𝑉𝐸𝑅𝐴𝐺𝐸 𝐸𝑅𝑅𝑂𝑅 = ∑ | 𝛽 𝑖,𝑛𝑢𝑚 𝛽 𝑖,𝑎𝑛 | 6 𝑖=1 6 (1.13)
Where indices 1 to 6 correspond to the six distinct spherical geometries.

First, an assessment of numerical convergence was tried for the different methods. For Monte Carlo methods numerical convergence depends on the number of rays cast, while for the projection methods it depends on the resolution (in total voxels) of the grid. The results of this test are depicted in Fig.

1.6.

Total voxels

Number of rays cast

Average error

Monte Carlo Projection method

Various Monte Carlo methodologies showed similar results, so they have been represented as a band in the chart. It can be seen that Monte Carlo methods stabilize between 10 5 and 10 6 rays cast, while projection method stabilizes around 10 7 and 10 9 total voxels. It must be however noted that the projection method has a higher average error. This may be due to sample size requirements, as will be discussed later.

Focusing specifically on Monte Carlo methods, we performed an assessment of the effects of different boundary conditions and different sample sizes on the results. To obtain this, in addition to the base geometries, other geometries constituted of a smaller number of spheres (1000, 500 and 250 spheres)

were generated and then processed. Monte Carlo surface methods using three different boundary conditions, as proposed in the relevant section, were considered. As can be seen in Fig. 1.7, while for all the considered boundary handling methods the average error is below 1% for the largest sample size, methods employing some degree of randomization perform better than symmetrical mirroring, and more specifically the periodic wrap around with position randomization method achieves the best results for all cases considered. For this reason, this was the boundary condition used in the rest of the study.

To show the interest of being able to manage rays exiting the domain, calculations were also carried out with rays shoot from the core of the structure, as proposed by various authors [22][24], using the largest geometries (2000 spheres in a 1000 3 cube). A region of 100 3 units, centered in the geometry, was chosen to cast the rays. Results are reported in Fig. 1.8 and compared with those of the periodic randomized wrap around method. It must be noted that, due to the size and optical thickness of the geometries generated, it was not possible to guarantee the condition of a very small portion of rays exiting the domain, typically required by authors using this kind of treatment of boundaries. As a result, the errors are very large, particularly for high porosity structures. This in turn means that opportune choice of the boundary Additionally, Monte Carlo methods with rays starting from the void or from the surface have been tested and compared with each other, to check if significant differences exist in convergence speed with medium size. The results are summarized in Fig. 1.9. In this case too, results are similar and overall satisfying, but one of the two strategies, surface casting, appears to slightly outperform the other.

The "best" Monte Carlo strategy (coupling surface casting with periodic randomized wrap around) was then compared to the projection method in terms of convergence with medium size. The results are compared in Fig. This finding can be explained as follows: projection algorithm starts from the boundary surfaces and only uses the three coordinate directions for projection, while Monte Carlo methods casts rays in a number of directions from many points within the domain. It is to be expected that the former would be faster but extract less information than the latter.

Voronoi structures 1.3.2.1 Geometry description

The porous materials analyzed in the previous sections are suitable for model validation. In an effort to mimic real solid foams, cellular structures generated on the basis of Voronoi partitions of space have also been considered. The final goal is to evaluate the suitability of the simple analytical approaches currently used for the evaluation of the radiative properties of solid open cell foams, evaluate their accuracy and range of validity and to propose new, more precise laws than the ones currently used.

Average error

Number of scatterers

Projection method Monte Carlo

It has been shown in prior studies that the Voronoi tessellation technique allows a faithful representation of the porous morphology not only for polymeric foams closed cell foams [34] but also for solid open-cell foams [START_REF] Randrianalisoa | Microscale direct calculation of solid phase conductivity of voronoi s foams[END_REF]. As previously described in [34][71], the cells are created from a distribution of nucleation points (centres). A so called Voronoi paving associates a region of the surrounding space to each nucleation point. The cell associated to a given particle corresponds to all of the space which is closer to that particle than any other. By positioning the nucleation points along a regular centred cubic lattice, the uniform polyhedrons generated are tetrakaedecahedron (the Kelvin cells) while a face centred cubic lattice leads to rhomboidal dodecahedron cells. It is also possible to generate volumes composed of cells with non-uniform shape and size distributions. The method used in the present work starts from a regular distribution of points and consists in allowing a random small displacement to each nucleus in a restricted area around their initial positions. The size of this area, normalised by the distance between two initial nuclei, defines a perturbation coefficient . By varying the value of , one can obtain cellular structures with a non-uniform distribution of the cell's size and of the cell's shape. Large values of  leads to wide distributions while for  = 0, one obtains the regular structures mentioned previously (tetrakaedecahedron or rhomboidal dodecahedron).

In the present study, we used the free, open source Voro++ software [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF]

Results and comments

In the light of the results obtained for the spherical geometries, additional convergence calculations have been carried out to test numerical models for the Voronoi geometries. The testing has been specifically focused on evaluating the convergence for the different methods (Projection method and Monte Carlo method with rays cast from surface and periodic randomized boundary wrapping).

Initially, sub-cuts of the base geometries have been generated, with volumes 2, 4 and 8 times smaller than the original. Then results obtained for these reduced cuts have been compared to the results obtained for the largest geometry, and the evolution of average relative deviations with increasing sample has been evaluated. Additionally, to exclude geometry sampling effects, results have been normalized using the respective Sv ratios. In Fig. 1.12, the average relative deviations are plotted as a function of relative sample volume i.e. the ratio of the sub-cut volume to the original volume.

Fig. 1.12 -Convergence of numerical methods as a function of geometry sub sampling.

As can be seen, while the Monte Carlo methods both exhibit good convergence properties, the projection method is still non converging at the geometry sizes considered. It may also be observed that results are fully consistent with those recovered for the spherical geometries: the Monte Carlo method performs significantly better than the projection method. For this reason the results obtained by projection method have been excluded from further analysis.

In Fig. 1.13 the computed values of the extinction coefficient are presented as a function of porosity for different cell size distribution dispersions () and strut shapes, along with values obtained using analytical relation (4), commonly used in literature for cellular materials.

Interestingly, it appears that cell size distribution dispersion has no strong effects on predicted extinction coefficients. Structures with triangular struts show significantly higher extinction coefficients than structures with circular struts for any given level of porosity. Higher specific surface for a given porosity seems to be the most likely justification of this behavior. The average value of the numerical results is persistently higher than the analytical result, in range from 7% to 30%, averaging 19% across the porosity range.

Average relative deviation

Relative sample volume

Conclusions

In the present paper, various numerical methodologies have been compared for different porous media. The geometries included random overlapping opaque sphere arrangements and Voronoi structures, which make possible to generate realistic foam structures.

Specific attention has been devoted to Monte Carlo methods. Different Monte Carlo approaches

proposed in literature to determine β (namely, RDFI method based least squares fitting identification and inverse mean free path method) have been shown to have a common basis. This prompts to favour the inverse mean free path methodology because of its simplicity. suggest that Monte-Carlo methods converge faster with smaller sized domains, so that Monte-Carlo methods are recommended when there is uncertainty about the sufficient size of the sample, while the projection methods could be used when there is need for faster computation.

Concerning Voronoi structures, it has been found that analytical relations such as those used in previous literature give results that are consistently lower than the numerical methodologies.

Chapter 2 -Indirect identification methods and semitransparent cellular media 2.1 Classification of methods for the identification of radiative properties

While direct methods for radiative properties identification have become prevalent in most recent literature thanks to their simplicity, some features of the inverse methods which were historically developed first make them still interesting. In this section we classify the methods presented in literature as direct, inverse or hybrid and present their general schemas and some of their strengths and weaknesses. In Figs. 2.1-2.3, the schemes of the three are illustrated. The methods are presented as illustrated flow diagrams, for the sake of readability. Finally, Fig. 2.3 gives a general scheme of hybrid direct-inverse methods, such as those seen in [START_REF] Coquard | Radiative characteristics of opaque spherical particles beds: a new method of prediction[END_REF][57] [START_REF] Coquard | Radiative Characteristics of Beds of Spheres Containing an Absorbing and Scattering Medium[END_REF] and the one that will be presented in the following. Direct Monte Carlo simulations in the complex medium are used to collect statistical information about the medium (both ray statistics such as free path distribution and physical quantities such as transmittance/reflectance). A part of this information is processed directly and allows to establish definite values for a part of the radiative properties. Another part of the data is fed as input to an inverse method, to determine the rest of the parameters. When judiciously applied, such methods allow to alleviate most parameter sensitivity problems typical of pure inverse methods, while preserving some advantages such as the capability to take into account multiple scattering and the direct comparison of real and equivalent medium.

Modeling radiation in cellular media with a semi-transparent solid phase.

When the solid material of the cellular medium itself is semi-transparent to radiation, some significant additional modeling challenges appear. In this case, the radiation propagates in the solid as well as in the fluid part of the structure. For semi-transparent, non-magnetic media, assuming smooth surfaces, the scattering of light at the interface is governed by Snell's law and Fresnel equations [43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF].

(2.1) Where n1 and n2 are the refraction indices of the two media, θi and θt are the incidence and transmitted angle (Fig. 2.4) and R is the reflectivity. An important thing to note with respect to these equations is that Snell's law has no possible solutions for θi > θc = sin -1 𝑛 1 𝑛 2

(2.2)
. In case of passage from a more dense to a less dense medium (n1 > n2), all rays incident above the critical angle θc will be totally internally reflected.

In general, in semi-transparent media, the propagation behavior in the two phases can be expected to be quite different, and inter-phase radiation transfer is expected to be significant. It must be considered that the underlying model of the Homogeneous Phase Approach (HPA) is that of disperse, non-interacting scatterers distributed uniformly in the domain. In foams with a semi-transparent solid phase, it is clear that these conditions are not achieved.

Notes on the Multi Phase Approach (MPA)

One recent development partially addressing these concerns is the so called Multi Phase Approach (MPA). This approach has been presented independently in Gusarov [38] and Lipinski et al. [37]. In this approach the medium is modeled as two coupled homogeneous phases, each one with its own Radiative Transfer Equation. The MPA is very useful in the case of Local Thermal Non Equilibrium (LTNE) between the two phases, but it is also generally interesting in cases where the solid phase is semi-transparent and its volume fraction is not negligible.

It can be useful recalling the basic equations and setup. In a standard MPA, such as presented in Gusarov [38], the solid phase and the fluid phase are homogenized separately, then interfacial interactions (and thus the coupling of the phases) are handled through four interfacial scattering coefficients σ00, σ01, σ10, σ11, and their respective phase functions Φ 00, Φ 01, Φ 10, Φ11. The intrinsic properties of the phases (e.g. absorption coefficient) are not homogenized, and each phase preserves its actual intrinsic properties.

Ω • ∇I 1 = -𝛽 1 𝐼 1 + 𝜎 11 4𝜋 ∫ 𝐼 1 (Ω′)Φ 11 (Ω ′ , Ω)𝑑Ω′ 4𝜋 + 𝜎 21 4𝜋 ∫ 𝐼 2 (Ω′)Φ 21 (Ω ′ , Ω)𝑑Ω′ 4𝜋 (2.3a) Ω • ∇I 2 = -𝛽 2 𝐼 2 + 𝜎 22 4𝜋 ∫ 𝐼 2 (Ω′)Φ 22 (Ω ′ , Ω)𝑑Ω′ 4𝜋 + 𝑓 1 𝑓 2 𝜎 12 4𝜋 ∫ 𝐼 1 (Ω′)Φ 12 (Ω ′ , Ω)𝑑Ω′ 4𝜋 (2.3b)
It is useful to add an emission contribution and change the algebraic form of the third term (the interphase coupling terms) we get:

Ω • ∇I 1 = 𝑛 1 2 𝛼 1 𝐵 1 0 -𝛽 1 𝐼 1 + 𝜎 11 4𝜋 ∫ 𝐼 1 (Ω′)Φ 11 (Ω ′ , Ω)𝑑Ω′ 4𝜋 + 𝑓 2 𝑓 1 𝜎 21 4𝜋 ∫ 𝐼 2 (Ω′)Φ 21 (Ω ′ , Ω)𝑑Ω′ 4𝜋 (2.4a) Ω • ∇I 2 = 𝑛 2 2 𝛼 2 𝐵 2 0 -𝛽 2 𝐼 2 + 𝜎 22 4𝜋 ∫ 𝐼 2 (Ω′)Φ 22 (Ω ′ , Ω)𝑑Ω′ 4𝜋 + 𝑓 1 𝑓 2 𝜎 12 4𝜋 ∫ 𝐼 1 (Ω′)Φ 12 (Ω ′ , Ω)𝑑Ω′ 4𝜋 (2.4b)
Where f are volume fractions and B 0 are blackbody intensities. We can write:

𝛽 1 = 𝛼 1 + 𝜎 11 + 𝜎 12
(2.5a)

𝛽 2 = 𝛼 2 + 𝜎 22 + 𝜎 21 (2.5b)
Where σ11 and σ22 can include the contribution of the intrinsic diffusion of phases 1 and 2 respectively.

𝜎 11 = 𝜎 11,𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 + 𝜎 1,𝑖𝑛𝑡𝑟𝑖𝑛𝑠𝑖𝑐 (2.6a
)

𝜎 22 = 𝜎 22,𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 + 𝜎 2,𝑖𝑛𝑡𝑟𝑖𝑛𝑠𝑖𝑐 (2.6b)
For a general overview of the physical phenomena involved and the corresponding coefficients, it can be useful to refer to Fig. 2.5. It's easy to verify that the chosen values verify the conservation of energy in the case of non-absorbing media (α1 = α2 = 0). Following the integration seen in [38], we get:

∇ • ∫ Ω(𝑓 1 𝐼 1 + 𝑓 2 𝐼 2 )𝑑Ω 4𝜋 = ∫ [(-𝛽 1 + 𝜎 11 + 𝜎 12 )𝑓 1 𝐼 1 + (-𝛽 2 + 𝜎 22 + 𝜎 21 )𝑓 2 𝐼 2 ] 4𝜋 𝑑Ω = 0 (2.7)
The usefulness of the algebraic substitution for the coupling term is seen here, as one can immediately see that the energy conservation is automatically satisfied. It's worthwhile observing that the model presented is equivalent to that illustrated in Lipinski et al. [37].

It is possible to obtain an additional constraint on σ12 and σ21, by requiring the conservation of energy in local thermodynamic equilibrium. In this condition we can posit:

𝑇 1 = 𝑇 2 = 𝑇 (2.8) 𝐵 1 0 = 𝐵 2 0 = 𝐵 0 (2.9) 4𝜋𝑛 1 2 𝛼 1 𝐵 0 = 𝛼 1 ∫ 𝐼 1 (Ω′)𝑑Ω′ 4𝜋 (2.10a) 4𝜋𝑛 2 2 𝛼 2 𝐵 0 = 𝛼 1 ∫ 𝐼 2 (Ω′)𝑑Ω′ 4𝜋 (2.10b)
In addition, each phase must be in radiative equilibrium:

𝑓 1 ∇ • ∫ Ω𝐼 1 𝑑Ω 4𝜋 = 4𝑓 1 𝜋𝑛 1 2 𝛼 1 𝐵 0 -(𝛽 1 -𝜎 11 )𝑓 1 ∫ 𝐼 1 4𝜋 𝑑Ω + 𝜎 21 𝑓 2 ∫ 𝐼 2 4𝜋 𝑑Ω = 0 (2.11a) 𝑓 2 ∇ • ∫ Ω𝐼 2 𝑑Ω 4𝜋 = 4𝑓 2 𝜋𝑛 2 2 𝛼 2 𝐵 0 -(𝛽 2 -𝜎 22 )𝑓 2 ∫ 𝐼 2 4𝜋 𝑑Ω + 𝜎 12 𝑓 1 ∫ 𝐼 1 4𝜋 𝑑Ω = 0 (2.11b)
From Eq. (2.11a), or identically Eq. (2.11), with Eq. (2.5a-b), (2.8), (2.9) and (2.10a-b), after simplification we obtain:

𝜎 10 𝑓 1 ∫ 𝐼 1 4𝜋 𝑑Ω = 𝜎 01 𝑓 0 ∫ 𝐼 0 4𝜋 𝑑Ω (2.12)
We then apply Eq. (2.10a-b) again and we simplify to obtain:

𝑓 1 𝑛 1 2 𝜎 12 = 𝑓 2 𝑛 2 2 𝜎 21 (2.13)
While this equation is trivially satisfied if one obtains the values of β and σ using the methods proposed in [38][74], if these values are obtained otherwise (e.g. with a Monte Carlo method), this may not be the case. It will reveal itself useful later in this work, when developing an improvement of the existing MPA model. It may be noted that Eq. (2.13) is also implied in Dauvois et al. [42], who follow a more complex derivation that is valid for non-isotropic media where 𝜎 𝑎𝑏 = 𝜎 𝑎𝑏 (Ω).

Other recent developments -Generalized Radiative Transfer Equation

Another significant issue of conventional RTE description of propagation of radiation in cellular media is the assumption of Beerian (i.e. exponential) extinction of radiation. Especially in the case of the solid phase, this assumption can turn out to be far from reality. A novel approach, called Generalized Radiative Transfer Equation (GRTE) [41] has been recently proposed to address this problem. In the GRTE approach, rather than using a simple exponential decay function (characterized by the extinction coefficient β) and a scattering albedo σ, the entire extinction and scattering distribution functions are plugged into the RTE:

(2.14)

In another very recent development, Dauvois et al. [42], building upon a GRTE model, have also addressed the possible effects of dependence of a ray's propagation behavior from its history. They achieved this by classifying their extinction length information according to ray history in microscale MCRT simulation, and then using the full hierarchy of extinction length distributions (organized on three tiers) in the homogenized simulation, according to ray history:

(2.15a)

(2.15b)
According to their formulation, a ray is assigned a different extinction/scattering distribution function depending on its previous history, with up to 3 phase transitions taken into account. They've also

shown that the results of their more accurate model differ quite significantly from those of a standard MPA [42].

General observations

Unfortunately, literature comparing these advanced homogenized models with reference solutions or full Monte Carlo simulations is scarce and the evidence is mixed [40][48]. Another limit of these models is that they introduce significant complication in the resolution of the homogeneous equations.

The MPA doubles the number of equations and requires the addition of coupling terms to the problem.

In the case of the GRTE, on top of the substantial complication of the equations themselves (especially when ray histories are accounted for), their form complicates substantially the adoption of conventional methods of resolution of the homogenized RTE, such as the Discrete Ordinates (DO) method, requiring Monte Carlo solution of the radiation field in the homogenized medium. In all cases this has significant impacts on computational costs, which are significantly higher than those of the HPA.

In the light of this state of affairs, we concentrated our efforts in the field of radiative modeling in semi-transparent foams on two main paths:

1) The development of more efficient and general reference methods (Section 3).

2) The development of new, more accurate, yet simple, homogenized methods (Section 4).

Chapter 3 -A general framework for the generation of realistic cellular morphologies. In the present work, a methodology is proposed that makes it possible to generate realistic Representative Volume Elements (RVE) of porous foam structures. This methodology, which will be explained in detail in the following, involves the generation of a pseudo-random periodic packing of seeding spheres, the generation of a Voronoi-Laguerre diagram of these spheres using voro++ [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF], the stabilization of the resulting cell structure using Surface Evolver [START_REF] Brakke | The surface evolver[END_REF], and the extraction of the skeleton from the stabilized cell structure. Three different types of porous morphologies can be obtained, covering a large variety of porous foam materials:

1. The packed spheres, after a process of further inflation, conversion to voxels and porosity fine-tuning and inversion, can finally be used as basis for meshing open cell porous structure such as carbon foams using the iso2mesh [START_REF] Fang | Tetrahedral mesh generation from volumetric binary and grayscale images[END_REF] package (Fig. 3.1).

2. Alternatively, the packed spheres are fed to in the voro++ [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF] software to create a Voronoi-Laguerre tessellation, that is then refined and stabilized in Surface Evolver [START_REF] Brakke | The surface evolver[END_REF]. The stabilized Voronoi-Laguerre cells, through a process of wall duplication-separation, are transformed into realistic 3D models of high porosity closed cell foams such as cellular plastics (Fig. 3.2).

3. In another variation, cell faces are eliminated and only the skeleton formed by struts is left.

Then, 3D polygonal struts of controlled shape are added following the skeleton and virtually welded at their intersections through a shrink-wrapping [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF] process, allowing to obtain realistic models of high porosity open cell foams such as metal foams (Fig. 3.3). The scheme of the morphology generation method is illustrated in Fig. 3.4, illustrating the sub processes involved, the inputs and the final outputs. In the following various parts are explained, that the obtained structures are periodic, that which is often useful for simulation purposes.

Significant features and innovation relevant to specific parts and specific morphologies will be highlighted in the respective sections.

Close random sphere packing generation.

As seen in the introduction, the process always starts with the generation of a random periodic packing of spheres in a unit cube. The inputs of the algorithm are the number of spheres and their relative diameter distribution. The number of spheres in the initial periodic packing controls the number of cells in the final structure and the diameter of the cells in the final structure is proportional to the diameter of the corresponding spheres. To achieve the final desired metric cell size distribution, the structures are scaled at the end of the generation process. The diameter distributions used are Gaussian unimodal distributions and lognormal unimodal distributions, of varying standard deviation, chosen on a case-by-case basis according to data collected by tomographic analysis. The dispersion of the diameter distribution is controlled through its coefficient of variation 𝐶𝑉 = 𝑠𝑡𝑑(𝑑)/𝑑 ̅ for the Gaussian unimodal distribution and through its geometric coefficient of variation 𝐺𝐶𝑉 = √ 𝑒 𝑣𝑎𝑟(ln(𝑑)) -1 for the lognormal unimodal distribution. During the process of packing, the diameter of each single sphere is varied multiple times, but the ratio 𝑑 𝑖 𝑑 ̅ ⁄ of the diameter of any given sphere di to the average diameter of all cells 𝑑 ̅ stays constant. Therefore, the coefficients of variation of the distribution also stay constant.

The spheres are initially positioned in space with a Random Sequential Absorption (RSA) algorithm [START_REF] Kraynik | Structure of random foam[END_REF] that drops the prescribed number of spheres into space according to the prescribed size distribution, enforcing non-overlapping condition, with a final packing density of 30%. The values of di are randomly picked from a normal or lognormal distribution of mean 1 and prescribed CV or GCV (respectively), then scaled to occupy 30% of the volume of a unit cube:

∑ 𝜋𝑑 𝑖 3 6 = 0.3.
After reaching a density of 30% by RSA, the spheres are packed through a sequential inflationpacking algorithm [START_REF] Zhang | Jamming transition in emulsions and granular materials[END_REF], briefly presented in the scheme of Fig. 3.5. The objective of the algorithm is to reach the maximum packing density of spheres possible without intersections among spheres. This is obtained iteratively. At each step, the spheres are inflated by a small amount (starting at 1%). All sphere pairs are then tested for contact, and for contacting sphere, adjustment vectors that allow to resolve the contact (i.e. make the spheres tangent) are also calculated. Even though the coordinates of the spheres are stored only for spheres whose centers sit inside the bounding box, all the contact pairs evaluation are conducted considering that the domain and the spheres wrap periodically around the 3 coordinate axes, so that the final output is periodic. Fig. 3. 6 shows a simplified 2D representation of the domain. Each sphere is then moved according to the largest adjustment vector associated to it. The contact pairs are then recalculated and the adjustment repeated until no contact pairs are left. Every 10 contact-resolution iterations a small random perturbation is added, to simulate shaking, in order to avoid low density jammed configurations. When after 10000 contact-resolution iterations there are still contacts, the last inflation step is rolled back and the inflation per iteration is reduced, until 0.001% inflation per iteration is reached. The algorithm then stops and the resulting structure is considered to be jammed/packed. The final output of the algorithm is the list of the centers of all spheres and their respective final diameters. First, the spheres are further inflated, keeping the relative radius distribution and the centers unchanged and allowing intersections (Fig. 3.8). The inflation is calculated as to make the sum of volume of all spheres equal to the required porosity (not accounting for intersections): this allow to crudely match the required value of porosity. The resulting geometrical structure is then converted into (binary) voxel format. The voxel structure allows to evaluate the porosity exactly, and to correct the error caused by neglecting the intersections. To fine-tune porosity to the required value, the structure is converted to a real-valued format, following which a 5x5x5 3D Gaussian Blur filter [START_REF] Nixon | Feature extraction & image processing for computer vision[END_REF] is applied. A threshold filter is then applied to transpose the structure again in binary voxel format, taking care to choose the threshold value in such a way as to obtain the exact required final porosity (Fig. 3.9). Finally the voxel structure is meshed using iso2mesh [START_REF] Fang | Tetrahedral mesh generation from volumetric binary and grayscale images[END_REF]. Similar approaches have already been tried with good results in [START_REF] Kırca | Computational modeling of micro-cellular carbon foams[END_REF][76] [START_REF] Chueh | Effective conductivity in random porous media with convex and non-convex porosity[END_REF]. Compared to existing approaches, ours is most similar to the one seen in [31], in that it allows to generate a structure that is periodic in all the three coordinate directions, while only requiring the cell size distribution as additional morphological input, and to obtain a defined value of porosity without having to perform additional manipulations such as pore deletions. However, we eschew the contact laws based approach of [31] in favor of a packing-and-inflation approach that, while less realistic, is considerably simpler and still allows to obtain fairly realistic output morphologies.

Voronoi-Laguerre diagram and Surface Evolver processing

To obtain highly porous closed (Fig. 3.

2) and open cell (Fig. 3.3) structures, such as those seen respectively in plastic and metal foams, further processing of the packed spheres (see Paragraph 3.2) is needed. A periodic 3D Voronoi-Laguerre diagram is generated using Voro++ [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in C++[END_REF], using the packed spheres as seeding spheres of the diagram. The Voronoi -Laguerre diagram partitions the space in polyhedral regions, one for each seeding sphere. For every point in a given region, its minimum distance to the external surface of that region's seeding sphere is lower than its minimum distance from the external surface of any other seeding sphere. It is important to note that Voronoi-Laguerre tessellation guarantees that each seeding sphere is fully contained inside its respective cell (Fig. 3.11). The final output of Voro++ is a list of cells defined in terms of their vertices (Fig. 3.

11).

A polyhedral cell can be obtained from its vertices by calculating the convex hull of said vertices. The resulting data structure is adapted into a Surface Evolver [START_REF] Brakke | The surface evolver[END_REF] input file. The structure is processed in Surface Evolver (Fig. 3.12) combining standard evolution and compression/traction annealing [START_REF] Kraynik | Structure of random foam[END_REF]. Periodicity of the structure is preserved during the Surface Evolver processing.

Treating the structure with Surface Evolver makes it possible to achieve two results:

-The resulting structure is an energetically stable structure rather than just a random structure, that which should make it more realistic. The Surface Evolver refines the initially simple Voronoi mesh, allowing to obtain curved faces [START_REF] Brakke | The surface evolver[END_REF][78] [START_REF] Kraynik | The structure of random foam[END_REF]. The difference in refinement between initial mesh and final mesh can be seen in Fig. 3.12. In addition, during the evolution process vertices and struts that are unstable with respect to Plateau's laws are exploded to obtain energetically stable configurations.

-The cell size distribution can be directly specified and enforced in Surface Evolver, rather than indirectly. In fact, while the Voronoi-Laguerre tessellation guarantees that the seeding spheres are included in the respective cell, it does not guarantee a fixed correspondence between sphere volume and cell volume (as can also be seen in Fig. 3.11). Then, to finely control the final cell size distribution, an additional step is required in Surface Evolver. Each cell is assigned a target volume 𝑉 𝑖 = 𝐾𝑑 𝑖 3 , where di is the diameter of the respective seeding Voronoi sphere, and K is chosen so that the sum of cell volumes equals the volume of the cubic periodic RVE: ∑ 𝑉 𝑖 = 𝑉 𝑅𝑉𝐸 .

The final output of Surface Evolver is a list of cells defined in terms of their faces, struts and vertices, that can be further processed to obtain realistic high porosity closed and open cell structures.

Closed cell plastic foams

To finalize the closed cell structure, first the level of mesh refinement is chosen and the geometry obtained by Surface Evolver (see Paragraph 3.4) is refined to the desired number of elementary triangular facets. Following the refinement step, cells are "separated", that is to say, every face is duplicated and each the two copies is grouped with other connected faces, each group constituting the boundaries of one of the original cells. Each cell-group can then be assigned with a displacement vector: one can see that with an opportunely chosen set of displacement vectors, a wall of chosen thickness between cells can be represented (Fig. 3.13). To obtain walls of uniform thickness, an opportune error function is constructed and minimized. The relevant magnitudes will be hereby presented with reference to Fig. 3.14. To construct the error function, let 𝐹 𝑛 (1) and 𝐹 𝑛 (2) be two corresponding duplicated facets, let then 𝐹𝐶 𝑛 (1) and 𝐹𝐶 𝑛 (2) be the centroids of 𝐹 𝑛 (1) and 𝐹 𝑛 (2) respectively after displacement, let

𝐹𝑁 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 𝐹𝑁 𝑛 (1) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 𝐹𝑁 𝑛 (2)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ be the unit normal of the facets 𝐹 𝑛 (1) and 𝐹 𝑛 (2) , so that |𝐹𝐶 𝑛 (1) 𝐹𝐶 𝑛

⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝐹𝑁 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ | is the distance between the facets measured along their normal, i.e. the local wall thickness, let finally THICK be the target thickness of the walls (chosen according to desired final porosity), the error function is written as:

𝐸𝑅𝑅 ({𝐹𝐶 𝑛 (1) 𝐹𝐶 𝑛 (2) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ } 𝑛=1…𝑁 𝑓 ) = ∑ (|𝐹𝐶 𝑛 (1) 𝐹𝐶 𝑛 (2) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝐹𝑁 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ | -𝑇𝐻𝐼𝐶𝐾) 2 𝑁 𝑓 (3.1)
Where Nf is the number of facets .Being defined as the sum of the squared differences between desired thickness and effective local wall thickness, this function trivially tends to 0 when all walls have the desired thickness, so by minimizing it we can find a set of vectors that allow to best approximate this result. However this error function has an extremely high number of independent vectors, equal to the number of facets. One can use the fact that facets in the same cell have the same displacement to drastically reduce this number. ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ be the displacement associated to the cell-groups that 𝐹 𝑛 (1) and 𝐹 𝑛 (2) are respectively grouped into and 𝐹𝐶 𝑖 be the original centroid of the facets, making reference again to Fig. 3.14 one can also write:

𝐹𝐶 𝑛 (1) 𝐹𝐶 𝑛 (2) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 𝐹𝐷 𝑛 (2) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ -𝐹𝐷 𝑛 (1) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ (3.2) 𝑅𝑅 ({𝐹𝐷 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗ } 𝑛=1…𝑁 𝑐 ) = ∑ (|(𝐹𝐷 𝑛 (2) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ -𝐹𝐷 𝑛 (1) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) • 𝐹𝑁 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ | -𝑇𝐻𝐼𝐶𝐾) 2 N 𝑓 (3.3)
Where Nf is the number of facets and Nc is the number of cells. By this substitution, the number of independent vectors is reduced to the number of cells, which is about 300 times smaller than the number of facets. This error function is minimized using an Interior Point Algorithm, with cell-group

displacement

Open cell metal foams

Fig. 3.17 -Cell structure after Surface Evolver processing and corresponding skeleton.

To finalize the open cell structure, edges with three neighboring cells, i.e. struts, are then isolated to obtain the structure skeleton, used in further processing (Fig. 3.17 where 0 ≤ x ≤ ls is the local abscissa and ls is the length of the strut. The ratio of minimum to maximum diameter, t = 𝑑 𝑠 𝑚𝑖𝑛 / 𝑑 𝑠 𝑚𝑎𝑥 can be controlled (Fig. 3.

19).

The resulting structure constituted by all the struts together is not a continuous mesh, but rather presents a number of self-intersections, namely at the junction points of the struts. The application of a "virtual welding" (shrink-wrapping) process at the junction makes it possible to obtain a continuous mesh, free of self-intersections. This process involves identification of intersecting facets, construction of the convex hull of said facets, then iterative refining and projection of said convex hull into the original polygonal structure, and smoothing of the resulting mesh [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF]. Fig. 3.20 depicts a typical junction before and after the process. As part of the skeleton edges extend to the borders of the bounding box (Fig. 3.17), the 3D struts that have these edges as axes partly extend outside the bounding box. As a final step, these struts are cut along the bounding box, to obtain a structure that is fully periodic and fully contained into the cubic RVE bounding box, with controlled porosity, cells of controlled size distribution and struts of controlled shape (Fig. 3.21).

This model has been employed for radiative heat transfer simulation and validated against tomographical data in recent work by Cunsolo et al. [START_REF] Cunsolo | Radiative properties of irregular open cell solid foams[END_REF] Compared to previous literature [34][35] [15], the process hereby presented allows to finely control the cell size distribution through Surface Evolver and to control the shape of the strut to a degree that had not been realized before in pseudo-random structures.

Introduction

Cellular foams are a key material for many technological applications. Their high porosity (or low relative density) and large specific surface area play an important role from the thermal point of view.

For example, high porosity closed cell polymer foams are used as efficient insulating materials

[1][2][3][4]
. Metal or ceramic foams are being employed in a variety of high temperature applications, such as volumetric solar energy receivers for CSP plants [11], compact heat exchangers [START_REF] Lu | Heat transfer in open-cell metal foams[END_REF], porous radiant burners [9][10] and fire barriers [START_REF] Lu | Thermal transport and fire retardance properties of cellular aluminium alloys[END_REF]. Accurate modeling of thermal properties is obviously highly desirable for the optimization of the performance in these applications. Considering the high porosity (typically in a range from 85% up to 98%), radiative heat transfer contribution can be significant, and in some cases even prevalent over other heat transfer modes [START_REF] Mendes | Measurement and simplified numerical prediction of effective thermal conductivity of open-cell ceramic foams at high temperature[END_REF]. For this reason, a large number of analytical and numerical approaches have been dedicated to the characterization of radiative heat transfer in cellular solid foams. Most studies focus on determining appropriate equivalent continuous medium properties.

Specifically focusing on radiation, the current baseline state-of-the-art approach is based on the utilization of the Radiative Transfer Equation (RTE). General lines on the usage of RTE for radiative transfer can be found for example in textbooks [12][43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF]. While alternative approaches exist, the RTE approach is usually considered sufficiently accurate for most practical cases, if the relevant coefficients (radiative properties) are correctly determined [39][66] [START_REF] Kamiuto | Correlated radiative transfer in packed-sphere systems[END_REF].

As such, most of the literature has been focused on finding efficient and reliable ways to determine radiative properties [13][70] [START_REF] Patel | Evaluation of radiative properties of a representative foam structure using blocked-off region approach integrated with finite volume method[END_REF]. In recent years, numerical methods based on Monte Carlo techniques for the determination of radiative properties are becoming established in order to study either real structures obtained from tomographic imaging or computer generated structures that closely mimic the microstructure of the real foams.

Tancrez and Taine [36] proposed to use the Radiative Distribution Function Identification (RDFI) model and determined radiative properties of spherical packed beds. Zeghondy et al. [22][23] and Petrasch et al. [24] applied the RDFI approach to tomographic data of cellular foam samples. Coquard et al. [25][26] [27] proposed to use an alternative Monte Carlo approach based on mean free path calculation. Cunsolo et al. [15] recently presented a review including an extensive discussion of numerical methods, including Monte Carlo techniques.

Techniques fully based on tomographic data provide satisfactory agreement with experimental data, but their dependence on high quality scans of existing foam samples makes them of limited utility for design purposes. To overcome these limitations, a number of recent studies have sought to digitally reproduce the foam structures using different approaches, including mathematical morphology operations applied on existing tomography data [29][30], simulation of the bubbling process [31],

regular [32][33] and irregular [34][35] Voronoi partitions. By computer generating a number of structures and running numerical simulations [29][30] [35] it is possible to obtain useful results for the optimization of energy transfer. Irregular 3D Voronoi structures seem to be particularly promising for this purpose as they can be described with a limited amount of parameters, are based on wellknown generation methods and approximate the structures of real foam reasonably well [START_REF] Randrianalisoa | Microscale direct calculation of solid phase conductivity of voronoi s foams[END_REF].

In the current paper, a novel methodology is presented that allows the generation of polygonal mesh to represent high porosity open cell foams with high control of a number of geometrical parameters.

The methodology is subsequently applied systematically with individually varying parameters to generate a number of structures. The generated structures are introduced into a Monte Carlo algorithm for the calculation of radiative properties, and especially of the extinction coefficient, β. For each structure, the specific surface area Sv and porosity ε are also calculated, and the normalized extinction [36] is deduced. This procedure makes it possible to determine which parameters have a significant effect on the extinction coefficient and which parameters have a negligible effect, thus defining the inputs required to calculate the extinction coefficient with a given accuracy. Based on this assessment, new analytical relations are given that fit the numerical results with a minimum number of parameters and more accurate than those present in the literature. These relations are expected to be useful for material design purposes.

coefficient β + = 4β•ε ⁄ Sv

Methodology

Digital generation methodology

Digital generation methodology -Presentation

In the present work, a methodology is proposed that makes it possible to generate realistic foam structures. The methodology involves the generation of a pseudo-random lattice of seeding points, the generation of a Voronoi diagram of these points, the stabilization of the resulting cell structure in Surface Evolver [START_REF] Brakke | The surface evolver[END_REF], the addition of polygonal struts along the resulting skeleton, and finally the virtual welding of the polygonal struts at their intersections through a shrink-wrapping [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF] process.

While similar generation methods based on Voronoi diagrams have already been presented in literature [15][34] [35], the current approach does present some distinct features, specifically the use of Surface Evolver and the capability to directly generate intersection-free triangular meshes.

The process requires initially generating a number of seeding points. The corresponding final structure will be a periodic structure containing as many cells as initial seeding points. For this study, a number of 128 initial seeding points, corresponding to 128 cells final structure has been used. This number of cells ensures convergence of Monte Carlo algorithms [15] and allows creating a cubic Kelvin foams with a whole number of Kelvin periodic units (made up by 2 cells), making comparison with regular structures simpler. The points are generated with a Random Sequential Absorption algorithm [START_REF] Feder | Random sequential adsorption[END_REF] that drops equal-sized spheres into space, enforcing non-overlapping condition, with a final packing density around 30%, the centers of the spheres being used as the seeding point. This serves to insure a minimum distance between any two seeding points (Fig. 4.1). -The resulting structure is an energetically stable structure rather than just a random structure, that which should make it more realistic.

-The cell size distribution can be directly specified and enforced in Surface Evolver, rather than indirectly as in approaches using perturbed regular structures [15][35]

The structure is refined and processed in Surface Evolver until a stabilized final structure is obtained.

The size distribution of cells used are Gaussian unimodal distributions of varying standard deviation.

The dispersion of the cell size distribution is controlled through its coefficient of variation 𝐶𝑉 = 𝑠𝑡𝑑(𝑑)/𝑑 ̅ . Edges with three neighboring cells, i.e. Plateau edges, are then isolated to obtain the structure's skeleton, used in further processing (Fig. 4.2). ] where 0 ≤ x ≤ 𝑙 𝑠 is the local abscissa and 𝑙 𝑠 is the length of the strut.

The ratio of minimum to maximum diameter, t = 𝑑 𝑠 𝑚𝑖𝑛 / 𝑑 𝑠 𝑚𝑎𝑥 can be controlled (Fig. 4.4).

The resulting structure constituted by all the struts together is not a continuous mesh, but rather presents a number of self-intersections, namely at the junction points of the struts. The application of a "virtual welding" (shrink-wrapping) process at the junction makes it possible to obtain a continuous mesh, free of self-intersections. This process involves identification of intersecting facets, construction of the convex hull of said facets, then iterative refining and projection of said convex hull into the original polygonal structure, and smoothing of the resulting mesh [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF]. Fig. 4.5 depicts a typical junction before and after the process. The present algorithm presents two distinct features:

1) It makes use of the software Surface Evolver to obtain a more realistic, stabilized cell structure with a finely controlled final cell size distribution.

2) It directly generates the triangular mesh without ever resorting to a voxel representation of the structure, that which allows representing very fine details in large structures.

The second feature of the algorithm is of paramount importance for the feasibility of the next part of the study, where to provide a more realistic description of the structure, both sub-pore scale geometrical features (e.g. strut cross section) and pore scale geometrical features (e.g. pore size distribution dispersion) are considered.

Digital generation methodology -Validation

To verify the capabilities of the proposed methodologies, digitally generated structures have been compared with data obtained through tomography of real metal foam samples. Tomographic data was acquired for four samples by means of micro computed tomography (μCT). General information about the samples is presented in The tomographic data has been analyzed using the free software iMorph [START_REF] Brun | IMorph: A 3D morphological tool to fully analyse all kind of cellular materials[END_REF] to extract the equivalent

Radiative properties calculation methodology

Monte Carlo algorithms for calculation of radiative properties transfer have been extensively investigated in literature. The basic premise of the method is casting a large number of rays inside a digitally represented structure, and applying statistical treatment to the histories of the rays to characterize the radiative properties of the structure. Namely, for a given structure Monte Carlo radiation methods allow to calculate the extinction coefficient β, the scattering coefficient σ and scattering albedo ω = σ / β, the scattering phase function Φ(θ). A comprehensive review of Monte Carlo methods for calculation of radiative properties in cellular media can be found in [15].

In the present paper, radiation propagation is assumed to obey the laws of Geometric Optics. This approximation is considered acceptable for metal and ceramic foams, where the strut diameter is typically in excess of 100µm. The solid phase is considered opaque. The rays are launched from random points in the fluid phase, and are considered extinct (by absorption or scattering) when they intercept the solid surface. Rays that traverse the bounding box without intersection with the solid are launched again from the opposite side of the bounding box, after randomization of position [15].

The extinction coefficient is calculated as the inverse of the mean free path length of the rays:

𝛽 = 1 ∑ 𝑠 𝑛 𝑁 𝑛 𝑁 = 1 𝑠 𝑎𝑣𝑔 (4.1)
It is worthwhile noting that such a characterization of the radiative properties of the medium is only valid if extinction in the medium follows a quasi-Beer-Lambert law. For all the geometries processed in this work, coherency of numerically calculated extinction with a Beer-Lambert law has been verified with an error of less than 1% across the entire range of values. In With respect to scattering, the reflection is assumed to be diffuse and, additionally, reflectivity is considered independent of incidence angle, as it is assumed in several previous works For these reasons, the following analysis focuses on the extinction coefficient, β.

Furthermore, the foams have been considered isotropic in the following analysis. Indeed, anisotropy effects have usually been found to be small in metal foams, with directional values within 10% of average values [25][28]. Furthermore, recent studies on aluminum foams have found anisotropy to be decreasing with increasing PPIs and practically undetectable in foams over 30 PPI [START_REF] Ranut | High resolution microtomography-based CFD simulation of flow and heat transfer in aluminum metal foams[END_REF] such as our samples.

The computations in the following have been executed on grids comprising a number of triangular mesh elements between 6x10 5 and 10 6 , with a memory occupation of 25-50MB per grid. Each computation involved the tracing of 10 6 rays, with an average running time of around 22 minutes on 2.70 Ghz dual core processor.

Results and discussion

Four parameters and their influence have been analyzed in this study: dispersion of cell size distribution (CV = {0%; 2.5%; 5%; 7.5%; 10%; 15%}), variation of strut diameter along the length (t = {0.33; 0.5; 0.66; 0.84; 1.0}), porosity (ε = {98%; 96%, 94%; 92%; 90%; 88%}) and normalized curvature of strut cross section (k = {-0.29; 0; 0.5; 1.0}). A base case (CV = 5%; t = 1.0; ε = 94%; k Our results show that not all parameters considered have comparable influence on radiation. Dispersion of cell size does not appear to affect the extinction coefficient and the normalized extinction coefficient greatly. In fact, all structures considered appear to predict a small interval spanning ±1.5% around the average value (Fig. 4.13).

The small influence of cell size dispersion prompts further investigation. Considering that a number of models in literature, both analytical [17] and numerical [32][33] have used fully regular structures as models to compute the properties of real structures, it is deemed worth of interest comparing a structure of this type to the irregular ones. Thus, an additional structure has been generated using a regular Voronoi partition based on a BCC lattice (Kelvin foam) (Fig. 4.14).

While this single data point is further away from the average, the differences remain small in absolute value: the regular structures show higher normalized extinction coefficient (+2%) and a lower extinction coefficient (-6%) than random structures. This can be justified considering that struts of the regular structure meet at tetrahedral angles, thus minimizing self-shadowing and increasing the former, while on the other hand the regular structure allows the surface to be minimized, thus reducing the latter. It must be noted that long-range ordering effects in regular structures, namely the existence of preferential directions that never meet extinction, are ignored due to the nature of the Monte Carlo algorithms used, due to the choice of randomization of rays exiting the bounding box. This result partially supports the practice of using fully regular structures to estimate radiative properties in cellular media while minimizing required computational power, if some caution is used to correct the small errors incurred. However, the difference is such that the Kelvin structure has been excluded from the following analysis. Strut cross section shape is also found to have large effects on the extinction coefficient and smaller, but not negligible effects on the normalized extinction coefficient. Overall, it is the second most influent parameter among those considered (Fig. 4.17).

Fig. 4.17 -Effects of the variations of k on β and β +

Simplified analytical relations 4.4.1 Analytical relations -Presentation

Based on the above findings, some simplified relations can be devised to calculate the extinction coefficient using a reduced set of parameters.

It can be noted that the normalized extinction coefficient, i.e. the ratio of extinction coefficient and specific surface area, varies in a quite restricted range and is sensitive only to porosity and strut cross section shape. Thus, one can devise a very precise estimation of the extinction coefficient using porosity (ε) (Fig. 4.16) and strut cross section shape (k) (Fig. 4.17) to estimate the normalized extinction coefficient with a small error, then multiplying it by the specific surface area to obtain the value of β. Using our data, the following equation was estimated:

𝛽 = 𝑆 𝑣 4𝜀 1.76[1+0.4(1-𝑘)] (4.4)
The relation is proposed in a form similar to that of Brewster's relations for spherical beds [START_REF] Coquard | Radiative characteristics of opaque spherical particles beds: a new method of prediction[END_REF]. The analytical solution for the independent scattering β = Sv ⁄ 4 is empirically adjusted with a factor depending from porosity and strut cross section shape Physical consistency is assured, as lim 𝜀→1 𝛽 = 𝑆 𝑣 /4 . The strut cross section shape appears as a (1k)

factor that can be thought of as "deviation from circular shape". The proposed relation fits all the data with a maximum error under 3% and a root mean square error under 1% (Fig. 4.18). The form of the relation has been chosen to generalize that typically seen in literature [9], adding dependence from strut cross section shape and variation of strut diameter along its length. Allowing the exponent associated with (1 -ε) to freely vary results in values very close to 0.5, so the square root form has been retained and the exponent imposed in the regression. Strut cross section shape and variation of strut diameter along its length appear as (1k) and (1t) factors that can be thought of as "deviation from circular shape" and "deviation from constant diameter" respectively. The exponent of two associated with (1k) and (1t) is chosen a priori, following from the shape of the curve on Figs.

Analytical relations -Validation

To validate the analytical relations (4.4) and (4.5), tomographic data from the same four samples presented in Section 2 have been used in direct Monte Carlo simulations, and the results have been compared with those obtained applying the analytical relations to morphological data from the same samples. The morphological data required as input have been obtained using iMorph and are listed in The results of the validation can be seen in Table 4.3.

Sample number

β [m -1 ] (MC) β [m -1 ] (4.4)
Relative error (4.4) As expected, relation (4.4) affords the best accuracy, with a maximum error below 2% and a root mean square error of 1.21%. Relation (4.5), on the other hand, affords a maximum error under 3.5% and a root mean square error of 2.11%. In the light of the wide variability of geometrical characteristics of the four samples considered, the agreement between analytical results and Monte Carlo simulations on tomographic data is considered very satisfying for both relations.

β [m -1 ] (4.
In addition, the accuracy of relations ( 1) and ( 2) has been compared with two reference relations from literature. The first is the one given by Brewster [START_REF] Brewster | Volume scattering of radiation in packed beds of large, opaque spheres[END_REF] for spherical beds and uses Sv as a variable:

𝛽 = 𝑆 𝑣 4𝜀 (4.6)
The second is the one given by Glicksman et al. [17] and uses the cell diameter as a variable:

𝛽 = 4.09 √1 -𝜀 𝑑 𝐺 (4.7a)
Considering that the cell representation used in [17] is that of an equivalent dodecahedron and that the mean diameter of the dodecahedron dG is used, an equal volume rescaling from equivalent sphere to dodecahedron has been adopted for the equivalent cell diameter (dG = 1.08d [START_REF] Russ | Practical stereology[END_REF]). Additionally, the calculations in [17] use a correction factor of √2/3 to account for the concavity of the strut cross section typically observed in plastic foams [START_REF] Glicksman | A study of radiative heat transfer through foam insulation[END_REF].

As the samples used in this study do present convex cross sections, this correction factor has been reversed. This ultimately yields:

𝛽 = 3.09 √1 -𝜀 𝑑 (4.7b)
The results of the comparison can be seen in Table 4.4.

Sv dc

Conclusions

In this study, a novel methodology to generate open cell structures along with Monte Carlo numerical methods have used to predict radiative properties of open cell solid foams. Triangular meshes were used to replicate real foams, with finely controlled cell size distribution and strut shape, and smooth, realistic strut intersections free of self-intersecting elements. To demonstrate the capabilities of the algorithm, the geometrical characteristics of the resulting structures have been compared to those of four real foam samples, with satisfying results.

The methodology was subsequently applied, individually varying four parameters (porosity, strut cross section, strut minimum to maximum diameter ratio, dispersion of cell size distribution) in the typical experimentally observed range, to produce a number of structures.

The application of Monte Carlo ray tracing algorithms to the resulting structures made it possible to reduce the number of parameters and to devise and propose two analytical relations to determine the value of the extinction coefficient, fitting all the computed data with small error.

The relations have been validated by comparison with results from direct Monte Carlo computations on four real foam samples with satisfying results, and then compared to existing reference relations from literature. Good absolute accuracy and significant accuracy improvement over existing relations is achieved.

Additionally, simplified forms of the proposed relations have been obtained and their accuracy compared with full form relations and with reference relations. The simplified relations are less accurate than the original ones, but they still represent a significant improvement over reference relations. These relations should be useful for the design of materials e.g. for energy efficiency.

Appendix 4A -Evaluating the normalized curvature (k)

To evaluate the normalized curvature k, it is necessary to establish a relation between k and some stereological measures. In this work, iMorph [START_REF] Brun | IMorph: A 3D morphological tool to fully analyse all kind of cellular materials[END_REF] was used as stereological tool, so in the following a relation is proposed to derive the normalized curvature from measures provided by the software.

The cross-section measurement function available in iMorph has been used. The function automatically identifies struts and then provides geometrical measures for the identified 2D cross sections, namely surface Ss, perimeter Ps, diameter of circumscribed circle do.

A relation between the normalized curvature and a dimensionless ratio of the measures obtained by the software has been established.

To this end, three generated geometries at set values of k = {-0.29, 0.35, 1} have been treated. In addition, this ratio can be easily related to the ratio of the area of the cross section to the area of a simple planar shape enclosing the cross section. With reference to the area of the circumscribed circle of the cross section:

𝜋 𝑆 𝑠 𝐴 𝑐𝑖𝑟𝑐𝑙𝑒 ⁄ = 1.11𝑘 + 1.52 (4A.2)
With reference to the area of an equilateral triangle having the same vertices as the cross section [START_REF] Singh | Modelling radiative heat transfer in packed beds[END_REF]:

3√3 4 𝑆 𝑠 𝐴 𝑡𝑟𝑖𝑎𝑛𝑔𝑙𝑒 ⁄ = 1.11𝑘 + 1.52 (4.A.3)
It must be stressed that, while relations (4A.1)-(4A.3) have been obtained using iMorph, they can be in principle applied to the same geometrical parameters measured using any stereological technique.

Appendix 4B -Evaluating the diameter ratio (t)

To evaluate the diameter ratio t, it is necessary to establish a relation between t and some stereological measures. It is worth noting that the diameter ratio can also be readily measured from SEM scans of the foam [START_REF] Baillis | Thermal radiation properties of dispersed media: theoretical prediction and experimental characterization[END_REF], but in this work, iMorph [START_REF] Brun | IMorph: A 3D morphological tool to fully analyse all kind of cellular materials[END_REF] was used as stereological tool, so in the following a relation is proposed to derive the diameter ratio from measures provided by the software.

The cross-section measurement function available in iMorph has been used, together with a granulometry and segmentation of the solid phase, which allowed identification and measurement of strut junctions. The relevant geometrical measures provided are:

-For the strut cross sections, the diameter of inscribed circle di and the diameter of circumscribed circle do.

-For the strut junctions, the maximum diameter dj,max.

A relation between the diameter ratio and a dimensionless ratio of the measures obtained by the software has been established. The main concern is stability of the measure with respect to variation of other geometrical parameters, namely porosity and strut cross section shape. Four geometries, with porosity and cross section shapes corresponding to those measured for the four samples, have been generated with a constant t = 0.5. Ideally, one would be able to find a ratio that stays constant over the four geometries. The dimensionless ratio (do+di)/(2dj,max) has been picked. It shows very low sensitivity to variations of porosity and strut cross section shape, with all values falling in a ±1% interval (Tab. 4B.1). To relate this ratio to the value of the diameter ratio, three geometries at set values of t = {0.2, 0.5, 1} have been treated, and the dimensionless ratio has been plotted as a function of t (Fig. 4B.1). Fig. 4B.1 -Variation of (do+di)/(2dj,max) with t.

Reference

The relationship is not linear, but a simple power law (plotted in the figure) fits the data points very well:

𝑑 𝑜 + 𝑑 𝑖 2𝑑 𝑗,𝑚𝑎𝑥 = 𝑡 2 3 (4B.1)
It must be stressed that, while relation (4B.1) has been obtained using iMorph, it can be in principle applied to the same geometrical parameters measured using any stereological technique.

Introduction

In literature, the few attempts at comparing results of direct homogenized approaches with results from Direct Monte Carlo simulations [40][48], in cellular media with a semi-transparent solid phase, have revealed substantial discrepancies. In this section, we present the results of our efforts to improve over the current state of the art. The efforts are organized along two main lines: on one hand establishing flexible and efficient reference methods, on the other hand, improving existing homogenized approaches.

In chapter 5, a purely numerical, Direct Monte Carlo Homogenization reference technique is proposed, based on a periodic REV, allowing to simulate radiation within arbitrary cavities with arbitrary boundary conditions and calculate macroscopic radiative quantities (such as transmittance, reflectance, configuration factors etc.) using ray-counting methods typically applied in MCRT for participating media. The main ideas driving the approach are spatial-directional decoupling of the morphological domain from the physical domain and simultaneous ray casting in the superposed domains. The technique is validated against full scale Monte Carlo simulations and compared to the existing HPA model.

In chapters 6 and 7, more accurate and robust versions of the existing Homogeneous Phase Approach and Multi Phase Approach are proposed, in a highly systematic and easily reproducible fashion. In the present work, a new two-pronged strategy is presented to improve the accuracy of conventional homogenized approaches, while retaining much of their simplicity:

1) On one hand, hybrid direct-inverse methods are used to identify radiative parameters. In the particular implementation hereby discussed, the extinction coefficient and scattering albedo are determined directly from free path distribution and relative occurrence of scattering events, while the scattering phase function is characterized as an Henyey-Greenstein and its asymmetry parameter gh is determined through an inverse method. These methods make it possible to determine parameters straightforwardly, while allowing to capture multiplescattering effects that are typically lost on direct methods. In addition, the methods that will be presented in this work use unbounded domains for all calculations, thus moderating the problems tied to the choice of spherical [START_REF] Coquard | Radiative characteristics of opaque spherical particles beds: a new method of prediction[END_REF][57] [START_REF] Coquard | Radiative Characteristics of Beds of Spheres Containing an Absorbing and Scattering Medium[END_REF] or slab-like [START_REF] Subramaniam | Solution of the inverse radiation problem for inhomogeneous and anisotropically scattering media using a Monte Carlo technique[END_REF][55] boundaries (e.g. dependence on boundary geometry).

2) On the other hand, the RTE is modified, with additional equations to take into account the history of the rays, namely their origin, a feature absent in previous models [37][66] [START_REF] Kamiuto | Correlated radiative transfer in packed-sphere systems[END_REF] with the exception of the recent GRTE-MPA [41][42] models, which however impose significant computational complexity and require ad hoc Monte Carlo codes to be solved. In contrast, our modified RTE equations retain a certain simplicity and can be solved by Discrete Ordinates or similar resolution schemes.

Chapter 5 -A Direct Monte Carlo Homogenization method for simulation in arbitrary enclosures.

Introduction

Most literature on radiative transfer in porous media is concerned with determination of appropriate homogeneous models and their relative properties [15][22][36] [37][41] [42]. However, the hypotheses underlying these models are often not fully achieved in the media, and the results of the models are often in disagreement. In addition the few comparisons in literature with direct simulations have given mixed results [40][48]. With the increasing computational power available today, the directly simulating the propagation of radiation at the microscopic scale becomes feasible even for macroscopic systems. Geometrically modeling a large volume of porous material at the micro-scale is prohibitive in most cases, however by using Representative Volume Element (RVE) theory the size of the radiation cavity can be dissociated from the size of the underlying geometrical model, i.e. it is possible to model a large volume of material using only an internal geometrical representation the size of a RVE. While some authors have presented attempts at direct Monte Carlo simulation of radiation propagation [40] that dissociate physical domain from computational domain for some specific configuration, these efforts have been punctual and not systematic nor general. In the following, we present a general Direct Monte Carlo Homogenization (DMCH) methodology that allows to directly simulate radiation heat transfer in an enclosure of any shape, filled with a porous medium, using as inputs only a RVE of the porous medium, the shape of the cavity and the relative boundary conditions.

Morphological domain and physical domain

The RVE can be naturally periodic (such is the case of computer-generated domains such as those of Section 2), or it can be made periodic through simple symmetry operations, as shown in Fig. 5.1 Lambertian point source on a surface. Importantly, the physical domain is described in terms of its axes x and y (for this 2D case, in a 3D case an additional z axis is required). In this case, the quantity of interest to evaluate can be the configuration factor from the point source to the walls, or the fraction of radiation absorbed in the medium, etc. It must be stressed that, while for the sake of clarity this particular configuration has been chosen, the method lends itself equally well to enclosures of any shape, with boundary conditions of any kind and to any directional distribution of radiative intensity.

In addition, it must be remarked that surface and volume sources can be simply modeled as In the ray casting process, a large number (typically ≥ 10 6 ) of points are chosen randomly inside the RVE (computational domain). In addition, for each point a random Cartesian reference system is defined: ignoring symmetry, this requires choosing a random vector in 2D space (the second being its perpendicular), or two random perpendicular vectors in 3D space (the third being their cross product). Additionally, for problems that present local azimuthal symmetry in 3D space (e.g. azimuthally symmetrical emission distribution, which turns out to be the case for many real problems), it is possible to choose just one random vector. The random reference system is the physical reference system for that ray, and thus all the calculations relative to the physical domain are performed within that system (Fig. 5.3), while the reference system of the computational domain stays fixed and unchanged at all times. Finally, the direction of the ray is chosen according to the prescribed directional distribution, while taking into account the individual reference system of each ray: this can be seen in Fig. 5.3, as the relative orientation of the Lambertian directional distribution with respect to the random physical reference system stays constant, while it changes with respect to the reference system of the computational domain.

Ray casting process

At each iteration, all rays are cast inside the RVE (computational domain), augmented with a parallelepiped bounding box and with the physical domain, rotated and translated according to its random reference system. It's important to stress that the coordinates of the source points of the rays with respect to the physical reference system are completely independent of their coordinates with respect to the computational domain. Treatment of interaction with the material is detailed in a number of works [15][22][36] [40]. At each iteration, the propagation of a ray can be terminated by different kinds of events:

-Absorption by the material. In this case, after additional calculations (e.g. energy exchange) are executed, the ray is removed.

-Scattering by the material. The new direction of the ray is determined, after which, if necessary, additional calculations (e.g. path length logging) are executed, then finally the ray is enqueued for the next iteration.

-Crossing of the bounding box. The ray is periodically transported on the other side of the RVE and enqueued for the next iteration (Fig. 5.4). It's important to note that the position of the ray is reset with respect to the computational domain, but not with respect to the physical reference system.

-Interaction with the physical domain: the events prescribed by the boundary conditions (reflection, absorption, etc.) are simulated, and if necessary additional calculations (e.g. energy propagation, path length logging, ray counting) are performed, then according to the physics the ray is removed (absorption) or enqueued for the next iteration (reflection).

The process is iterated until all rays have been removed. Then, according to the specific problem, the logged quantities and distributions can be used to calculate the desired variables according to standard MCRT techniques, e.g. ray counting can be used to calculate configuration factors. Fig. 5.5

synthesizes the evolution of a single ray from inception to extinction. It is important to note that this method makes it possible to execute the entire ray casting process while keeping a single copy of the computational domain in memory. Like other Monte Carlo methods, it also lends itself very well to parallelization.

Validation and comparison with Homogeneous Phase Approach

The proposed Direct Monte Carlo Homogenization (DMCH) methodology is validated by crosscomparing it to full MCRT simulation in a fully represented volume of porous medium and with a conventional HPA method calibrated on the RVE.

The hemispherical transmittance Trans, hemispherical reflectance Refl and total absorption Abs through a slab are evaluated (Fig. 5.6). Where N = 10 6 is the total number of incident rays cast and NT, NR and NA are the number of rays that traverse the slab, are reflected towards the incoming direction or are absorbed in the slab respectively.

It is worth noting that NT + NR + NA = N.

Due to its simplicity in terms of representation, a medium constituted by Semi-Transparent Identical Non Overlapping Spheres immersed in a transparent medium is chosen as reference. Two distinct levels of porosity ε, 98% (Fig. 5.7) and 75% (Fig. 5.8), are employed. For each level of porosity, two volumes are generated, a 3x3x3 periodic RVE containing 100 spheres and a 30x30x12 non-periodic RVE containing 40000 spheres. In both cases the spheres are added to the volume by Random Sequential Absorption (RSA). The relatively small size of the RVE is necessary to make the size of the mesh for the full domain Monte Carlo simulation manageable.

To minimize differences due to sampling, 10 slabs of 30x30x6 slabs are cut out of the 30x30x12 volume. For each slab, a radiation flux is imposed on a 6x6 center region on the bottom (Fig. 5.9).

Picking this center region allows to minimize the number of rays escaping the domain from the sides, which is contained below 0.5% of the total for all the simulations. For each domain, 10 5 rays, for a grand total of 10 6 rays are cast from this region and reflectance, transmittance and absorbance are calculated by simple ray counting (Eq. 5.1a-c). The results are averaged across the 10 slabs. For the new DMCH method, a simulation is run using the 100-sphere RVE as RVE and a 30x30x6 parallelepiped as reference physical domain, with casting and counting of 10 6 rays according to the same criteria seen above (Eq. 5.1a-c).

For the Homogenous Monte Carlo method, first the homogeneous parameters are determined using a conventional approach outlined in [40], with 10 6 rays used in the parameter identification process and using the 100-sphere RVE as input. Then a 30x30x6 parallelepiped domain is considered, with casting and counting of 10 6 rays according to the same criteria seen above (Eq. 5. 1a-c).

In all the simulations, possible scattering effects due to the phase transition at the bottom or top boundary are not considered, i.e. phase continuity at the top and bottom boundaries is assumed.

Finally, the error of the homogenized methods is defined as:

ERROR = |TRANS -trans| + |REFL -refl| + |ABS -abs| (5.2)
Where the uppercase represents quantities obtained with the full domain Monte Carlo and the lowercase represents quantities obtained with the homogenization methods.

Finally, various values of the refractive index n2 and the absorption coefficient α2 of the semitransparent phase are considered. The results are summed up in Two main observation can be made in the light of the results shown in the table:

1) The Direct Monte Carlo homogenization method makes it possible to obtain deviations consistently below 2% with respect to the full domain MCRT. Residual errors can be attributed to sampling and to the small size of the RVE.

2) At very high porosities, the hypotheses of the HPA are satisfied (namely, that the medium can be represented as dispersed, point-like scatterers far from each other), so both homogenization approaches turn out to be satisfying. However, at low porosities, HPA hypotheses fail and its errors can be very large, while the errors of the Direct Method stay consistently low.

Computational aspects

The final mesh used for the 3x3x3 RVE containing 100 spheres is composed of 10000 triangular elements. The mesh used for the 30x30x12 full volume is composed of 4000000 triangular elements, with an average of 2000000 triangular elements in each 30x30x6 slab. The simulations are run on a single i7 4790K CPU with 8GB of RAM.

The computation times for the various configurations are reported in A few observations can be made. The proposed method produces a significant decrease of computational time, up to about 5-fold, compared to a Full Domain Monte Carlo. This can be attributed to the sensible reduction of the size of the mesh and the consequent increase in efficiency of the ray casting process. Additionally, the memory footprint is reduced 200-fold. In fact, the Direct Monte Carlo runs in fixed memory irrespective of domain size: this makes it very interesting for larger domains that cannot be represented directly in memory and also makes massive parallelization of the algorithm very simple.

Compared the Direct Monte Carlo Homogenization to the Homogeneous Phase approach using effective radiative properties, we observe that the latter is about 50 to 100 times faster in execution, but if the time for parameter identification is included, the difference in speed is reduced to a factor of 3 to 4. Considering additionally that the absolute values of the computation times are quite small, this makes the Direct Monte Carlo homogenization a viable alternative to the Homogeneous Phase Approach in a number of applications.

Chapter 6 -Improved Homogeneous Phase Approach (HPA+)

Introduction

In the following, we'll be explaining the bases of an Improved Homogeneous Phase approach. First, we'll discuss and example elucidating why the origin of the ray can be of significance with respect to its scattering extinction "history", then we'll propose a modification of the classic HPA approach to take into account such effects, allowing to increase accuracy while retaining simplicity, finally we will flesh out our model in detail, explaining the process of determination of its parameters and testing it.

Ray history effects and general setup

Let's consider for our example an element of phase 2 immersed in phase 1, with n2 > n1. Let this element be a shell of fixed thickness (Fig. 6.1): for some materials, such as closed cell plastic foams this is actually a representative element of the structure. With reference to Fig. 6.1, for rays coming from phase 1 (in dark red),

𝑛 2 𝑛 1
sin 𝛼′ 1 = sin 𝛼 1 ≤ 1, so all rays coming from phase 1 and traversing phase 2, also traverse the 2->1 interface.

On the contrary, for rays emitted inside phase 2 (or incident therein, due to boundary conditions),

we'll have:

• A fraction where

𝑛 2 𝑛 1
sin 𝛼 3 ≤ 1 (in green), that pass the 2->1 boundary with a deviation due to refraction.

• A fraction where

𝑛 2 𝑛 1
sin 𝛼 2 > 1 (in blue), that will be reflected multiple times until they're eventually absorbed. This fraction can be calculated as

𝐹 = √1 - 𝑛 1 2 𝑛 2 2
⁄ of the locally incident radiation [43], so it can be deduced that it can be significant (over 50% of the radiation) for values of the refraction index typically found in semi-transparent plastics and ceramics.

Similar observations can be made observing a sphere, or other simple geometrical forms. It is also worth noting that the effects described will be less pronounced in domains with highly irregular boundaries: as incidence angles are less correlated and more random, the difference between rays coming from inside and rays coming from outside is reduced.

For HPA radiation homogenization models, this can have important consequences if porosity is below 95%. In a standard HPA model all rays are mixed together, losing specificity tied to their origin.

However, as seen in the example just shown, the interaction of rays originating in phase 2 can be dramatically different than for the rays originating in phase 1.

We hereby propose an extension of the basic HPA model that allows to take into account a large part of these effects, decomposing them in two sub-effects that enrich the basic HPA model: The aforementioned coefficients intervene to modify the structure of the standard HPA model. The modification only pertains to the fraction of radiation that is emitted in the solid phase or thereby enters the domain from outside, while the rest of the radiation is treated as usual.

Indeed, in addition to the standard homogenized phase, the enriched model has an additional phase, labeled trapped phase, with a separate RTE, that takes into account the trapping and absorption effects (a). The RTE of the trapped phase is fully characterized by its absorption coefficient αt, i.e. the trapped phase is characterized as a purely absorbing phase, for reasons we'll see in the next. A fraction C of the radiation entering the domain in the solid phase will be assigned to the trapped phase and propagate according to its separate RTE. The rest is assigned to the standard homogenized phase. It is important to note the standard homogenized phase and the trapped phase are not coupled, that which makes the model significantly simpler than other two equations models such as the MPA. To take into account the additional scattering effects in (b), all the radiation entering the domain from outside in the solid phase undergoes a transformation of its directional distribution, according to a Henyey-Greenstein phase function ΦHG(θ) that is fully characterized by its asymmetry coefficient gHG. To model this effect, it is necessary to modify the form of the boundary conditions. This method allows to approximate the multiple scattering events that can take place before absorption in the solid or scattering in the fluid as a single scattering event (see Fig. 6.2). This is a useful approximation for high porosity media because the mean free path associated with propagation in the solid phase is much smaller than the homogenized mean free path. It is also important to note that the volumetric emission term for radiation is isotropic, so that modifying it by applying the directional distribution transformation has no net effect: for this reason the terms associated to the transformation will only appear in boundary conditions, such as imposed directional intensity or emitting/reflecting walls. As part of the radiation entering the domain in the solid phase is assigned to the standard non-trapped homogenized phase, this contribution too (characterized by coefficient 1 -C) will have modified boundary conditions. Exiting radiation (green) is approximated as entering the homogeneous phase directly with a distribution of angles modified using a Henyey-Greenstein function, with asymmetry factor gHG chosen so that the distribution of angles θ from the discrete scale is preserved. Absorbed radiation (blue) is approximated as entering the trapped phase with a distribution of directions modified in the same fashion, while the absorption coefficient of the trapped phase αt is chosen to preserve the distribution of the distance from entry to absorption s, so that the possibly tortuous path of the rays is approximated as a linear path from point of entry to point of absorption. It must be noted that while the distributions of θ for the two radiation fractions may differ in principle, they're treated as one and the same for the sake of simplicity. The modification of direction distribution of radiation in the boundary conditions justifies our choice of treating the trapped phase itself as purely absorbing.

The approximation should be satisfying for mid-high levels of porosity (85-95%), and afford a significant improvement over conventional HPA models, while preserving much of their simplicity. Calculate ghg (6.8) Obtain modified HPA equations and boundary conditions (6.9a-c, 6.11a-b, 6.12, 6.14a-b, 6.15)

Full presentation of the HPA+

6.3.1

Inverse fitting

6.3.3

Calculate C and αt (6.26, 6.27)

6.3.2

Obtain compatibility conditions on ghg, C and αt (6.20, 6.

The novel proposed method can be articulated in three parts, which are labeled (6.3.1-3) (Fig. 6.3) and can be summarized as follows:

6.3.1) Determination of HPA coefficients βh, ωh, gh, with rays cast from the fluid phase.

6.3.2) Determination of adjustment coefficient gHG for scattering effect in solid phase 6.3.3) Modification of equations and calculation of adjustment coefficients C and αt for trapping effect in solid phase.

It is worth underlining that, during the determination of the coefficients of the homogenoues phase, rays are originated only from points in the fluid phase, as in [40].

In the course of the process, heavy use of Monte Carlo Ray Tracing simulations is made. Further details and discussion on MCRT for radiation heat transfer can be found in [15], some adaptations relevant to semi-transparent media are discussed in [40] and in Chapter 3. With respect to design choices discussed in [15], periodic wrapping at the boundaries, which was shown to have the best convergence, is used for the simulation, without necessity of randomizing the positon upon wrapping, because the domains are periodic themselves, thus ensuring phase continuity.

The three phases (6.3.1-3) just summarized and schematized in Fig. 6.3 are detailed in the following. The process of determining suitable equivalent properties of the homogeneous medium starts with a Monte Carlo Ray Tracing simulation. As implied before, all rays for this initial simulation are cast from the fluid phase. Similarly to what discussed in [40], with reference to Fig. 6.4, absorption and scattering events are stored, namely the path length associated with each event is added to an accumulator Lext, while the number of events are counted in counter Nabs and Nscatt. This allows to calculate βh and ωh simply as:

Ray tracing from the fluid phasehybrid determination of coefficients

𝛽 ℎ = 𝑁 𝑎𝑏𝑠 + 𝑁 𝑠𝑐𝑎𝑡𝑡 𝐿 𝑒𝑥𝑡 (6.1) ω ℎ = 𝑁 𝑠𝑐𝑎𝑡𝑡 𝑁 𝑎𝑏𝑠 + 𝑁 𝑠𝑐𝑎𝑡𝑡 (6.2)
A specific adaptation of our algorithm compared to those already presented in literature is that each ray is tracked until it is finally absorbed by the medium: this means that each ray can give rise to multiple scattering logging events.

Additionally, for each absorption event (once for each ray), the distance from the source point to the point of final absorption is also logged in the distribution sabs. This is necessary for the inverse analysis.

Another specific feature is that all path lengths in the case of non-absorption inside the solid phase (scattering events) are logged at the points of exit of the ray from the solid phase.

A further adaptation is that, in case a scattering event results in a scattering angle lower than 2°, the scattering path length is not logged: this allows to implicitly take into account scattering with strong forward peaks (e.g. in case of parallel walls), without having to resort to using ad hoc phase functions such as the Delta-Eddington phase function. The scattering and extinction coefficient obtained in this way are the same that one would obtain applying the transport approximation of the scattering phase function [START_REF] Dombrovsky | Radiation heat transfer in disperse systems[END_REF]. A total of N = 10 6 rays are launched and tracked.

Having determined βh and ωh using the direct approach just outlined, an alternative inverse approach is hereby proposed for the determination of the scattering phase function. In our approach, a simple Henyey-Greenstein, characterized by a single asymmetry coefficient gh phase function is postulated to be sufficient to correctly capture the scattering. To determine the asymmetry coefficient of said phase function, an inverse fitting approach is used. A number of Monte Carlo HPA simulations are run with βh and ωh already obtained and different values of gh. As in the case of Monte Carlo micro structural simulations, each ray is tracked up to its absorption, and the distance from the source point to the point of absorption sabs is logged.

The distribution of sabs is used to calculate the locally absorbed fraction of intensity at a certain distance s from the source point:

𝐼 𝑎𝑏𝑠 (𝑠) 𝐼 0 = ∑ [𝑠 - ∆𝑠 2 ≤ 𝑠 𝑎𝑏𝑠 𝑛 ≤ 𝑠 + ∆𝑠 2 ] 𝑁 𝑖=1 𝑁∆𝑠 (6.3)
Where N is the total number of rays cast, Δs is a discretization interval opportunely chosen to guarantee acceptable smoothness and [] are Iverson brackets:

[x] = { 1 0 𝑖𝑓 𝑥 𝑖𝑠 𝑡𝑟𝑢𝑒; 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. (6.4)
This magnitude is compared between the MC micro structural simulation and the HPA, then a particle swarm algorithm [START_REF] Qi | Inverse radiation analysis of a one-dimensional participating slab by stochastic particle swarm optimizer algorithm[END_REF] is used to minimize squared differences (i.e. to obtain best curve fit), with gh as independent variable (Fig. 6.5).

The error function to minimize can be written:

𝐸𝑅𝑅(𝑔 ℎ ) = ∫ ( 𝐼 𝑎𝑏𝑠 * (𝑠) 𝐼 0 - 𝐼 𝑎𝑏𝑠 𝐻𝑃𝐴 𝑔 ℎ (𝑠) 𝐼 0 ) 2 𝑑𝑠 ∞ 0 (6.5)
Where 𝐼 𝑎𝑏𝑠 𝐻𝑃𝐴 𝑔 ℎ (𝑠) indicates the locally absorbed fraction calculated using the particular determination of the HPA model that is given by a certain value of gh and 𝐼 𝑎𝑏𝑠 * (𝑠) indicates the locally absorbed fraction calculated by direct Monte Carlo analysis. Starting from the HPA model presented in the previous section, we develop our approximation. As already noted in the standard HPA, the rays used to calculate βh, ωh, gh, originate from the fluid phase.

Ray tracing from the solid phaseevaluation of coefficient ghg

As such, the resulting HPA model will satisfactorily capture the behavior of rays emitted in the fluid phase, but not of rays emitted in the dense phase. These differences in behavior are captured in our model through an additional trapping effect and an additional scattering effect. To evaluate the coefficients associated to these effects, we use Direct MC simulations at the micro scale.

Another MC micro structural simulation is executed with N = 10 6 rays starting exclusively from the solid phase.

With respect to Monte Carlo assumptions, they are the same used before, but all rays start from the solid phase. The rays are followed up to their final absorption or their exit from the solid domain, whichever comes first (Fig. 6.6). This makes it possible to obtain information about ray directions at scales comparable to the scattering scale of the solid phase. With reference to Fig. 6.6, for each ray, an angle θ is logged. For rays passing into the fluid, θ is the angle between its final direction vector and the initial direction vector 𝑖 ⃗ . For absorbed rays, θ is the angle between the final direction vector and the vector 𝑠 that goes from the source point to the point of absorption.

We underline once again that, under our approximation, the complex path taken by the ray at the microscale is approximated as a straight line between its source point and its final absorption point (Fig. 6.2).

The scattering effects for both trapped and non-trapped rays are taken into account into gHG, which is also calculated using the information collected in the simulation. Using the collected values of θ a scattering distribution P*(θ) can be simply calculated as:

𝑃 * (𝜃) = ∑ [𝜃 - ∆𝜃 2 ≤ 𝜃 𝑛 ≤ 𝜃 + ∆𝜃 2 ] 𝑁 𝑛=1 𝑁∆𝜃 (6.6)
Where N is the total number of rays cast, [] are Iverson brackets and Δθ is a discretization interval opportunely chosen to guarantee acceptable smoothness. We observe that ∫ 𝑃 * (𝜃)𝑑𝜃 ≅ 1 2π . We note once again that, as the emission term is isotropic, its directional distribution is invariant with respect to scattering, so the scattering effect will only be applied to boundary conditions. It is worthwhile remembering that the mean free path associated with the additional scattering effect in the solid is much smaller than the scattering mean free path associated with the standard HPA (Fig. 6.7), so we can consider it to be negligible. The additional scattering effect can then be completely decoupled from the usual scattering in the standard homogeneous phase and approximated as a scattering event taking place when radiation enters the domain (Fig. 6.7).

We approximate the measured scattering distribution through a Henyey-Greenstein phase function.

𝛷 𝐻𝐺 (𝜃) ≅ 2π𝑃 * (𝜃) (6.7)

By choosing an asymmetry factor gHG equal to that of the directly simulated scattering distribution:

𝑔 𝐻𝐺 = 1 2 ∫ 𝑃 * (cos 𝜃) cos 𝜃 1 -1
𝑑(cos 𝜃) (6.8) Fig. 6.8 shows the computed distribution and the corresponding Heyney-Greenstein phase function. ). In the next section, it will be shown that it is possible to obtain more rigorous values of C and αt by imposing energy conservation conditions on the equations. In our tests, the value of C and αt calculated directly at the micro-scale never differed more than 10% from the values obtained through the energy conservation condition, which comforts us about the consistency of our definitions.

Modified equations and calculation of C and αt.

Finally the HPA equations can be rewritten to take into account these observations. First, the constitutive equations for both the standard homogenized phase and the trapped phase will be presented. Then, the necessary modifications to the boundary conditions will be defined.

As for the constitutive equations, without losing generality, on can postulate n2 > n1 and write:

Ω • ∇I ℎ = [ 𝑓 1 𝑛 1 2 𝛼 1 ⏞ 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 1 + 𝑓 2 (1 -𝐶)𝑛 2 2 𝛼 2 ⏞ 𝑁𝑜𝑛-𝑡𝑟𝑎𝑝𝑝𝑒𝑑 𝑒𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2 ] 𝐵 0 ⏞ 𝑇𝑜𝑡𝑎𝑙 𝑒𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑖𝑛 ℎ𝑜𝑚𝑜𝑔𝑒𝑛𝑒𝑜𝑢𝑠 𝑝ℎ𝑎𝑠𝑒 -𝛽 ℎ 𝐼 ℎ + 𝜎 ℎ 4𝜋 ∫ I ℎ (Ω′)Φ ℎ (Ω ′ , Ω)𝑑Ω′ 4𝜋 (6.9a) Ω • ∇I 𝑡 = 𝑓 2 𝐶𝑛 2 2 𝛼 2 𝐵 𝑣 ⏞ 𝑇𝑟𝑎𝑝𝑝𝑒𝑑 𝑒𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2
-𝛼 𝑡 𝐼 𝑡 (6.9b)

I = I ℎ + I 𝑡 (6.9c)
Here, the subscripts 1 and 2 respectively indicate properties of the two phases, f1 (fluid) and f2 (solid)

indicating the volumetric fractions, α1 and α2 indicating the absorption coefficients. βh, ωh, Φh are the properties of the homogeneous medium, calculated according to section 6.3.1. C and αt are the adjustment coefficients for the trapping effect, which as anticipated will be calculated in the next.

To more easily make sense of the equations, it is useful to refer to synthetic Fig. 6.9. As one can see from Eq. (6.9c), the total local intensity is given as the sum of homogenized intensity Ih, with contributions for both the solid dense phase f1 and the fluid phase f2, and trapped intensity It, only in the solid phase f2. The homogenized intensity equation (Eq. 6.9a) is written in a form very similar to standard HPA (Eq. 1.1), the only difference being the emission term, which is directly calculated by summing the emission contributions of the fluid phase in its entirety and of the non-trapped fraction of the solid phase.

As anticipated, the trapped intensity equation (Eq. 6.9b) is simpler than standard HPA: the phase itself is modeled as purely absorbing (without scattering), and actual scattering effects taking place are taken into account by modifying the boundary conditions with the additional scattering effect modeled through the coefficient gHG, as explained in subsection 6.3.1. Its emission term is proportional to the total emission in the solid phase 𝑓 2 𝑛 2 2 𝛼 2 𝐵 𝑣 and to the trapping factor C. An important point is that in most cases αt >> αh, so It can be ignored at long distance, considerably simplifying the equations. It is interesting to observe that the sum of emission of the two equations equals exactly the effective emission taking place in the medium.

Fig. 6.9 -Synthetic figure depicting the effects on the discrete scale and the corresponding approximations at the homogenized scale. In black, the radiation entering the domain in the fluid phase, assigned to the homogenized phase. In blue and green, the trapped and non-trapped fractions of the radiation entering the domain in the solid phase.

As already discussed, the additional scattering effect and its associated phase function ΦHG do not appear explicitly in any of the constitutive equations, because the volumetric emission terms are directionally isotropic they are thus not affected by the transformation of the directional distribution.

The additional scattering effect will appear only in boundary conditions when the radiation source term has a directional component.

To be able to apply the new equations to real problems, one must be able to write boundary conditions accordingly. In the following, the two most commonly used boundary conditions (namely, incident directional intensity on an open boundary and diffusely emitting and reflecting wall) will be rewritten for the new equations. Neglecting possible reflection effects due to solid-fluid interfaces at the boundary, it is written:

𝐼(𝑟 𝑤 + , Ω 𝑤 ) ⏞ 𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦, 𝑖𝑛𝑠𝑖𝑑𝑒 = 𝐼(𝑟 𝑤 -, Ω 𝑤 ) ⏞ 𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦, 𝑜𝑢𝑡𝑠𝑖𝑑𝑒 = 𝐼 𝑤 (6.10)
In our formulation it becomes:

𝐼(𝑟 𝑤 + , Ω) ⏞ 𝐻𝑜𝑚𝑜𝑔𝑒𝑛𝑒𝑜𝑢𝑠 𝑝ℎ𝑎𝑠𝑒 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦, 𝑖𝑛𝑠𝑖𝑑𝑒 = 𝐼 𝑤 𝑓 1 [Ω = Ω 𝑤 ] ⏞ 𝑃ℎ𝑎𝑠𝑒 1 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 + 𝐼 𝑤 𝑓 2 (1 -𝐶) ⏞ 𝑃ℎ𝑎𝑠𝑒 2 𝑛𝑜𝑛-𝑡𝑟𝑎𝑝𝑝𝑒𝑑 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 Φ 𝐻𝐺 (Ω 𝑤 , Ω) 4𝜋 ⏞ 𝐴𝑑𝑑𝑖𝑡𝑖𝑜𝑛𝑎𝑙 𝑠𝑐𝑎𝑡𝑡𝑒𝑟𝑖𝑛𝑔 (𝑖𝑛𝑠𝑖𝑑𝑒) {𝑛 𝑤 • Ω > 0} (6.11a) 𝐼(𝑟 𝑤 + , Ω) ⏞ 𝑇𝑟𝑎𝑝𝑝𝑒𝑑 𝑝ℎ𝑎𝑠𝑒 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦,𝑖𝑛𝑠𝑖𝑑𝑒 = 𝐼 𝑤 𝑓 2 𝐶 ⏞ 𝑃ℎ𝑎𝑠𝑒 2 𝑡𝑟𝑎𝑝𝑝𝑒𝑑 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 Φ 𝐻𝐺 (Ω 𝑤 , Ω) 4𝜋 ⏞ 𝐴𝑑𝑑𝑖𝑡𝑖𝑜𝑛𝑎𝑙 𝑠𝑐𝑎𝑡𝑡𝑒𝑟𝑖𝑛𝑔 (𝑖𝑛𝑠𝑖𝑑𝑒) {𝑛 𝑤 • Ω > 0} (6.11b) 𝐼(𝑟 𝑤 -, Ω) ⏞ 𝑇𝑜𝑡𝑎𝑙 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦, 𝑜𝑢𝑡𝑠𝑖𝑑𝑒 = 𝐼(𝑟 𝑤 + , Ω) ⏞ 𝑇𝑜𝑡𝑎𝑙 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦, 𝑖𝑛𝑠𝑖𝑑𝑒 + 𝑓 2 𝐼 𝑤 ⏞ 𝑇𝑜𝑡𝑎𝑙 𝑝ℎ𝑎𝑠𝑒 2 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 Φ 𝐻𝐺 (Ω 𝑤 , Ω) 4𝜋 ⏞ 𝐴𝑑𝑑𝑖𝑡𝑖𝑜𝑛𝑎𝑙 𝑠𝑐𝑎𝑡𝑡𝑒𝑟𝑖𝑛𝑔 (𝑜𝑢𝑡𝑠𝑖𝑑𝑒)
{𝑛 𝑤 • Ω < 0} (6.12)

Where [] are Iverson brackets, while I(rw -, Ω) and I(rw + , Ω) represent the left and right (external and internal) limit values of total intensity (Eq. 6.9) at the boundary.

In this case, as the radiation is directional, so additional scattering effects are significant (Fig. 6.10)

and are modeled by the diffusion term ΦHG. The entirety of radiation incident on the solid phase (proportional to the volumetric fraction f2) undergoes the additional scattering effect (Fig. 6.10, blue arrows). In the homogenized phase (Eq. 6.11a), only a part of the radiation undergoes the additional scattering effect, proportional to the non-trapped fraction of the solid phase, while the rest, proportional to the volumetric fraction of the fluid phase, is not scattered (Fig. 6.10, red arrows). In the trapped phase (Eq. 6.11b) the entirety of radiation undergoes the additional scattering effect. One can easily observe that for gHG < 1 the scattering implies a reduction of the energy injected in the domain. As Eq. (6.12) shows, the scattering creates a discontinuity in the values of the intensity around the open boundary, with different values for 𝐼(𝑟 𝑤 -, Ω) and 𝐼(𝑟 𝑤 + , Ω), due to partial reflection of the incoming radiation outside the domain (Fig. 6.11, green arrows) summing itself with the radiation coming out of the domain (Fig. 6.10, pink arrows). It is worth noting that this phenomenon is not due to reflection of radiation at the boundary, but rather to reflections inside the domain at very short distance, and approximated in our treatment as if they were at the boundary.

The second boundary condition discussed is a diffusely emitting and reflecting surface [43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF] of emissivity εs and reflectivity ρs at position rw and with local normal nw, which is usually written: In our formulation it becomes:

𝐼(𝑟 𝑤 + , Ω) = 𝜀 𝑠 𝐵 0 ⏞ 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑓𝑟𝑜𝑚 𝑤𝑎𝑙𝑙 + 𝜌 𝑠 𝜋 ∫ 𝐼(𝑟 𝑤 + , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′ 𝑛 𝑤 •Ω ′ <0 ⏞ 𝑅𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚 𝑤𝑎𝑙𝑙 {𝑛 𝑤 • Ω > 0}
I ℎ (𝑟 𝑤 + , Ω) = 𝑓 1 (𝜀 𝑠 𝑛 1 2 𝐵 0 + 𝜌 𝑠 𝜋 ∫ 𝐼(𝑟 𝑤 , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′ 𝑛 𝑤 •Ω ′ <0 ) ⏞ 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚 𝑤𝑎𝑙𝑙 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 1 + 𝑓 2 (1 -𝐶) ⏞ 𝑁𝑜𝑛-𝑡𝑟𝑎𝑝𝑝𝑒𝑑 𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐼 𝑠 (𝑟 𝑤 , Ω) {𝑛 𝑤 • Ω > 0} (6.14a)
I 𝑡 (𝑟 𝑤 + , Ω) = 𝑓 2 𝐶 ⏞ 𝐼 𝑠 (𝑟 𝑤 , Ω) 𝑇𝑟𝑎𝑝𝑝𝑒𝑑 𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 {𝑛 𝑤 • Ω > 0} (6.14b)
Where, for the sake of compactness:

𝐼 𝑠 (𝑟 𝑤 + , Ω) ⏞ 𝑆𝑢𝑝𝑝𝑜𝑟𝑡 𝑡𝑒𝑟𝑚 𝑓𝑜𝑟 𝑒𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2 = (𝜀 𝑠 𝑛 2 2 𝐵 0 + 𝜌 𝑠 𝜋 ∫ 𝐼(𝑟 𝑤 , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′ 𝑛 𝑤 •Ω ′ <0 ) ⏞ 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚 𝑤𝑎𝑙𝑙 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2 × 1 4𝜋 ∫ Φ 𝐻𝐺 (Ω ′ , Ω)𝑑Ω ′ 𝑛 𝑤 •Ω ′ >0 1 - 𝜌 𝑠 4𝜋 2 ∫ ∫ Φ 𝐻𝐺 (Ω ′ , Ω′′)|𝑛 𝑤 • Ω ′′ |𝑑Ω ′ 𝑑Ω ′′ 𝑛 𝑤 •Ω ′ >0 𝑛 𝑤 •Ω ′′ <0 ⏟ 𝑀𝑢𝑙𝑡𝑖𝑝𝑙𝑒 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛𝑠 𝑎𝑡 𝑡ℎ𝑒 𝑤𝑎𝑙𝑙 ⏞ 𝐴𝑑𝑑𝑖𝑡𝑖𝑜𝑛𝑎𝑙 𝑠𝑐𝑎𝑡𝑡𝑒𝑟𝑖𝑛𝑔 {𝑛 𝑤 • Ω > 0} (6.15)
In this case too, as the radiation is directional (being uniform only on the positive hemisphere with respect to the surface), additional scattering effects are significant (Fig. 6.11) and are modeled by the scattering term ΦHG and I is the total incident intensity (Eq. 6.9).

The support term I s (Eq. 6.15), allows to model the additional scattering effect associated with the interface between the wall and the solid phase. As in the case of a diffusely emitting and reflecting surface the intensity is uniform for all the positive hemisphere, the factor that transforms the directional distribution through the scattering term ΦHG can be taken out of the integral. Just as in the case of open boundary, the scattering transformation through ΦHG implies a reduction of the energy injected in the domain, with a part of the radiation (proportional to

1 4𝜋 ∫ Φ 𝐻𝐺 (Ω′, Ω′′)dΩ′ 𝑛 𝑤 •Ω ′ >0
for any direction Ω ′′ where 𝑛 𝑤 • Ω ′′ < 0) being scattered towards the wall. However, as the wall can be reflective, in this case a part (proportional to

𝜌 𝑠 𝜋 ∫ 𝐼(𝑟 𝑤 , Ω ′′ )|𝑛 𝑓 • Ω ′′ |𝑑Ω ′′ 𝑛 𝑤 •Ω ′′ <0
) of the energy backscattered towards it goes back into the domain again. This gives rise to a geometric series of ratio

𝜌 𝑠 4𝜋 2 ∫ ∫ Φ 𝐻𝐺 (Ω ′ , Ω′′)|𝑛 𝑤 • Ω ′′ |𝑑Ω ′ 𝑑Ω ′′ 𝑛 𝑤 •Ω ′ >0 𝑛 𝑤 •Ω ′′ <0
, the sum of which gives the adjustment factor to account for multiple reflections. Fig. 6.12 -Illustration of scattering at the boundary for a diffusely emitting and reflecting wall at the discrete (left) and homogeneous scale (right), up to the second reflection. Coloring is coherent between the two sides to help understand discrete-continuum correspondence. On the homogenous scale, the directional distribution transformation sequence is also reported.

As the boundary backscattering phenomenon has been referenced to in the description of both boundary conditions and its justification may appear unclear, a more detailed explanation is in order. Fig. 6.12 illustrates the correspondence between discrete and continuum scale for scattering at the boundary for the specific case of a diffusely emitting and reflecting wall, but the general principle is valid for any boundary condition. Refraction and reflection at the short range can give rise to rays (in blue in the figure) which are oriented towards the boundary itself. As this short range scattering is approximated as a point-wise scattering in the continuum model, it corresponds to the presence of negative lobes in the scattering phase function ΦHG. In the particularly complex case of a partially reflecting wall, this gives rise to multiple reflections (Eq. 6.15). In the case of an open boundary, seen previously, it gives rise to the discontinuity of the values of intensity at the border (Eq. 6.12).

The term I s is used to calculate the homogenized (6.14a) and trapped (6.14b) intensity: for the former, it is weighed by the non-trapped fraction of the solid phase and summed with a standard intensity for the fluid phase, while for the latter it represents the entirety of the radiation, proportional to the trapped fraction of the solid phase.

To further clarify how the "homogenized phase", "trapped phase" and "additional scattering" treatments are applied to different fractions of radiations, it is useful to refer to Table 6 As was already mentioned, it is very useful to introduce some constraints on the values of the coefficients and on the form of the equations. This can be done by applying Kirchhoff's law of radiation [43][44] [START_REF] Brewster | Thermal radiative transfer and properties[END_REF]. Applying Kirchhoff's law to Eq. (6.9a-b) gives:

𝑓 1 𝑛 1 2 𝛼 1 + 𝑓 2 (1 -𝐶)𝑛 2 2 𝛼 2 ≅ 𝛼 ℎ [𝑓 1 𝑛 1 2 + 𝑓 2 (1 -𝐶)𝑛 2 2 ] (6.16a) 𝛼 2 ≅ 𝛼 𝑡 (6.16b)
Where the ≅ sign takes into account the fact that boundary conditions for the reformulated problem are not exactly the same as a conventional model of radiation in participating media. These equations can be developed more rigorously using the exact boundary conditions on the wall to obtain compatibility conditions on the additional terms C, αt and gHG, ultimately allowing to calculate values of C and αt as anticipated.

Starting from the trapped phase, we consider a black wall boundary and Eq. (6.15b), obtaining:

I 𝑡 (𝑟 𝑤 , Ω) = 1 4𝜋 𝑓 2 𝑛 2 2 𝐵 0 𝐶 ∫ Φ 𝐻𝐺 (Ω ′ , Ω)𝑑Ω ′ 𝑛 𝑤 •Ω ′ >0
{𝑛 𝑤 • Ω > 0} (6.17)

And then the total emitted energy:

𝑄 𝑡 𝑒𝑚𝑖𝑡𝑡𝑒𝑑 = 𝑓 2 𝑛 2 2 𝐵 0 𝐶 1 4𝜋 ∫ ∫ Φ 𝐻𝐺 (Ω ′ , Ω)|𝑛 𝑤 • Ω|𝑑Ω ′ 𝑛 𝑤 •Ω ′ >0 𝑛 𝑤 •Ω>0 𝑑Ω = 𝜋𝑓 2 𝑛 2 2 𝐵 𝑣 𝐶 𝑔 𝐻𝐺 + 1 2 (6.18)
Where

1 4𝜋 ∫ ∫ Φ 𝐻𝐺 (Ω ′ , Ω)|𝑛 𝑤 • Ω|𝑑Ω ′ 𝑛 𝑤 •Ω ′ >0 𝑛 𝑤 •Ω>0 𝑑Ω = 𝜋 𝑔 𝐻𝐺 +1
2 is a numerical solution, that can be easily verified for g = {-1, 0, 1}. And considering Eq. (6.9b) for a surface delimiting a semi-infinite plate of isothermal non-scattering medium, the incident radiation is From which one can see that:

𝑄 𝑡 𝑖𝑛𝑐𝑖𝑑𝑒𝑛𝑡 = ∫ ∫ 𝑓 2 𝑛
𝛼 2 = 𝛼 𝑡 𝑔 𝐻𝐺 + 1 2 (6.20)
For the homogenized phase, considering a black wall boundary and Eq (6.15a), we obtain:

I ℎ (𝑟 𝑤 , Ω) = 𝑓 1 𝑛 1 2 𝐵 0 + 1 4𝜋 𝑓 2 𝑛 2 2 𝐵 0 (1 -𝐶) ∫ Φ 𝐻𝐺 (Ω ′ , Ω)𝑑Ω ′ 𝑛 𝑤 •Ω ′ >0 {𝑛 𝑓 • Ω > 0} (6.21)
From which the total emitted energy can be derived:

𝑄 𝑡 𝑒𝑚𝑖𝑡𝑡𝑒𝑑 = 𝜋𝑓 1 𝑛 1 2 𝐵 0 + 𝜋𝑓 2 𝑛 2 2 𝐵 0 (1 -𝐶) 𝑔 𝐻𝐺 + 1 2 (6.22)
And considering Eq. (6.9a) for a surface delimiting a semi-infinite plate of isothermal non-scattering medium, the incident radiation is:

𝑄 𝑡 𝑖𝑛𝑐𝑖𝑑𝑒𝑛𝑡 = 𝜋 [𝑓 1 𝑛 1 2 𝛼 1 + 𝑓 2 𝑛 2 2 𝛼 2 (1 -𝐶)]
𝛼 ℎ 𝐵 0 (6.23)

Equating incident (Eq. 6.22) and emitted (Eq. 2.21) energy we finally obtain:

𝑓 1 𝑛 1 2 𝛼 1 + 𝑓 2 𝑛 2 2 𝛼 2 (1 -𝐶) = 𝛼 ℎ [𝑓 1 𝑛 1 2 + 𝑓 2 𝑛 2 2 (1 -𝐶) 𝑔 𝐻𝐺 + 1 2 ] (6.24)
While it would be possible to determine all coefficients numerically, numerical incertitude can bring about values that do not fully respect Eqs. (6.20) and (6.24). We prefer to use a numerically calculated value of gHG and use Eqs. (6.20) and (6.24) to obtain C and αt.

There are multiple reasons for this. On one hand, gHG is the only additional term that appears in both compatibility conditions. On the other hand, one can see that performing the adjustment this way, it is unlikely to obtain physically meaningless values of C and αt, while the reverse is not true.

From Eq. (6.20) one can see that the only case where the value of αt cannot be determined from gHG is when the latter is exactly equal to -1. On the contrary, any value of αt smaller than α2 gives rise to a physically meaningless gHG > 1.

From Eq. (6.24), differentiating with respect to C and gHG and rearranging we obtain:

𝑑𝐶 𝑑𝑔 𝐻𝐺 = 𝛼 ℎ (1 -𝐶)

Numerical results

Testing is conducted on for 3 different physical configuration and 5 different porous morphologies, which will be presented in detail in the following. The geometries are characterized by a single geometrical parameter, a length L, which is a characteristic length from which all geometrical measurements of the virtual sample can be deduced (see Figs. 6.13-15). The structure is considered to be constituted by a non-interacting fluid phase and a semi-transparent solid phase. The solid material is characterized in terms of its index of refraction n2, and its coefficient of absorption α2: In the determination of properties, we chose to cast rays from both the solid and the fluid phase according to their volumetric fraction, because this marginally improves the results compared to the benchmark.

2) HPA-INV: a standard HPA model, using Eq. (1.1) with standard BCs (Eqs. 6.10, 6.13) and determination of properties through the hybrid inverse-direct coefficient fitting discussed in section 2.4A.

3) HPA+: an extended model (HPA+), using Eqs. (6.9a-c) and modified BCs for open boundaries (Eqs. 6.11a-b, 6.12) and opaque walls (Eqs. 6.14, 6.11a-c), with properties determined following the entire process discussed in Paragraph 6.3

For each of the 5 morphologies, the three models are compared in terms of maximum and average error compared to the benchmark, with averages taken separately for each physical configuration over the 9 triplets characterizing the solid material optical properties. Each physical configuration has a specific definition of error, depending on the boundary conditions and the main quantities of interest, which will be illustrated in the following (Eqs. 6.30, 6.32, 6.38). This error is also used to check for numerical convergence of each method with respect to the number of rays launched, for both the Direct Monte Carlo method and the homogenized methods. Numerical convergence is considered achieved when, using the defined errors to compare the results of two distinct runs of the same methodology, they're repeatably within 1% from each other.

The three different physical configurations are hereby described:

1) Simulation of collimated radiation on an unbounded slab of thickness L and evaluation of hemispherical transmittance Trans, hemispherical reflectance Refl and total absorption Abs (Fig. 6.13). Where uppercase identifies Full Direct Monte Carlo and lowercase identifies homogenized methods. Where N = 10 6 is the total number of rays cast, Nn is the number of rays impacting on the n th wall (Fig. Where uppercase identifies Full Direct Monte Carlo and lowercase identifies homogenized methods.

3) Simulation of radiative heat exchange between two black walls at fixed temperatures TA = 300K

and TB = 600K and evaluation of heat flux Q (Fig. 6.15). The hypothesis that the temperature is a function of linear position T = T(x) is made. The domain is divided in 50 equal-thickness layers, and temperature is assumed to be uniform in each layer.

Additionally, to accelerate convergence, rather than using a direct energy packet exchange calculation of temperature [43] + 𝑁 𝑛,𝐴 + 𝑁 𝑛,𝐵 = 𝑁 𝑛 .

The 5 morphologies considered cover a relatively vast variety of porous material, and have been digitally generated making use of the techniques illustrated in Section 2. All the morphologies are periodic in the three coordinate directions and are made up of 128 spheres/cells. This number of cells has been considered to satisfactorily represent the random structure while not imposing very high computational loads [15][85]. A general presentation of all the seven morphologies is given in Table 6.2. The generation parameters relative to each morphology are also given.

Comments

The results show a trend, fairly consistent for all the porous morphologies and physical configurations considered, where the Inverse HPA method improves over the Classic HPA method, and the HPA+ method improves over the Inverse HPA method.

The adoption of the Inverse fitting of the scattering phase function appears justified in the results: in all cases the Inverse Method outperforms the Classic Method in terms of error, and the dramatic improvement seen in pretty much all morphologies (and more markedly morphology #5) appears sufficient to justify the adoption of such a method over the standard one in most cases, also considering that no modification of the equations is needed to do so.

Results also show that, if further accuracy is required, modification of the equations according to what has been presented as HPA+ can provide a further dramatic improvement of accuracy. The HPA+ method makes it possible to reduce errors up to 8 times, compared to the Classic HPA.

The error reduction is especially marked, both in absolute and relative terms, when the porosity is lower, which is consistent with the idea that the traditional HPA methods neglect some effects tied to the presence of rays incident on the solid phase and that the HPA+ correction considers these effects.

Importantly, the HPA+ method appears to be able to consistently guarantee errors below 10% for all the considered morphologies and set of properties on Test #3, which considers a very directly engineering-relevant metric, i.e. the heat flux through a wall filled with the porous medium, while the other proposed methods give errors up to 30%.

Given these results, it would appear advisable to use Inverse fitting of HPA scattering phase function in all cases, considering that minimal additional calculations are required, compared to the direct method. In addition, when lower (<90%) levels of porosity are involved, and high accuracy is required, the proposed HPA+ method can provide a very significant improvement over existing methods, at the price of a small complication in the equations.

Chapter 7 -Improved Multi Phase Approach (MPA+) 7.1 Introduction

In the following, we'll be explaining the bases of an Improved Multi Phase Approach (for information on MPA, see Subparagraph 2.2.1). Following the scheme established with the HPA+, we'll first discuss the significance of ray histories, with reference to specific problems that arise in the context of MPA. Then we'll propose a modification of the classic MPA to take into account such effects, allowing to increase accuracy while retaining simplicity, finally we will flesh out our model in detail, explaining the process of determination of its parameters and testing it.

Ray history effects and general setup of MPA+

Recalling Fig. 6.1 from Chapter 6, most of the same considerations apply. Indeed, even if in a MPA model each phase is characterized separately, the same core problem illustrated before present themselves more or less in the same fashion, namely, rays in the solid phase behave very differently according to their origin (in the solid phase itself, or from the fluid phase). However, the fact that the solid phase is treated explicitly has some impact. The phenomenon of "trapping" identified previously is still relevant, and must be considered along with the presence of a corresponding "trapped phase. In Fig. 7.1 above, an extreme case is depicted, where a large fraction of rays originating in the solid phase are "trapped", while most rays coming from outside pass through. This can be modeled by assigning to the "trapped phase" a part of radiation usually assigned to the solid phase. The trapped phase is modeled with an additional constitutive equation, in the same fashion as was done for the Improved HPA, but in this case it is a scattering and absorbing phase. The trapped phase is then fully characterized by the trapped fraction C1 and the radiative coefficients αt, σt and asymmetry factor gt.

Concerning the scattering effect, on the other hand, it is already taken into account explicitly in the standard MPA. Indeed, the MPA model already includes a scattering phase function Φ 21 (Ω ′ , Ω), and such phase function can be expected to not depend significantly on the origin of the ray. Looking at Fig. 7.1 it is clear that the distribution of angles (α3⇒α3') associated with the non-trapped fraction must be similar to the distribution (α1'⇒α1) associated with the fraction coming from the fluid phase, while rays intersecting the interface at angles beyond the critical angle will be totally internally reflected in both cases, thus having no effect on the scattering phase function Φ 21 (Ω ′ , Ω). Not having to explicitly model the additional scattering in the two main phases, it is convenient to model the trapped phase as a scattering-absorbing phase: this simplifies considerably the writing of the boundary conditions for the modified equations.

Finally, looking at Fig. 7.1 it is easy to highlight another effect that is specific of the MPA context.

In MPA, a unique scattering path length characteristic of the solid phase must be calculated.

Comparing the intensity coming from the fluid dense phase (red) with the non-trapped intensity originating in the solid phase (green), one can easily see that some geometries (such as the one in figure) can cause a systematic difference between the average path lengths associated with these intensities (in this case, the latter's path length is on average half the former's) in the solid phase. As only one set of coefficients (and a corresponding path length) can be associated to the solid phase, it is necessary to find a way to model this difference without altering the coefficients. We choose to add an additional coefficient C2. The non-trapped fraction (1-C1) of radiation originating in the solid phase is distributed between the two homogenized phases according to C2, in such a way that the mean transmission path length in the solid is preserved in the transition from discrete to homogenous scale, according to the following principle:

• Radiation assigned to the homogenized solid phase will have a mean transmission path length

𝑠 2 = 1 𝛼 2 +𝜎 21
in that phase, with 𝜎 21 determined using rays coming from the fluid, possibly longer than the actual path length for rays originating in the solid.

• Radiation assigned to the homogenized fluid phase will have a mean transmission path length in the solid s2 = 0.

• The actual path length 𝑠 2 * can then be approximated as a weighted average of the standard

homogenized path length 𝑠 2 = 1 𝛼 2 +𝜎 2
and the zero path length s2 = 0, according to a conveniently chosen weighing coefficient C2, so that:

𝑠 2 * = 𝑠 2 (1 -𝐶 2 ) (7.1)
It is worth noting that this effect does not depend on directionality, so it will intervene to modify both the emission terms of the constitutive equations and the boundary conditions. To take into account the variability in mean path length in the solid according to the origin, the non-trapped fraction of radiation is not entirely assigned to the solid phase, but it is distributed between the two standard homogenized phases according to the coefficient C2, as explained above.

As was the case for the HPA, the aforementioned coefficients intervene to modify the structure of the standard MPA model. The modification only pertains to the fraction of radiation originating in the solid phase, while the fraction of radiation originating in the fluid phase is treated according to standard MPA convention.

The approximation should be satisfying for a large interval of porosities, and afford a significant improvement over conventional MPA models, while preserving much of their simplicity. The novel proposed method can be articulated in four parts, which are labeled (7.3.1-4) (Fig. 7.4) and can be summarized as follows: It is worth underlining that, during the determination of the coefficients of the homogenoues phases, rays are originated only from points in the fluid phase, similarly to [40].

Full presentation of theMPA+

In the course of the process, heavy use of Monte Carlo Ray Tracing simulations is made. Further details and discussion on MCRT for radiation heat transfer can be found in [15], some adaptations relevant to semi-transparent media are discussed in [40] and Chapter 2. With respect to design choices discussed in [15], periodic wrapping at the boundaries is used for the simulation, without randomization of position.

The four phases (7.3.1-4) just summarized and schematized in Fig. 7.4 are detailed in the following. A specific adaptation of our algorithm compared to those already presented in literature is that each ray is tracked until it is finally absorbed by the medium: this means that each ray can give rise to multiple scattering logging events.

Additionally, for each absorption event (once for each ray), the distance sabs from the source point to the point of final absorption and the angle θabs between its original direction vector and the vector going from the source point to the point of absorption are stored (Fig. 7.5). This is necessary for the inverse analysis.

Having determined σij using the direct approach just outlined, an alternative inverse approach is hereby proposed for the determination of the scattering phase functions. In our approach, the scattering phase function are defined as Henyey-Greenstein phase function, each one , characterized by a single asymmetry coefficient, for a total of 4 coefficients g11, g12, g21, g22 (collectively gij). It can also be noted that, by optical reversibility, g12 = g21 [38], which reduces to 3 the number of coefficients to be found.

A number of Monte Carlo MPA simulations by numerically solving by means of Monte Carlo simulations Eqs. (2.4a-b) using the σij already obtained and different values of gij. As in the case of Monte Carlo micro structural simulations, each ray is originated in the fluid phase, and each ray is tracked up to its final absorption, logging the distance from the source point to the point of final absorption sabs and the angle θabs between its original direction vector and the vector going from the source point to the absorption point.

The distribution of sabs can be reshaped to represent the polar distribution of absorbed radiation. This is simply obtained by counting all the rays falling in a given interval:

𝐼 𝑎𝑏𝑠 (𝑠, 𝜃) 𝐼 0 = ∑ [𝑠 - ∆𝑠 2 ≤ 𝑠 𝑎𝑏𝑠 𝑛 ≤ 𝑠 + ∆𝑠 2 ∧ 𝜃 - ∆𝜃 2 ≤ 𝜃 𝒆𝒙𝒕 𝑛 ≤ 𝜃 + ∆𝜃 2 ] 𝑁 𝑛=1 𝑁∆𝑠∆𝜃 (7.4)
Where N is the total number of rays cast, [] are Iverson brackets and Δs and Δθ are discretization intervals opportunely chosen to guarantee acceptable smoothness.

The resulting values of 𝐼 𝑎𝑏𝑠 (𝑠,𝜃) 𝐼 0 are compared between the MC micro structural simulation and the MPA-MC simulation, then a particle swarm algorithm [START_REF] Qi | Inverse radiation analysis of a one-dimensional participating slab by stochastic particle swarm optimizer algorithm[END_REF] is used to minimize squared differences (i.e. to obtain best curve fit), with gij as independent variables (see Fig. The average path of a ray in the solid up to passage in the fluid s* can be then calculated as:

𝑠 * = 𝐿 𝑠𝑐𝑎𝑡𝑡 𝑁 𝑜𝑢𝑡 (7.6) In this case, the calculation of the coefficients is made significantly less straightforward due to the presence of the explicit treatment of the dense phase in the standard MPA equations. However, it is still possible to obtain the coefficient in a fairly simple manner.

Writing C2 as the weighing coefficient that preserves the transmission mean path of non-trapped rays:

𝑠 * = (1 -𝐶 2 )𝑠 2 (7.7)
Where s2 is the path length in the solid associated with homogeneous phase 2:

𝑠 2 = 1 𝛼 2 + 𝜎 21 (7.8) Then C2 can be obtained as:

𝐶 2 = 1 - 𝑠 * 𝛼 2 + 𝜎 21 (7.9)
It would be possible to obtain C1 from its phenomenological definition by observing its relation with the fraction of rays that ever reach the fluid phase:

𝑁 𝑜𝑢𝑡 𝑁 = (1 -𝐶 1 ) [𝐶 2 + (1 -𝐶 2 ) 𝜎 21 𝛼 2 + 𝜎 21 ] (7.10)
Where first factor reflects the fact that none of the trapped radiation ever reaches the fluid phase and the second factor reflects the fact that a fraction C2 of the radiation immediately reaches the fluid phase and a fraction (1 -C2) propagates through the solid phase.

However, in subsection 7.3.4, it will be shown that it is possible to obtain a more rigorous value of C1 by imposing an energy conservation condition on the equations. In our tests, the value of C1 calculated by using Eq. (7.10) never differed more than 5% from the value obtained through the energy conservation condition, which comforts us about the consistency of our definitions. Rather than considering the trapped phase a purely absorbing one as in the HPA+, in this case we characterize it as an absorbing and scattering phase. The absorption coefficient, on the base of purely physical reasoning, can be set equal to that of the solid phase, αt = α2. To determine the other properties, once again a hybrid inverse approach is applied, quite similar to the one used to determine the properties of the homogeneous phase in the HPA+. A MC micro structural simulation is executed with N = 10 6 rays starting exclusively from the solid phase. In this case, the rays are followed up to their final absorption in solid phase, while rays exiting the solid phase are not considered. Internal scattering events are stored, namely the path length associated with each event is added to an accumulator Lscatt while the number of events is counted in a counter Nscattt. This allows to calculate σt simply as:

𝜎 𝑡 = 𝑁 𝑠𝑐𝑎𝑡𝑡 𝐿 𝑠𝑐𝑎𝑡𝑡 (7.11) Additionally, for each absorption event (once for each ray), the distance from the source point to the point of final absorption is also logged in the distribution sabs. The scattering phase function of the trapped phase is then characterized as a Heyney-Greenstein phase function and its asymmetry coefficient gt is calculated according to the method already illustrated in Subparagraph 6.3.1, i.e. by inverse fitting between the micro-structural simulation and a number of HPA-MC simulations.

The distribution of sabs is used to calculate

𝐼 𝑎𝑏𝑠 (𝑠) 𝐼 0 = ∑ [𝑠 - ∆𝑠 2 ≤ 𝑠 𝑎𝑏𝑠 𝑛 ≤ 𝑠 + ∆𝑠 2 ] 𝑁 𝑖=1 𝑁∆𝑠 (7.12)
Where N is the total number of rays cast, [] are Iverson brackets and Δs is a discretization interval opportunely chosen to guarantee acceptable smoothness.

This magnitude is compared between the MC micro structural simulation and the HPA, then a particle swarm algorithm [START_REF] Qi | Inverse radiation analysis of a one-dimensional participating slab by stochastic particle swarm optimizer algorithm[END_REF] is used to minimize squared differences (i.e. to obtain best curve fit), with gt as independent variable.

The error function to minimize can be written:

𝐸𝑅𝑅(𝑔 ℎ ) = ∫ ( 𝐼 𝑎𝑏𝑠 * (𝑠) 𝐼 0 - 𝐼 𝑎𝑏𝑠 𝐻𝑃𝐴 𝑔 𝑡 (𝑠) 𝐼 0 ) 2 𝑑𝑠 ∞ 0 (7.13)
Where 𝐼 𝑎𝑏𝑠 𝐻𝑃𝐴 𝑔 𝑡 (𝑠) indicates the locally absorbed fraction calculated using the particular determination of the HPA model that is given by a certain value of gt and 𝐼 𝑎𝑏𝑠 * (𝑠) indicates the locally absorbed fraction calculated directly using the Monte Carlo simulation at the discrete scale.

Modified equations and calculation of C1

Similarly to what has been done for the HPA+, we'll now present the modified equations. First, the constitutive equations for both the standard homogenized phases and the trapped phase will be presented. Then, the necessary modifications to the boundary conditions will be defined.

As for the constitutive equations, without losing generality, on can postulate n2 > n1 and write: It is useful to explain the equations starting from the bottom.

Ω • ∇I 1 = 𝑛 1 2 𝛼 1 𝐵 1 0 + 𝑓 2 (1 -𝐶 1 ) 𝑓 1 𝐶 2 𝑛
In Eq. (7.14d) we can clearly see the total intensity being given by the weighted sum of the three phases, and the trapped phase being weighed against the solid phase through the coefficient C1.

Indeed, the non-trapped solid phase is reduced to a virtual volume fraction of f2(1-C1).

Eq. (7.14c) is in this case a conventional HPA equation. It is interesting to observe that the sum of emission of the three equations equals exactly the effective emission taking place in the medium.

The first two terms represent the modified forms of the conventional MPA terms. There are two distinct modification with respect to the conventional form (see also Paragraph 7.2):

1) The terms that account for the partial assignment of radiation in phase 1 rather than phase 2.

2) The terms that account for the presence of the trapped phase. Fig. 7.2 can be useful for guidance when interpreting the equations.

By applying the treatment seen in Subparagraph 2.2.2 to Eq. (7.14a) or, identically, to Eq. ( 7.14b),

we can obtain a compatibility condition that, as previously mentioned, can be used to calculated the value of C1. Recalling the definitions:

𝛽 1 = 𝛼 1 + 𝜎 11 + 𝜎 12 (7.15a) 𝛽 2 = 𝛼 2 + 𝜎 22 + 𝜎 21 (7.15b)
And assuming the thermodynamic equilibrium conditions:

𝑇 1 = 𝑇 2 = 𝑇 (7.16) 𝐵 1 0 = 𝐵 2 0 = 𝐵 0 (7.17) 4𝜋𝑛 1 2 𝛼 1 𝐵 0 = 𝛼 1 ∫ 𝐼 1 (Ω′)𝑑Ω′ 4𝜋 (7.18a) 4𝜋𝑛 2 2 𝛼 2 𝐵 0 = 𝛼 2 ∫ 𝐼 2 (Ω′)𝑑Ω′ 4𝜋 (7.18b)
And writing the radiative equilibrium condition, for Eq. (7.14a): From which we can directly obtain C1 as:

𝑓 1 ∇ • ∫ Ω𝐼 1 𝑑Ω 4𝜋 = 4𝜋𝑓 1 𝑛 1 2 𝛼 1 𝐵 0 + 4𝜋𝑓
𝐶 1 = 1 - 𝜎 12 𝑓 1 𝑛 1 2 𝑓 2 𝑛 2 2 (𝛼 2 𝐶 2 + 𝜎 21 ) (7.22)
The choice of obtaining C1 from C2 and not the contrary can be justified by observing that to calculate C1 directly (Eq. 7.10) C2 is required, while the contrary is not true (Eq. 7.9). Additionally, differentiating (Eq. 7.21) and rearranging we obtain: , so that adjusting C1 we run less risks of falling outside the physically meaningful range.

Having chosen to model the trapped phase as an absorbing and scattering phase, there's no additional scattering effects to be modeled as modified boundary conditions, therefore the writing of equations is relatively straightforward. In the following, the two most commonly used boundary conditions (namely, incident directional intensity on an open boundary and diffusely emitting and reflecting wall) will be rewritten for the new equations. These boundary conditions are useful to compute the numerical results in the following. Note that the factor C1 for the trapped phase does not appear explicitly in (7.24b-c), because it is already present in Eq. (7.14d). When reflected at the boundary, radiation either in the trapped phase or in phase 2 is considered to be re-emitted randomly in the solid phase, and thus undergoes the subdivision in trapped phase, phase 2 and in phase 1 already discussed, hence why the support term I2+t appears in all of the equation.

To further clarify how the "homogenized phase", "trapped phase" and "assignment to phase 1" treatments are applied to different fractions of radiations, it is useful to refer to Table 7

Numerical results

Testing is conducted on with conditions identical to those proposed in Paragaph 6.4, the only difference being the addition of a sixth morphology, characterized by open cell structure and low porosity (ε = 0.7)

As a benchmark, the radiation problems are solved using the Direct Monte Carlo method illustrated in Chapter 5. Using this benchmark as a base, five different homogenized models are compared:

1) HPA: a standard HPA model [22][40], using Eq. (1.1) with standard BCs for open boundaries and opaque walls (Eqs. 6.10, 6.13), and determination of properties according to the method shown in [40], i.e. determination of extinction coefficient β and scattering albedo ω with Eqs.

(6.1-2) and direct calculation of scattering phase function [40]. It must underlined that, in the property determination, we chose to cast rays from both the solid and the fluid phase according to their volumetric fraction, because this appears to marginally improve the results.

2) HPA+: an extended homogeneous phase model, using Eqs. (6.9a-c) and modified BCs for open boundaries (Eqs. 6.11a-b, 6.12) and opaque walls (Eqs. 6.14, 6.11a-c), with properties determined following the entire process discussed in Paragraph 6.3.

3) MPA: a standard MPA model [37][38], using Eqs. When this has been found to be the case, the most favorable result has been chosen.

Comments

The results show a trend, fairly consistent for all the porous morphologies and physical configurations considered, where the standard HPA and MPA method are outperformed by the more advanced alternatives.

With respect to the classical methods, it is noted that the HPA results are really acceptable only for a handful of morphologies, namely very low density morphologies (#2, #6), while the MPA fails very badly in morphologies with strong multiple scattering effects, namely "Closed cell foam" type morphologies (#5, #6). In these structures, , due to the fact that the solid walls are almost everywhere parallel, the scattering angle of a given solid⇒fluid transmission event is very strongly correlated to the scattering angle of the previous fluid⇒solid transmission event for the same ray, an effect that the MPA model fails to capture. The MPA also produces somewhat inaccurate results for the other morphologies.

The MPA-GRTE consistently outperforms the MPA method, is quite accurate overall in most cases and appears to be the most accurate method overall for low-porosity morphologies of the type "Open cell foam" (#3, #4). However, it inherits the weakness of the standard MPA with respect to morphologies where multiple scattering effects are very strong, namely "Closed cell foam" type morphologies (#4, #5). This effect is so strong that even the standard HPA outperforms the MPA-GRTE for these morphologies.

Carrying from the previous results, HPA+ consistently improves over HPA. Comparing it to MPA or MPA-GRTE, one notes that the accuracy of HPA+ is comparable or better than MPA in all morphologies considered, barring the "Open cell foam" (#3, #4) type, but worse than MPA-GRTE, barring the Closed cell foam" type morphologies (#5, #6). More generally, there appears to a correlation between porosity and accuracy of the HPA+ method, with less porous morphologies yielding less and less accurate results: as such, this method should be mainly used for high porosity morphologies.

Finally, MPA+, as can be seen by the results, consistently improves over both MPA and HPA+, showing the highest accuracy in all morphologies considered, barring the "Open cell foam" (#3, #4) type. Even for the latter, the accuracy is deemed satisfying and a significant improvement over the standard MPA. Importantly, the MPA+ is the only one of the methods considered that appears to produce fully satisfying results (maximum error below 10%) for all the morphologies proposed, regardless of porosity or multiple scattering phenomena. In addition, while the coefficient calculation process is relatively less straightforward than the one necessary for the standard MPA, requiring the use of hybrid inverse methods, the final result is a relatively compact collection of coefficients and the equations are comparable in complexity to those of the standard MPA, and lend themselves to solution by standard techniques thus making the MPA+ a promising alternative to standard MPA.

General conclusion and future prospects

To correctly model radiation heat transfer in porous cellular foam, both accurate radiation models and accurate morphological models of the structure are necessary. In this work, we tried to advance knowledge on both these fronts.

We presented a framework for the parametric digital generation of realistic 3D foam morphologies, based on a novel combination of tools such as sphere packings, Voronoi-Laguerre diagrams and Surface Evolver. Overall our framework allows finer control of the morphological parameters than others previously seen in the literature, and the generation of more realistic structures. For open cell high porosity structures, we validated our results by comparing cell size distribution and cell connectivity distribution with experimental tomographic data. In the future, it will be interesting to extend the experimental validation to all the morphologies that the method can produce. Additionally, the generation capabilities can be made useful in other fields, distinct from radiation heat transfer, where a good representation of microscopic morphology is necessary. Applications to thermal conduction and solid mechanics are already underway, and we hope to further enlarge the scope of application.

Concerning radiation heat transfer, we've been mainly focused on presenting analyses and methods that maximize simplicity while not sacrificing effectiveness.

For cellular media with an opaque solid phase, significant effort was devoted to determination of the most appropriate methodologies among those available in literature and to the development of simplified analytical relations. Using increasingly accurate morphological models, we've been able to propose more accurate relations that make it possible to calculate the extinction coefficient directly from the knowledge of easily measurable morphological parameters, with average deviation below 2.5% when compared to direct tomographic analysis.

For cellular media with a semi-transparent solid phase, the sparsity and inconsistency of literature comparing homogenized methods with reference solutions pushed us to propose a new generalized Direct Monte Carlo Homogenization (DMCH) reference method with vastly improved characteristics in terms of memory occupation and parallelizability, applicable to any kind of reference problem/geometry. It is of particular interest that the computational costs of our direct method are somewhat comparable to those of the Homogeneous Phase method, with computational times of the order of minutes. With the increasing availability computational power and massive CPU/GPU parallelization, we can imagine such direct methods becoming more and more useful not just as reference methods but also as standalone tools.

Subsequent analysis using the newly developed reference method has revealed significant problems with existing standard Homogeneous Phase (HPA) and Multi Phase (MPA) methods in semitransparent media, with errors in excess of 30% for very simple problems, and pushed us to develop our own Improved Homogeneous Phase (HPA+) and Improved Multi-Phase (MPA+) methods. Once again, with maximum simplicity as one of our targets, the improved methods are based on the systematic use of hybrid direct-inverse parameter identification techniques and on minor alterations of the existing homogenized equations, by the addition of a "trapped" non-interacting phase and boundary condition adjustment. The form of the homogenized equations stays fairly recognizable, and they can be solved with the entire array of techniques already known and used for the standard HPA/MPA. In spite of their simplicity, we have shown that the Improved methods allow to realize a very significant reduction of error and to consistently achieve acceptable errors (<10%) over a vast range of cellular morphologies and physical configurations, being very competitive even when compared with significantly more complex homogenized methods from the most recent literature. This work also opens up some possibilities for future evolutions and applications:

-The morphology generation algorithms and the radiative models could be modified to make it possible to take into account anisotropy. -The development of the DMCH as a general purpose radiation simulation tool and its integration with existing radiation/heat transfer simulation toolchains could be pursued.

-The HPA+ and MPA+ models could be applied and validated on non-cellular porous media with a semi-transparent phase, e.g. fibrous media.

FOLIO

Fig. 1 -Fig. 3 -

 13 Fig. 1 -Mousse polymère à pores fermés Fig. 2 -Mousse de carbone à pores ouverts

Fig. 5 -

 5 Fig. 5 -Isolation d'un bâtiment avec mousse polymère

Fig. 9 -

 9 Fig. 9 -L'approche multi-échelle. De gauche à droite: (A) Le milieu et sa microstructure. (B) Le modèle hétérogène du milieu (C) Le modelé homogène équivalent du milieu. Les flèches en (B) et (C) représentent les flux locaux des quantités physiques (par ex. l'énergie), qui sont fortement discontinus dans le milieu hétérogène mais continus dans le milieu homogène équivalent.

Fig. 10 -

 10 Fig. 10 -Résultats numériques pour le transfert radiatif de la section 1.Enfin, le cas plus complexe d'une phase solide semi-transparente induisant des interactions

Fig. 11 -Fig. 12 -Fig. 13 -

 111213 Fig. 11 -Mousse à pores ouverts à faible porosité (ε = 70%)

Fig. 14 -Fig. 15 -

 1415 Fig. 14 -Distributions du diamètre équivalent des cellules pour les 4 échantillons et pour 4 structures générées. L'aire en commun est remplie.

Fig. 17 -

 17 Fig. 17 -Lancement des rayons dans la méthode d'Homogénéisation Monte Carlo Directe (DMCH)

-Fig. 18 -

 18 Fig. 18 -Example of inverse coefficient fitting.

Fig. 20 -Fig. 21 -Fig. 22 -

 202122 Fig. 20 -Transmittance/Réflectance d'un échantillon.

3. 3 4 Radiative

 34 Open cell porous structures .............................................................................................. 3.4 Voronoi-Laguerre diagram and Surface Evolver processing ........................................... 3.5 Closed cell plastic foams .................................................................................................. 3.6 Open cell metal foams ...................................................................................................... Chapter Properties of Irregular Open Cell Solid Foams 4.1 Introduction ...................................................................................................................... 4.2 Methodology ...................................................................................................................

  HPA+ Nv Number of scatterers per unit volume (m -3 ) t Trapped phase of HPA+/MPA+ Ns Number of tomography slices ij From phase i ∈ {1, 2} to phase j ∈ {1, 2}, Standard deviation of the distribution of x var(x)

  0.1), carbon (Fig. 0.2), ceramic (Fig. 0.3) and metal foams (Fig. 0.4) constitute a large and varied class of modern engineering materials. The peculiar characteristics they derive from their cellular structure make them interesting for a number of applications in the field of thermal engineering.

Fig. 0. 1 -

 1 Fig. 0.1 -Closed cell plastic foam Fig. 0.2 -Open cell carbon foam

Fig. 0. 5 -

 5 Fig. 0.5 -Plastic foam insulation of buildings.

Fig. 0. 6 -

 6 Fig. 0.6 -Foam heat exchanger Fig. 0.7 -Porous burner Fig. 0.8 -Jülich solar tower [11]

Fig. 0. 9 -

 9 Fig. 0.9 -The multi scale approach. From left to right: (A) The real medium and its microstructure. (B) The heterogeneous medium model (C) The equivalent homogeneous medium model. The arrows in (B) and (C) represent the local flows of physical quantities (e.g. energy), which are highly discontinuous in the heterogeneous medium but continuous in the equivalent homogenous medium.Over the course of the years, a number of numerical methods to determine effective radiative

  commonly occurring types of foam structures, covering a large range of real materials: low-porosity open-cell structures (Fig. 0.11), high-porosity closed cell foams (Fig. 0.12), high-porosity open cell foams (Fig. 0.13).

Fig. 0. 11 -

 11 Fig. 0.11 -Open cell porous structure (ε=70%) Fig. 0.12 -High porosity closed cell foam (ε=85%)

  Efficient and flexible reference methods are required to validate the results obtained by homogenization models. A purely numerical, Direct Monte Carlo Homogenization reference technique, based on a periodic REV, allows to simulate radiation within arbitrary cavities with arbitrary boundary conditions and calculate macroscopic radiative quantities (such as transmittance, reflectance, configuration factors etc.) using ray-counting methods typically applied in MCRT for participating media. The technique is validated against full scale Monte Carlo simulations and compared to the existing HPA model.Recent models have been proposed in literature to increase the accuracy of homogenized methods, but they come with significant increase in complexity. On the basis of numerical and phenomenological observation of the sources of discrepancy between homogenized models and underlying discrete physical phenomena, we propose improvements of the existing Homogeneous Phase Approach and Multi Phase Approach, specifically targeted at the case of foams with a semitransparent solid phase, with the objective to increase accuracy while introducing as little additional complexity as possible. The two main ideas driving the modeling effort are the use of hybrid directinverse coefficient determination approaches and the simplified modeling of ray history effects. The resulting Improved Homogenized Approaches (HPA+ and MPA+) are extensively tested by comparing them with Direct Monte Carlo Homogenization simulations and existing homogenized models, on a varied set of morphologies making full use of the previously developed digital generation techniques, comprising closed cell foams and low and high porosity open cell foams, in 3 typical radiation heat transfer configurations.

Fig. 1 . 1 -

 11 Fig. 1.1 -The multi scale approach. From left to right: (A) The real medium and its microstructure. (B) The heterogeneous medium model (C) The equivalent homogeneous medium model. The arrows in (B) and (C) represent the local flows of physical quantities (e.g. energy), which are highly discontinuous in the heterogeneous medium but continuous in the equivalent homogenous medium.

[ 49 ]

 49 [START_REF] Baillis | Measurements of radiation characteristics of fused quartz containing bubbles[END_REF][51][52][53]. In recent years, numerical methods based on Monte Carlo techniques for the determination of radiative properties are becoming established in order to study either real structures obtained from tomographic imaging or computer-generated structures that closely mimic the microstructure of the real foam. Monte Carlo methods are based on direct simulation of the propagation of a large number of photons (generally at least 10 6

1 . 1 )

 11 As explained in the introduction, most of the research focuses on the determination of the relevant parameters of the RTE. These coefficients are the extinction coefficient 𝛽, which represent the ability of the media to interact in with radiation by means of scattering or absorption, the scattering albedo 𝜔 = 𝜎 𝛽 , which represents the prevalence of scattering events, and the scattering phase function Φ(Ω ′ , Ω), which captures the patterns of scattering itself. It is worthwhile to focus our attention now on the specific assumptions used for the calculation of the radiative properties (extinction coefficient, absorption coefficient and phase function) required in the RTE. If the size parameter 𝜒 = 𝜋𝐿 𝜆

Fig. 1 . 2 -

 12 Fig. 1.2 -Ray casting methods: (A) Casting from void region. (B) Casting from solid/fluid interface

Fig. 1 . 3 -

 13 Fig. 1.3 -Boundary handling methods -(A) Symmetrical reflection; (B) Periodic wrapping with position randomization (C) Recasting from random point in void region.

Fig. 1 . 4 -

 14 Fig. 1.4 -Typical IOOS geometry -(a) Overlapping spheres; (b) Spheres cut to fit bounding box.

Fig. 1 . 5 -

 15 Fig. 1.5 -Extinction coefficient as a function of porosity: numerical methods and analytical results.

Fig. 1 . 6 -

 16 Fig. 1.6 -Average error of numerical methods as a function of number of rays cast.

Fig. 1 . 7 -

 17 Fig. 1.7 -Average error of Monte Carlo methods as a function of number of scatterers for different boundary handling. The results can be observed in Fig. 1.7. A number of 10 6 rays, cast from the solid-fluid interface, was

Fig. 1 . 8 -

 18 Fig. 1.8 -Error of Monte Carlo methods as a function of porosity with or without boundary handling.

Fig. 1 . 9 -

 19 Fig. 1.9 -Average error of Monte Carlo as a function of number of scatterers for different ray origins.

Fig. 1 . 10 -

 110 Fig. 1.10 -Average error of various numerical methods as a function of number of scatterers.

  Fig. 1.11 for a circular cross section with  = 1 and ε = 93.73%.

Fig. 1 . 11 -

 111 Fig. 1.11 -Voronoi geometry. (a) Geometry cut to the bounding box (b) Detail of the junctions.

  value of the perturbation coefficient, each row reports the struts' cross sections, the porosity and Sv ratios for a single geometry. The values  = 0.5 and  = 1 have been chosen for the perturbation coefficient since they lead to cell size distributions with standard deviations close to real cellular geometries commonly encountered for open cell metal foams[START_REF] Randrianalisoa | Microscale direct calculation of solid phase conductivity of voronoi s foams[END_REF].

Fig. 1 . 13 -

 113 Fig. 1.13 -Extinction coefficient for Voronoi structures with triangular and circular strut as a function of porosity for  = 0.5 and  = 1.0.

  have been found to be nearly insensitive to choice of ray origination, with some marginal effects on convergence speed, which favour the choice of rays cast form the surface. The effects of boundary handling in Monte Carlo methods have been analysed. The choice of boundary handling does not appear to significantly affect final result, but can affect convergence speed. In particular, periodic wrap-around of rays with position randomization appears to improve convergence speed over all other methods. It has also been shown that such a method is vastly superior to the casting of rays from the core and discarding of exiting rays, in cases where sufficiently high optical thickness cannot be guaranteed.Monte Carlo methods have been additionally compared to projection based methods. The results

Fig. 2 .

 2 Fig.2.1 gives a general scheme of a direct methods, such as those seen in[15][22][32][36][37]

Fig. 2 . 1 -

 21 Fig. 2.1 -Flow diagram scheme of direct identification methods for radiative properties

Fig. 2 . 2 -

 22 Fig. 2.2 -Flow diagram of inverse identification methods for radiative properties

Fig. 2 . 3 -

 23 Fig. 2.3 -Flow diagram of hybrid direct-inverse identification method for radiative properties

Fig. 2 . 4 -

 24 Fig. 2.4 -Variables in Snell's law and Fresnel equations.

Fig. 2 . 5 -

 25 Fig. 2.5 -Microscale phenomena in the Multi Phase Approach and corresponding homogenized coefficients.

Fig. 3 . 1 -

 31 Fig. 3.1 -Open cell porous structure (ε = 70%)

Fig. 3 . 5 -

 35 Fig. 3.5 -Scheme of the inflation-packing algorithm

Fig. 3 .

 3 Fig. 3.6 -(A) Explicitly represented domain (unwrapped spheres). (B) Virtual domain used for contact calculation with periodic wrapping. (C) Domain cut to periodic bounding box. Note that both (A) and (C) are periodic units that can be used to tessellate space.

Fig. 3 .

 3 Fig. 3.7 presents a graphic overview of the typical course of the algorithm. The algorithm allows to

Fig. 3 . 8 -

 38 Fig. 3.8 -Packed spheres (left) and inflated intersecting spheres (right).

Fig. 3 .Fig. 3 . 9 -

 339 Fig. 3.9 -The three steps of porosity fine-tuning process from left to right. (A) Original voxels (B) After Gaussian Blur (C) After thresholding. The final structure is more porous, while preserving general topology.Blur is exaggerated for clarity.

Fig. 3 .

 3 Fig. 3.11 -Voronoi-Laguerre tessellation in 2D (left) and 3D (right) [46]. Included spheres in pink, cell vertices in dark blue. For the sake of visibility, no cell faces and only a part of cell vertices has been represented in the 3D picture.

Fig. 3 .

 3 Fig. 3.12 -(Left) Initial cell structure (Voronoi-Laguerre diagram). (Right) Processed cell structure (after Surface Evolver). For display purposes, the periodic structure is unwrapped and each cell is displayed only one time.The use of Voronoi diagrams and Surface Evolver to generate realistic closed cell structures has been

Fig. 3 .

 3 Fig. 3.13 -Original (left) and exploded (right) cell structure. Colors identify cells. Walls are pattern filled.

Fig. 3 .

 3 Fig. 3.14 -Detail of explosion of one pair of cells, with indication of magnitudes involved in calculation.

  vectors 𝐹𝐷 𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗ being the independent variables. The displaced resulting cells are shown in Fig. 3.15. The displaced cells are then used to carve holes in a cubic bounding box thus resulting in the finalized constant thickness closed cell structure (Fig. 3.16). The porosity of the final structure can be controlled by varying the target thickness of the walls THICK. The final output is a periodic structure of controlled cell distribution, porosity and mesh refinement. Compared with existing literature [78][79][80], the main novelty of this method is the full 3D explicit representation of the structure, which has been usually eschewed in previous efforts in favor of 2D shell representations. This should allow to improve accuracy of results, especially at relatively low (i.e. < 95%) levels of porosity.

Fig. 3 .

 3 Fig. 3.15 -Full cell structure (left) and close up (left) showing the space between cells (wall thickness) in 3D.

Fig. 3 .

 3 Fig. 3.16 -Finalized high porosity closed cell structure

  Fig. 3.18 -Cross section of the strut Fig. 3.19 -Longitudinal profile of the strut

Fig. 3 .

 3 Fig. 3.20 -Structure before (a) and after (b) the "virtual welding" process.

Fig. 3 .Chapter 4 -

 34 Fig. 3.21 -Example of finalized structure (CV = 5%; t = 1; ε = 94%; k = 1)

Fig. 4 . 1 -

 41 Fig. 4.1 -Progressive dropping of spheres into space at 5% (a), 16% (b) and 30% (c) packing

Fig. 4 .Fig. 4 . 3 -

 443 Fig. 4.2 -(a) Initial cell structure (Voronoi diagram) (b) Processed cell structure (after Surface

Fig. 4 . 5 -Fig. 4 . 6 -

 4546 Fig. 4.5 -Structure before (a) and after (b) the "virtual welding" process.

Fig. 4 . 7 -

 47 Fig. 4.7 -Equivalent cell diameter distributions of the 4 samples, as calculated by iMorph.

Fig. 4 . 9 -

 49 Fig. 4.9 -Equivalent cell diameter distributions of the 4 samples and respective cell size distribution for matched generated structures, as calculated by iMorph. Shared areas are filled.

Fig. 4 .

 4 Fig. 4.10 -Average cell connectivity distribution of the 4 samples and average cell connectivity distribution for matched generated structures, as calculated by iMorph. Shared area is filled.

Fig. 4 .

 4 Fig. 4.11 -Computed extinction PDFs and exponential fit.

[ 15 ] 3 )

 153 Fig. 4.12 -Computed scattering phase functions and analytical relation (4.3).

= 1 .

 1 0) is generated and then parameters are varied individually to evaluate the effects. The values of parameters have been chosen to span the range of typical values observed in real foams. All the results are computed for an average cell size of 1 mm. It is worth noting that, operating under geometric optics approximation, the average mean free path is directly proportional to the average cell diameter, i.e. β = 1 ⁄ savg ∝ 1 ⁄ dc.

Fig. 4 . 13 -

 413 Fig. 4.13 -Effects of the variations of CV on β and β + Fig. 4.14 -Kelvin foam

Fig. 4 . 15 -

 415 Fig. 4.15 -Effects of the variations of t on β and β + Variation of strut diameter along the length has moderate effects on the extinction coefficient. All

Fig. 4 .

 4 Fig. 4.18 -Comparison of analytical law (4.4) and Monte Carlo results on generated structures

  Fig. 4.19 -Comparison of analytical law (4.5) and Monte Carlo results on generated structures.Dimensionless parameters k and t, necessary for the computation of (4.4) and (4.5), can be easily

Fig. 4 .

 4 Fig. 4.20 -Typical tomography slice for each of the four samples.

Fig. 4A. 1 -

 1 Fig. 4A.1 -Variation of 4Ss/do 2 with k.

Finally, relation (4A. 1 )

 1 was used as reference to attribute values of normalized curvature to the tomographic samples (Fig. 4A.2).

Fig. 4A. 2 -

 2 Fig. 4A.2 -Determination of k for the four samples.

Finally, relation (4B. 1 )

 1 is used as reference to determine the value of the diameter ratio for the four samples (Fig.4B.2).

Fig. 4B. 2 -

 2 Fig. 4B.2 -Determination of t for the four samples.

Fig. 5 . 1 .

 51 Fig. 5.1. Creation of a periodic computational domain from non-periodic input data.

Fig. 5 .

 5 Fig. 5.2 schematizes an example on input data required, showing a RVE (which is used as the actual

Fig. 5 . 2

 52 Fig. 5.2 Input data of the algorithm. (A) Computational domain (REV). (B) Physical domain. The REV is also shown in the physical domain to give an idea of the difference in absolute size.

Fig. 5 . 3

 53 Fig.5.[START_REF] Coquard | Modeling of heat transfer in low-density EPS foams[END_REF] The ray casting process for three rays from three different source points in the computational domain. Periodic transition events are indicated with crosses. The origin points in the computational domain are randomized and the physical reference system is randomized.

Fig. 5 .

 5 Fig. 5.3 schematizes the actual ray casting process for four rays (in black, red, blue and green). In

Fig. 5 . 4 .

 54 Fig. 5.4. Ray trajectory in (A) physical reference system VS (B) computational domain.

Fig. 5 . 6 .

 56 Fig. 5.6. Schematization of physical configuration 1 and relevant magnitudes.

Fig. 5 . 7 -Fig. 5 . 9 -

 5759 Fig. 5.7 -The 98% porosity RVE Fig. 5.8 -The 75% porosity RVE

Fig. 6 . 1 -

 61 Fig. 6.1 -Example geometry to illustrate the effect of ray origin on ray propagation.

a) A trapping

  and absorption effect for a fraction of the energy emitted in the dense phase or entering it from outside (blue fraction in figure). This effect will be captured by defining a trapped fraction C and its relative absorption coefficient αt. b) A scattering effect for all the energy emitted in the dense phase or entering it from outside (blue and green fractions in figure). This effect will be captured by defining an additional scattering coefficient gHG.

Fig. 6 . 2 -

 62 Fig. 6.2 -Synthetic figure depicting the effects on the discrete scale and the corresponding approximations at the homogenized scale. In blue and green, the trapped and non-trapped rays respectively, represented on the left at the discrete scale and on the right as the approximated homogenized scale equivalent. In red the preserved distributions, s and θ. At the right side, the Heyney-Greenstein distribution of θ.

Fig. 6 .

 6 Fig. 6.2 synthetically displays the effects on the discrete scale and the corresponding approximations

Fig. 6 . 3 -

 63 Fig. 6.3 -Parameter identification flow chart for Improved Homogeneous Phase Approach

Fig. 6 . 4 -

 64 Fig. 6.4 -Illustration of ray casting from fluid phase for hybrid determination of radiative properties and logged quantities.

Fig. 6 . 5 -

 65 Fig. 6.5 -Illustration of the inverse fitting process (95% porosity open cell foam).

Fig. 6 . 6 -

 66 Fig. 6.6 -Illustration of ray casting from solid phase for evaluation of αt and ghg

Fig. 6 . 7 -

 67 Fig. 6.7 -Comparison of scattering from Direct Monte Carlo simulation (top), HPA-MC simulation (bottom left) and HPA-MC simulation (bottom right) augmented with Heyney-Greenstein scattering at the moment of entry of the rays into the domain.. The short range effects associated to the passage from solid to fluid give raise to a stronger diffusion in the Direct Monte Carlo. The differences are especially remarkable at short distances. The directional distribution transformation approximates well the short range effects.

Fig. 6 . 8 -

 68 Fig. 6.8 -Comparison of numerically calculated distribution P* and respective Heyney-Greenstein function ΦHG for a high-porosity open cell foam (ε= 0.89) One may observe that would be possible to straightforwardly determine C and αt from their phenomenological definitions by counting the proportion of absorbed rays (i.e. 𝐶 = 𝑁 𝑠 𝑁 ⁄ ) and their mean length (i.e. 𝛼 𝑡 = 1 ‖𝑠 ‖ ̅̅̅̅̅ ⁄). In the next section, it will be shown that it is possible to obtain more

Fig. 6 .

 6 Fig. 6.10 -Illustration of the prescribed intensity open boundary condition. Note that the treatment of radiation entering the solid phase differs from radiation entering the fluid phase, but treatment for radiation exiting the domain is the same in all cases.The first boundary condition discussed is a prescribed incident directional intensity Iw[43][44][START_REF] Brewster | Thermal radiative transfer and properties[END_REF]in

(6. 13 )Fig. 6 . 11 -

 13611 Fig. 6.11 -Illustration of the diffusely emitting and reflecting wall boundary condition. The treatment of radiation entering the solid phase differs from radiation entering the fluid phase.

  these are chosen to represent typical values of semi-transparent materials used in porous structures in the visible to the intermediate infrared region, such as semi-transparent plastics and ceramics[98][99][100][101][102][103]. The region of values where the behavior of the solid is more markedly semi-transparent is chosen for analysis (as opposed to "almost transparent" or "almost opaque" behavior). Nine different triplets of (α2, L, n2) are considered for each physical configuration / morphology combination (for a total of 135 tests for each resolution method).As a benchmark, the radiation problems are solved using the Direct Monte Carlo Homogenization method illustrated in Chapter 5. Using this benchmark as a base, three different homogenized models are compared: 1) HPA: a standard HPA model[22][40], using Eq. (1.1) with standard BCs for open boundaries and opaque walls (Eqs. 6.10, 6.13), and determination of properties according to the method shown in[40], i.e. determination of extinction coefficient β and scattering albedo ω with Eqs. (6.1-2) and direct calculation of scattering phase function by counting[40]:(6.28)

Fig. 6 . 13 -

 613 Fig. 6.13 -Schematization of physical configuration 1 and relevant magnitudes.

2 )

 2 Simulation of Lambertian emission for an eccentric surface element in a cubical domain of side L and calculation of view factors Cfi (i = 1…6) with respect to the domain's walls and total absorption Abs (Fig. 6.14).

Fig. 6 .

 6 Fig. 6.14 -Schematization of physical configuration 2 and relevant magnitudes. Note the eccentric position of the source point. The quantities are directly evaluated by counting the rays [43][44][104], i.e.: 𝐶𝑓 𝑛 = 𝑁 𝑛 𝑁 ; 𝐴𝑏𝑠 = 𝑁 𝐴 𝑁

  6.14) and NA is the number of rays absorbed in the domain. It is worth noting that ∑ 𝑁 𝑖 6 𝑖=1 + 𝑁 𝐴 = 𝑁. The error is quantified as: 𝑬𝑹𝑹𝑶𝑹 = ∑|𝑪𝑭 𝒏 -𝒄𝒇 𝒏 | 𝟔 𝒏=𝟏 + |𝑨𝑩𝑺 -𝒂𝒃𝒔| (6.32)

Fig. 6 .

 6 Fig. 6.15 -Schematization of physical configuration 3 and relevant magnitudes. Note that the temperature only depends on linear position.

Fig. 7 . 1 -

 71 Fig. 7.1 -Example geometry to illustrate the effect of ray origin on ray propagation.

  Fig. 7.2 -Two different path length distributions with the same average path length. The distribution on the right approximates the same path length as the distribution on the left but only use borders of the solid domain as starting point.

Fig. 7 .

 7 Fig. 7.3 synthetically displays the effects at the discrete scale and the corresponding approximations

Fig. 7 . 3 -

 73 Fig. 7.3 -Synthetic figure depicting the effects on the discrete scale and the corresponding approximations at the homogenized scale.In blue and green, the trapped and non-trapped rays respectively, represented on the left at the discrete scale and on the right as the approximated homogenized scale equivalent. In black the rays in the fluid phase. The homogenized trapped phase scatters and absorbs. For the non-trapped rays, the mean path length is preserved.

Fig. 7 . 4 -

 74 Fig. 7.4 -Parameter identification flow chart for Improved Multi Phase Approach

7. 3 . 1 .

 31 Ray tracing from the fluid phasehybrid determination of coefficients.

Fig. 7 . 5 -

 75 Fig. 7.5 -Illustration of ray casting from fluid phase for hybrid determination of radiative properties and

  𝑀𝑃𝐴 𝑔 𝑖𝑗 (𝑠) indicates the locally absorbed fraction calculated using the MPA with scattering certain set of asymmetry factors gij and 𝐼 𝑎𝑏𝑠 * (𝑠) indicates the locally absorbed fraction calculated for the micro scale direct Monte Carlo.

Fig. 7 . 6 - 7 . 3 . 2 .Fig. 7 . 7 -

 7673277 Fig. 7.6 -Illustration of the fitting process. In blue, the direct numerical polar distribution of absorbed radiation. In red, yellow and green, the error surfaces (

7. 3 . 3 .Fig. 7 . 8 -

 3378 Fig. 7.8 -Illustration of ray casting from solid phase for evaluation of trapped phase coefficients.

Fig. 7 . 10 -

 710 Fig. 7.10 -Illustration of the physical phenomena involved in the prescribed intensity open boundary condition. Note that the treatment of rays incident on the solid differs from rays incident on the fluid phase.With reference to Fig.7.10, the first boundary condition discussed is a prescribed incident directional

Fig. 7 .

 7 Fig. 7.11 -Illustration of the physical phenomena involved in the diffusely emitting and reflecting boundary condition. The treatment of rays incident on the solid phase differs from rays incident on the fluid phase. Note: the angular distribution of intensity is the same for the three phases for the radiation emitted and reflected in phase 2, the differently colored arrows represent the partitioning of intensity across the 3 phases. With reference to Fig. 7.11, the second boundary condition discussed is a diffusely emitting and

  (2.4a-b) with standard BCs for open boundaries and opaque walls (Eqs. 7.23, 7.25) and determination of properties according to the methods shown in[22][40]. It must be noted that the two references cited share the characteristic of considering a single scattering event for each ray, but differ significantly in terms of choice of initial points and directions: in some cases the results are quite different.

4 ) 5 )

 45 MPA-GRTE: a Multi-Phase Approach based on Generalized Radiative Transfer Equation (Eqs. 2.15a-b) and tracking of ray histories up to 3 levels of depth, based on the methods shown in [42]. MPA+: an extended multi-phase model, using Eqs. (7.14a-d) and modified BCs for open boundaries (Eqs. 7.24a-c) and opaque walls (Eqs. 7.26a-c, 7.27), with properties determined following the entire process discussed in Paragraph 7.3. Evaluation of deviation with respect to the Direct Monte Carlo simulation is conducted according to the methods described in Paragraph 6.4.

  Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2018LYSEI005/these.pdf © [S. Cunsolo], [2018], INSA Lyon, tous droits réservés -The digitally generated morphologies could be used to model other phenomena, such as fluid flow, convection or chemical reactions.

  and multi-RTE approaches. Journal of Quantitative Spectroscopy and Radiative Transfer, 202, 220-232. [44] Zhang, H. P., & Makse, H. A. (2005). Jamming transition in emulsions and granular materials. Physical Review E, 72(1), 011301. [45] Fang, Q., & Boas, D. A. (2009, June). Tetrahedral mesh generation from volumetric binary and grayscale images. In Biomedical Imaging: From Nano to Macro, 2009. ISBI'09. IEEE International Symposium on (pp. 1142-1145). IEEE. [46] Rycroft, C. (2009). Voro++: A three-dimensional Voronoi cell library in C++. Lawrence Berkeley National Laboratory. [47] Brakke, K. A.

Chapter 1 Radiative properties modeling of open cell solid opaque foam

  The resulting Improved Homogenized Approaches are extensively tested by comparing them with Direct Monte Carlo Homogenization simulations and existing homogenized models, on a varied set of morphologies making full use of the previously developed digital generation techniques. The improved models consistently outperform existing homogenized models. Introduction ........................................................................................................................................ 8 Results and comments ................................................................................................ 31 1.4 Conclusions ...................................................................................................................... 33

1.1 Introduction ...................................................................................................................... 10 1.2 Radiative properties modeling ......................................................................................... 13

1.2.1 Analytical method .............................................................................................. 14 1.2.2 Projection method ............................................................................................. 16 1.2.1 Monte Carlo method .......................................................................................... 17

1.2.3.1 Variations of ray casting method -Ray origin ........................................................... 17 1.2.3.2 Variations of ray casting method -Boundary handling.............................................. 18 1.2.3.3 Calculation of radiative properties .............................................................................. 20 1.3 Numerical results and comments ..................................................................................... 21 1.3.1 Spheres .............................................................................................................. 21 1.3.1.1 Geometry description.................................................................................................. 21 1.3.1.2 Results and comments ................................................................................................ 23 1.3.1 Voronoi structures ............................................................................................. 28 1.3.1.1 Geometry description.................................................................................................. 28 1.3.1.2

Chapter 2 Indirect identification methods and semi-transparent cellular media 2

  .1 Classification of methods for the identification of radiative properties ........................... 2.2 Modeling radiation in cellular media with a semi-transparent solid phase ...................... 2.2.1 Notes on the Multi Phase Approach (MPA) ...................................................... 2.2.2 Other recent developments -GRTE .................................................................. 2.2.3 General observations ........................................................................................

	SECTION 2
	DIGITAL GENERATION OF REALISTIC CELLULAR MORPHOLOGIES
	Chapter 3
	A

general framework for the generation of realistic cellular morphologies

3.1 Introduction ...................................................................................................................... 3.2 Close random sphere packing generation.........................................................................

  .5 Conclusions ...................................................................................................................... 4.A Evaluating the normalized curvature (k) .........................................................................

		Nomenclature
	Latin symbols	Greek symbols
	B 0	Blackbody intensity (W/m•sr)	α	Absorption coefficient (m -1 )
	C	Trapped fraction coefficient, HPA+	β	Extinction coefficient
	C1	Trapped fraction coefficient, MPA+	βI	Extinction coefficient -Least squares fitting (m -1 )
	C2	Mean path length adjustment coefficient, MPA+	βII	Extinction coefficient -Inverse average length (m -1 )
	CV	Coefficient of variation	βc	extinction coefficient -Corrected (units -1 )
	d	Cell / sphere diameter	βan	Extinction coefficient -Analytical (units -1 )
	ds	Strut diameter	βnum	Extinction coefficient -Numerical (units -1 )
	f	Volume fraction	β +	Extinction coefficient -nondimensional
	fc(ε)	Extinction coefficient correction function	ε	Porosity
	g	Phase function asymmetry factor	εs	Solid surface emissivity
	G	Area projected by a single scatterer	λ	Wavelength (m)
	Ge(s)	Extinction distribution function	μ	Scattering angle cosine
	GCV	Geometric coefficient of variation	μn	Scattering angle cosine of the n-th realization
	I	Radiation intensity (W/m•sr)	Φ(Ω,Ω')	Scattering phase function
	I(s)	Intensity at distance s from source (W/m•sr)	Φ(θ) Φ(μ)	θ = cos -1 (Ω•Ω') μ = Ω•Ω'
	I(r, Ω)	Intensity at position r in direction Ω (W/m•sr)	θ ρs	Angle between two directions Solid surface reflectivity
	I0	Intensity at source (W/m•sr)	σ	Scattering coefficient (m -1 )
	Iabs(s)	Intensity absorbed locally at distance s (W/m•sr)	χ ω	Size parameter Scattering albedo
	Iabs(s, θ)	{s, θ} (W/m•sr) Intensity absorbed locally at polar location	Ω	Direction vector
	k	Normalized curvature of strut cross section

4.2.1 Digital generation methodology ........................................................................ 4.2.1.1 Digital generation methodology -Presentation .......................................................... 4.2.1.2 Digital generation methodology -Validation ............................................................. 4.2.2 Radiative properties calculation methodology .................................................. 4.3 Results and discussion ..................................................................................................... 4.4 Simplified analytical relations ......................................................................................... 4.4.1 Analytical relations -Presentation .................................................................... 4.4.2 Analytical relations -Validation .......................................................................

44.B Evaluating the diameter ratio (t) .....................................................................................

ε Radius [units] S/V [1 / units]

	99.5%	8.44	0.001784
	99%	10.74	0.002852
	98%	13.80	0.004651
	95%	18.45	0.008365
	90%	23.32	0.01357
	85%	26.47	0.01765
	Table

[units -1 ] ε

	5.0E-03
	4.5E-03
	4.0E-03
	3.5E-03
	3.0E-03
	2.5E-03
	2.0E-03
	5.0E-04 1.0E-03 1.5E-03 84% β Projection method 86% 88%	90%	92% Monte Carlo 94%	96% Analytic solution 98%	100%	102%

Table 1 .

 1 2 -Summary of the characteristics of the Voronoi geometries.

	94.04%	0.01137	93.73%	0.01149
	92.12%	0.01273	91.69%	0.01273
	90.02%	0.01388	89.45%	0.01398
	85.08%	0.01607	84.41%	0.01596
	97.55%	0.01031	97.41%	0.01052
	94.69%	0.01486	94.36%	0.01495
	90.83%	0.01860	90.34%	0.01858
	88.12%	0.02042	87.66%	0.02043

Table 1 .

 1

2 reports the range of relevant geometric parameters as follows: the geometries are grouped into two columns according to the value of perturbation coefficient. In any single column,

  As already discussed in the General Introduction and recalled in Section 1, the current state of the art Monte Carlo techniques for the determination of radiative properties of porous media depend on the accurate knowledge of the porous morphology at the micro-scale. The object of analysis is the so

	reproduce	the	foam	structures	using	different	approaches

3.1. Introduction called Representative Volume Element (RVE): a limited amount of porous material fully representing the morphological characteristics of the material at a microscopic scale, thus allowing to characterize it. Techniques fully based on tomographic data provide satisfactory agreement with experimental data, but their dependence on high quality scans of existing foam samples makes them of limited utility for design purposes. To overcome these limitations, a number of recent studies have sought to digitally [15][29][30][31][32][33][34][35][71][75][76][77][78][79][80].

respect to the methodology as a whole, it is worthwhile to remark that, compared to methods existing in literature, the peculiarity of the present approach is that it allows to generate a number of morphologically diverse porous cellular structures with high control of the main morphological parameters in the context of a unified coherent process. Another significant characteristic common to the entire method is

			Start		Drops spheres from largest
						to smallest avoiding
			Relative	Number of cells		overlaps. 30% packing fraction.
			diameter
			distribution
			Drop spheres in bounding		Packing fraction of at
			box by RSA algorithm		least 62%
					2
			Packs spheres by inflation,
			contact force calculation and
	Each sphere is the		random shaking		Porosity
	included sphere of its
	Voronoi-Laguerre cell				4	Inflate spheres	3
			Calculate Voronoi-		allowing
			Laguerre tessellation		intersection
	Standard evolution +		with voro++ [46]
	Annealing [78]					Convert geometry
			Stabilize cellular		to voxels
	Mesh size	5	structure with Surface Evolver [47]		Fine-tune porosity
	of triangles Refine structure up to desired number		Extract skeleton	6	by Gaussian blurring and thresholding
	Separate cells		Porosity	Strut shape		Invert voxels
	(duplicate walls)			Mesh size		Mesh size
	Porosity	Add 3D struts
						Mesh the voxel
	Displace cells to					structure with
	obtain uniform		Weld struts together		iso2mesh [81]
	wall thickness		at the junctions
						Cut to
	Cut to bounding box		Cut to bounding box		bounding box
	End		End			End
	High porosity closed cell foam	High porosity open cell foam		Open cell porous structure
	Legend	Notes	Endpoints	Processes	Inputs	Subsections 2-6

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2018LYSEI005/these.pdf © [S. Cunsolo], [2018], INSA Lyon, tous droits réservés grouped by number (2-6) according to Fig. 3.4. With

Table 4

 4 .1. Samples 1 and 2 have been used in previous work of the authors, while samples 3 and 4 are new experimental material, introduced to have a better coverage of the typical range of porosities of metal foams.

	Sample number	Material	Nominal porosity	Nominal PPI	μCT resolution
	1 [32]	Al-NiP	90%	60	30 μm
	2 [7]	Al	90%	40	44 μm
	3	Al	94%	40	22 μm
	4	Al	97%	40	22 μm

Table 4 .

 4

1 -Characteristics of the samples considered.

Table 4 .

 4

	2. For further reference, Fig. 4.20 also shows a typical tomography slice for each of the four
	samples.
	Sample number	ε	dc [um]	Sv [m -1 ]	k (Appendix A)	t (Appendix B)
	1	87.2%	1842	1596	0.42	0.52
	2	89.3%	2431	1066	0.57	0.4
	3	93.6%	2892	742	0.36	0.34
	4	96.3%	2725	756	0.16	0.53

Table 4 .

 4 2 -Measured morphological characteristics of the four samples.

Table 4 . 3

 43

		Relative error
	5)	(4.5)

-Comparison of direct Monte Carlo computations on the four samples, and relations (4.4) and (4.5). Morphological data from Tab. 2.

Table 5

 5

	.1.

Table 5 .

 5 1 -Error of Direct Monte Carlo Homogenization and Homogeneous Phase approach with respect to Full Domain Monte Carlo simulations.

Table 5

 5

	.2:

Table 5 .

 5 2 -Computation time of different methods in various configurations.

Table 6 .

 6 .1: 1 -Different treatments applied to different fractions of radiation.

			Homogenized phase	Trapped phase	Additional
					scattering
	Fluid phase [f1]	X
	Solid phase [f2]	Trapped [f2 * C] Non trapped [f2 * (1 -C)]	X	X	X X

  Considering that αh is typically a fraction of α2, one can surmise that the sensitivity dC/dgHG is significantly smaller than the reciprocal dgHG/dC for admissible values of the coefficients. Both values are subject to numerical incertitude due to the nature of MC simulation and the calculation method itself, which results in small deviations from Eq. (6.23). However, the difference in sensivity implies that, if gHG is adjusted to restore the exact equivalence, large adjustment may be necessary, possibly causing gHG to fall outside of the physically meaningful interval [-1;1].

	Ω • ∇I ℎ = [𝑓 1 𝑛 1 2 𝛼 1 + 𝑓 2 𝑛 2 2 𝛼 2 (1 -𝐶)]𝐵 0 -𝛽 ℎ 𝐼 ℎ + Ω • ∇I 𝑡 = 𝑓 2 𝑛 2 2 𝛼 2 𝐵 𝑣 𝐶 -𝛼 𝑡 𝐼 𝑡 𝜎 ℎ 4𝜋	∫ I ℎ (Ω′)Φ ℎ (Ω ′ , Ω)𝑑Ω′ 4𝜋	(6.9a) (6.9b)
					I = I ℎ + I 𝑡	(6.9c)
	-Prescribed directional intensity at boundary:
	𝐼(𝑟 𝑤 + , Ω) = 𝐼 𝑤 𝑓 1 [Ω = Ω 𝑤 ] + 𝐼 𝑤 𝑓 2 (1 -𝐶) 𝐼(𝑟 𝑤 + , Ω) = 𝐼 𝑤 𝑓 2 𝐶 Φ 𝐻𝐺 (Ω 𝑤 , Ω) 4𝜋 𝐼(𝑟 𝑤 -, Ω) ⏞ = 𝐼(𝑟 𝑤 + , Ω) ⏞ + 𝑓 2 𝐼 𝑤 ⏞ Φ 𝐻𝐺 (Ω 𝑤 , Ω) Φ 𝐻𝐺 (Ω 𝑤 , Ω) 4𝜋 {𝑛 𝑤 • Ω > 0} {𝑛 𝑤 • Ω > 0} 4𝜋 {𝑛 𝑤 • Ω < 0}	(6.11a) (6.11b) (6.12)
	-Diffusely emitting/reflecting wall boundary:
	I ℎ (𝑟 𝑤 + , Ω) = 𝑓 1 (𝜀 𝑠 𝑛 1 2 𝐵 0 +	𝜌 𝑠 𝜋	𝑛 𝑤 •Ω ′ <0 ∫ 𝐼(𝑟 𝑤 , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′	) + 𝑓 2 (1 -𝐶)𝐼 𝑠 (𝑟 𝑤 , Ω) {𝑛 𝑤 • Ω > 0} (6.14a)
		I 𝑡 (𝑟 𝑤 + , Ω) = 𝑓 2 𝐶𝐼 𝑠 (𝑟 𝑤 , Ω) {𝑛 𝑤 • Ω > 0}	(6.14b)
	𝐼 𝑠 (𝑟 𝑤 + , Ω) = (𝜀 𝑠 𝑛 2 2 𝐵 0 + × 1 -	2(𝛼 ℎ 4𝜋 ∫ 𝑔 𝐻𝐺 + 1 2 ∫ Φ 𝐻𝐺 (Ω ′ , Ω′′)|𝑛 𝑤 • Ω ′′ |𝑑Ω ′ 𝑑Ω ′′ -𝛼 2 ) 𝑛 𝑤 •Ω ′ >0 𝑛 𝑤 •Ω ′′ <0 4𝜋 2 ∫ 𝜌 𝑠 𝑛 𝑤 •Ω ′ >0 Φ 𝐻𝐺 (Ω ′ , Ω)𝑑Ω ′ 1 𝑛 𝑤 •Ω ′ <0 𝜋 ∫ 𝐼(𝑟 𝑤 , Ω ′ )|𝑛 𝑓 • Ω ′ |𝑑Ω ′ ) 𝜌 𝑠	(6.25) (6.15)
	-Compatibility conditions:
					𝛼 2 𝛼 𝑡	=	2 𝑔 𝐻𝐺 + 1	(6.20)
	𝑓 1 𝑛 1 2 𝛼 1 + 𝑓 2 𝑛 2 2 𝛼 2 (1 -𝐶) = 𝛼 ℎ (𝑓 1 𝑛 1 2 + 𝑓 2 𝑛 2 2 (1 -𝐶)	𝑔 𝐻𝐺 + 1 2 )	(6.24)
	We set C and αt to fully enforce Eqs. (6.20) and (6.24):
		𝐶 = 1 -	2 (𝛼 ℎ -𝛼 1 )𝑓 1 𝑛 1 𝑓 2 𝑛 2 2 (𝛼 2 -𝛼 ℎ 𝑔 𝐻𝐺 + 1 2 )	(6..26)
					𝛼 𝑡 =	2𝛼 2 𝑔 𝐻𝐺 + 1	(6.27)

  [START_REF] Modest | Radiative heat transfer[END_REF][95][START_REF] Cengel | Heat transfer a practical approach[END_REF], the zonal method is used, in which radiative models are used to calculate configuration factors Cfn,m among layers, which are then used to obtain the temperatures by

	inversion of the heat transfer equilibrium equation matrix [43][44][95][105]. The configuration
	factors are evaluated directly by counting rays [43][44][104], i.e.:
			𝐶𝑓 𝑛,𝑚 =	𝑁 𝑛,𝑚 𝑁 𝑛	;	𝐶𝑓 𝑛,𝐴 =	𝑁 𝑛,𝐴 𝑁 𝑛	;	𝐶𝑓 𝑛,𝐵 =	𝑁 𝑛,𝐵 𝑁 𝑛	;	(6.33a-c)
			𝐶𝑓 𝐴,𝑚 =	𝑁 𝐴,𝑚 𝑁 𝐴	;	𝐶𝑓 𝐴,𝐴 =	𝑁 𝐴,𝐴 𝑁 𝐴	;	𝐶𝑓 𝐴,𝐵 =	𝑁 𝐴,𝐵 𝑁 𝐴	;	(6.33d-f)
	With reference to (6.33a-c), Nn = 10 5 is the total number of rays cast from random points in the n th
	layer (for a grand total of 5•10 6 rays) and Nn,m, Nn,A and Nn,B are the number of rays cast from the n th
	layer that are absorbed respectively in the m th slice, at boundary A or at boundary B. It is worth noting
	that ∑	50 𝑚=1	𝑁 𝑛,𝑚

Table 6 .

 6 2 -General presentation of the 5 porous morphologies considered.

			Test summary
			α2	n2	L
			0.69	1.40	5
			0.69	1.70	5
			0.69	2.00	5
			1.37	1.40	3
			1.37	1.70	3
			1.37	2.00	3
			2.74	1.40	2
			2.74	1.70	2
			2.74	2.00	2
			ε = 0.89
			Normal diameter distribution
			CVd = 0.10
			t = 0.5
			k = 0.16
		Result synthesis
		AVERAGE ERROR
	HPA	INV HPA		HPA+
	Test	8.52%	5.18%		3.85%
	Test	17.95%	13.00%		6.93%
	Test	20.46%	12.26%		2.99%
		MAXIMUM ERROR
	HPA	INV HPA		HPA+
	Test	14.32%	11.84%		6.67%
	Test	35.83%	31.93%		10.49%
	Test	37.63%	30.07%		9.24%
	Table 6.3 -Data and results for Morphology #1
	(High porosity open cell foam -Ref. Paragraphs 3.2, 3.4, 3.6)

  And, given the typical values of the variables on the right side, we can expect

		𝑑𝐶 2 𝑑𝐶 1	=	(𝛼 2 𝐶 2 + 𝜎 21 ) (1 -𝐶 1 )𝛼 2	(7.23)
				𝑑𝐶 1 𝑑𝐶 2	to be smaller than
	its reciprocal	𝑑𝐶 2 𝑑𝐶 1

  Where for the sake of compactness I2+t is defined as emitted and reflected radiation associated with phase 2 and the trapped phase, already taking into account the weighted sum of their intensities, as follows:

	Becomes:
									𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚
	I 1 (𝑟 𝑤 + , Ω) = (𝜀 𝑠 𝑛 1 2 𝐵 0 + 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚 𝑡ℎ𝑒 𝑤𝑎𝑙𝑙 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 1 𝜌 𝑠 𝜋 ∫ I 1 (𝑟 𝑤 , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′ 𝑛 𝑤 •Ω ′ <0 ⏞	)	+	𝑡ℎ𝑒 𝑤𝑎𝑙𝑙 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2 𝑎𝑠𝑠𝑖𝑔𝑛𝑒𝑑 𝑡𝑜 𝑝ℎ𝑎𝑠𝑒 1 (1 -𝐶 1 )𝑓 2 𝐶 2 𝑓 1 𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω) ⏞	{𝑛 𝑤 • Ω > 0}	(7.27a)
						𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑎𝑛𝑑 𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑖𝑜𝑛 𝑓𝑟𝑜𝑚
						𝑡ℎ𝑒 𝑤𝑎𝑙𝑙 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2
						𝑎𝑠𝑠𝑖𝑔𝑛𝑒𝑑 𝑡𝑜 𝑝ℎ𝑎𝑠𝑒 2
		I 2 (𝑟 𝑤 + , Ω) = (1 -𝐶 2 )𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω) ⏞	{𝑛 𝑤 • Ω > 0}	(7.27b)
					I 𝑡 (𝑟 𝑤 + , Ω) = 𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω)		{𝑛 𝑤 • Ω > 0}	(7.27c)
	𝑇𝑜𝑡𝑎𝑙 𝑟𝑎𝑑𝑖𝑎𝑡𝑖𝑜𝑛
	𝑒𝑚𝑖𝑡𝑡𝑒𝑑 𝑎𝑛𝑑
	𝑟𝑒𝑓𝑙𝑒𝑐𝑡𝑒𝑑
	𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω) ⏞ 𝑖𝑛 𝑝ℎ𝑎𝑠𝑒 2	= 𝜀 𝑠 𝑛 2 2 𝐵 0 +	𝜌 𝑠 𝜋	𝑛 𝑤 •Ω ′ <0 ∫ [(1 -𝐶 1 )𝐼 2 (𝑟 𝑤 + , Ω′) ⏞ 𝑁𝑜𝑛-𝑡𝑟𝑎𝑝𝑝𝑒𝑑 𝑝ℎ𝑎𝑠𝑒 2	+ 𝐶 1 𝐼 𝑡 (𝑟 𝑤 + , Ω ′ ) ⏞ 𝑇𝑟𝑎𝑝𝑝𝑒𝑑 𝑝ℎ𝑎𝑠𝑒	] |𝑛 𝑤 • Ω ′ |𝑑Ω ′	{𝑛 𝑤 • Ω > 0}	(7.28)
	𝐼 𝑖 (𝑟 𝑤 + , Ω) = 𝜀 𝑠 𝑛 𝑖 2 𝐵 0 +	𝜌 𝑠 𝜋	𝑛 𝑤 •Ω ′ <0 ∫ 𝐼 𝑖 (𝑟 𝑤 + , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′	{𝑛 𝑤 • Ω > 0, 𝑖 = 1, 2}	(7.26)

Table 7 .

 7 1. Different treatments applied to different fractions of radiation. = I 1 𝑓 1 + I 2 𝑓 2 (1 -𝐶 1 ) + I 𝑡 𝑓 2 𝐶 1 𝑓 2 (1 -𝐶 1 )𝑛 2 2 (𝛼 2 𝐶 2 + 𝜎 21 ) = 𝜎 12 𝑓 1 𝑛 1

	Ω • ∇I 1 = 𝑛 1 2 𝛼 1 𝐵 1 0 +	𝑓 2 (1 -𝐶 1 ) 𝑓 1	𝐶 2 𝑛 2 2 𝛼 2 𝐵 2 0 -𝛽 1 𝐼 1 +	𝜎 11 4𝜋	4𝜋 ∫ 𝐼 1 (Ω ′ )Φ 11 (Ω ′ , Ω)𝑑Ω ′	(7.14a)
					+	𝑓 2 (1 -𝐶 1 ) 𝑓 1	𝜎 21 4𝜋	∫ 𝐼 2 (Ω′)Φ 21 (Ω ′ , Ω)𝑑Ω′ 4𝜋
	Ω • ∇I 2 = 𝑛 2 2 𝛼 2 𝐵 2 0 (1 -𝐶 2 ) -𝛽 2 𝐼 2 +	𝜎 22 4𝜋	∫ 𝐼 2 (Ω′)Φ 22 (Ω ′ , Ω)𝑑Ω′ 4𝜋
								(7.14b)
							+	𝑓 1 𝑓 2 (1 -𝐶 1 )	𝜎 12 4𝜋	∫ 𝐼 1 (Ω′)Φ 12 (Ω ′ , Ω)𝑑Ω′ 4𝜋
	Ω • ∇I 𝑡 = 𝛼 2 𝐵 2 0 -(𝛼 2 + 𝜎 𝑡 )I 𝑡 +	4𝜋 𝜎 𝑡	4𝜋 ∫ I 𝑡 (Ω′)Φ 𝑡 (Ω ′ , Ω)𝑑Ω′	(7.14c)
								(7.14d)
	-Prescribed directional intensity at boundary:
					I 1 (𝑟 𝑤 + , Ω 𝑤 ) = 𝐼 𝑤 [1 + 𝐶2	𝑓 2 (1 -𝐶1) 𝑓1	]	(7.25a)
							I 2 (𝑟 𝑤 + , Ω 𝑤 ) = 𝐼 𝑤 (1 -𝐶2)	(7.25b)
								I 𝑡 (𝑟 𝑤 + , Ω) = 𝐼 𝑤	(7.25c)
	-Diffusely emitting/reflecting wall boundary:
	I 1 (𝑟 𝑤 + , Ω) = (𝜀 𝑠 𝑛 1 2 𝐵 𝑣 +	𝜌 𝑠 𝜋	∫ I 1 (𝑟 𝑤 , Ω ′ )|𝑛 𝑤 • Ω ′ |𝑑Ω ′ 𝑛 𝑤 •Ω ′ <0	) +	(1 -𝐶 1 )𝑓 2 𝐶 2 𝑓 1	𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω) {𝑛 𝑤 • Ω > 0}	(7.27a)
			I 2 (𝑟 𝑤 + , Ω) = (1 -𝐶2)𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω)	{𝑛 𝑤 • Ω > 0}	(7.27b)
						I 𝑡 (𝑟 𝑤 + , Ω) = 𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω)	{𝑛 𝑤 • Ω > 0}	(7.27c)
	𝐼 2+t 𝑤 (𝑟 𝑤 + , Ω) = 𝜀 𝑠 𝑛 2 2 𝐵 𝑣 +	𝜌 𝑠 𝜋	𝑛 𝑤 •Ω ′ <0 ∫ [(1 -𝐶 1 )𝐼 2 (𝑟 𝑤 + , Ω′) + 𝐶 1 𝐼 𝑡 (𝑟 𝑤 + , Ω ′ )]|𝑛 𝑤 • Ω ′ |𝑑Ω ′	{𝑛 𝑤 • Ω > 0}	(7.28)
	-Compatibility condition:

I 2 (7.21) Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2018LYSEI005/these.pdf © [S. Cunsolo], [2018], INSA Lyon, tous droits réservés

Table 7 .

 7 2. General presentation of the 6 porous morphologies considered.

					Test summary
				α2	n2	L
					0.69	1.40	5
					0.69	1.70	5
					0.69	2.00	5
					1.37	1.40	3
					1.37	1.70	3
					1.37	2.00	3
					2.74	1.40	2
					2.74	1.70	2
					2.74	2.00	2
				ε = 0.89
				Normal diameter distribution
				CVd = 0.10
				t = 0.5
				k = 0.16
			Result synthesis
		AVERAGE ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	8.63%	3.85%	11.18%	5.57%		2.23%
	Test	17.86%	6.93%	13.35%	5.62%		2.06%
	Test	20.55%	2.99%	10.54%	7.30%		1.50%
		MAXIMUM ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	14.12%	6.67%	14.72%	7.31%		5.43%
	Test	35.37%	10.49%	15.41%	8.13%		3.79%
	Test	34.10%	9.24%	12.94%	10.92%		2.96%

Table 7 .

 7 3 -Data and results for Morphology #1 (High porosity open cell foam -Ref. Paragraphs 3.2, 3.4, 3.6)

	Test summary
	α2	n2	L
	2.01	1.40	5
	2.01	1.70	5
	2.01	2.00	5
	4.02	1.40	3
	4.02	1.70	3
	4.02	2.00	3
	8.05	1.40	2
	8.05	1.70	2
	8.05	2.00	2
	ε = 0.96
	Normal diameter distribution
	CVd = 0.10
	t = 0.4
	k = 0.6

Table 7 .

 7 4 -Data and results for Morphology #2 (High porosity open cell foam -Ref. Paragraphs 3.2, 3.4, 3.6)

					Test summary
				α2		n2	L
				0.48	1.40	5
				0.48	1.70	5
				0.48	2.00	5
				0.95	1.40	3
				0.95	1.70	3
				0.95	2.00	3
				1.90	1.40	2
				1.90	1.70	2
				1.90	2.00	2
				ε = 0.85
				Lognormal diameter distribution
				GCVd = 0.3
			Result synthesis
		AVERAGE ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	10.83%	3.67%	4.03%	1.89%		2.95%
	Test	21.91%	5.92%	5.87%	2.42%		2.52%
	Test	25.81%	2.59%	3.47%	2.13%		2.38%
		MAXIMUM ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	18.85%	7.22%	5.92%	2.42%		3.60%
	Test	45.35%	11.87%	7.05%	3.79%		4.35%
	Test	50.75%	6.95%	5.13%	2.56%		3.16%

Table 7 .

 7 5 -Data and results for Morphology #3 (Open cell foam -Ref. Paragraphs 3.2, 3.3)

					Test summary
				α2		n2	L
				0.48	1.40	5
				0.48	1.70	5
				0.48	2.00	5
				0.95	1.40	3
				0.95	1.70	3
				0.95	2.00	3
				1.90	1.40	2
				1.90	1.70	2
				1.90	2.00	2
				ε = 0.70
				Lognormal diameter distribution
				GCVd = 0.3
			Result synthesis
		AVERAGE ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	20.37%	4.74%	6.32%	0.88%		3.45%
	Test	39.56%	15.07%	14.33%	0.89%		3.38%
	Test	40.23%	7.63%	11.13%	0.80%		2.27%
		MAXIMUM ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	32.90%	12.58%	10.81%	1.53%		4.38%
	Test	77.55%	22.96%	17.34%	1.37%		4.63%
	Test	73.85%	18.70%	13.84%	1.32%		3.75%

Table 7 .

 7 6 -Data and results for Morphology #4 (Open cell foam -Ref. Paragraphs 3.2, 3.3)

					Test summary
				α2		n2	L
				0.34	1.40	5
				0.34	1.70	5
				0.34	2.00	5
				0.68	1.40	3
				0.68	1.70	3
				0.68	2.00	3
				1.36	1.40	2
				1.36	1.70	2
				1.36	2.00	2
				ε = 0.85
				Lognormal diameter distribution
				GCVd = 0.3
			Result synthesis
		AVERAGE ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	10.22%	9.80%	48.88%	18.63%		3.31%
	Test	23.53%	10.12%	73.97%	22.74%		2.69%
	Test	24.84%	3.92%	61.90%	24.11%		2.90%
		MAXIMUM ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	18.15%	11.56%	71.17%	25.91%		6.85%
	Test	41.29%	16.25%	78.40%	36.84%		4.92%
	Test	41.28%	8.50%	68.88%	34.37%		8.09%

Table 7 .

 7 7 -Data and results for Morphology #5 (Closed cell foam -Ref. Paragraphs 3.2, 3.3)

					Test summary
				α2		n2	L
				4.00	1.40	5
				4.00	1.70	5
				4.00	2.00	5
				8.00	1.40	3
				8.00	1.70	3
				8.00	2.00	3
				16.00	1.40	2
				16.00	1.70	2
				16.00	2.00	2
				ε = 0.98
				Lognormal diameter distribution
				GCVd = 0.3
			Result synthesis
		AVERAGE ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	7.23%	3.50%	88.44%	25.88%		2.97%
	Test	9.05%	2.10%	103.65%	15.86%		1.80%
	Test	11.06%	2.61%	67.38%	25.15%		1.05%
		MAXIMUM ERROR
	HPA	HPA+	MPA	MPA GRTE	MPA+
	Test	11.67%	5.02%	101.14%	35.41%		6.64%
	Test	15.56%	4.35%	112.76%	24.35%		2.69%
	Test	23.34%	6.76%	74.62%	35.46%		2.58%

Table 7 .

 7 7 -Data and results for Morphology #6 (Closed cell foam -Ref. Paragraphs 3.2, 3.4, 3.6)

  ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON Les transferts thermiques par rayonnement dans des mousses sont modélisés à partir de la morphologie du matériau et des propriétés de sa phase solide. Dans ce travail de thèse, une attention particulière est portée sur les modèles radiatifs de Monte Carlo. Les différentes approches d'homogénéisation telles que « Homogeneous Phase » (HPA) and « Multi Phase » (MPA) sont discutées et comparées. Des développements novateurs sont proposés pour améliorer la précision des résultats.Nos avancées permettent de générer numériquement trois types de mousses périodiques couvrant un large domaine de matériaux cellulaires: mousse à pores fermés à haute porosité, mousse à cellules ouvertes à basse et haute porosité. Pour ces dernières, des comparaisons morphologiques avec des données expérimentales tomographiques, montrent des résultats satisfaisants et tendent à valider nos modèles de génération. Des mousses dont la phase solide est opaque ont tout d'abord été étudiées. Nos simulations ont permis de trouver de nouvelles lois analytiques précises permettant d'estimer les propriétés radiatives de ces matériaux à partir de données morphologiques. Ensuite, nous avons considéré des mousses, dont la phase solide est semi transparente. La modélisation du transfert radiatif au sein de ces milieux cellulaires est plus complexe. Les méthodes de modélisation des propriétés radiatives de la littérature, HPA et MPA, sont testées. Des simulations de Monte Carlo directes dans les matériaux ont permis de mettre en exergue les limites de ces modèles. Des modèles novateurs ont été proposés afin d'améliorer la précision des simulations. Ils sont basés sur une méthode hybride directe-inverse et une modification de l'équation de transfert radiatif classique. Ces nouveaux modèles (HPA+ et MPA+) ont été testés sur un ensemble varié de morphologies générées numériquement. Les modèles améliorés sont systématiquement plus précis que les modèles existants MOTS-CLÉS : Rayonnement, Tomographie, Monte Carlo, Morphologie, Numérique, Multi-Phase Approach

	NOM : CUNSOLO	DATE de SOUTENANCE : 23/01/2018
	(avec précision du nom de jeune fille, le cas échéant)
	Prénoms : Salvatore
	TITRE : RADIATIVE PROPERTIES COMPUTATIONAL MODELING OF POROUS CELLULAR MATERIALS
	NATURE : Doctorat	Numéro d'ordre : 2018LYSEI005
	Ecole doctorale : ED 162 -MEGA
	Spécialité : Thermique et Energétique
	RESUME : Laboratoire (s) de recherche: LaMCoS
	Directeur de thèse: Dominique Baillis
	Président de jury :
	Composition du jury :
	Dombrovsky, Leonid	Chief Researcher JIHT	Rapporteur
	Minea, Alina Adriana	Prof. TUIASI	Rapporteuse
	Enguehard, Franck	Prof. CentraleSupélec	Examinateur
	Rosato, Antonio	Prof. UNICAMPANIA	Examinateur
	Baillis, Dominique	Prof. INSA-Lyon	Directrice de thèse
	Bianco, Nicola	Prof. UNINA	Co-directeur de thèse

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2018LYSEI005/these.pdf © [S. Cunsolo], [2018], INSA Lyon, tous droits réservés

Acknowledgments

This thesis has been financed by a grant from Region Rhone Alpes and also by a grant from Université Franco Italienne / Università Italo Francese (UFI/UIF).

I would also like to extend my gratitude to

With reference to (6.33d-f), NA = 10 6 is the total number of rays cast from boundary A and NA,m, NA,A and NA,B are the number of rays cast from boundary A that are absorbed respectively in the m th layer, Solving this system, temperatures T1…50 can be determined. The highly automatized optimized linear system solver in MATLAB [106] is used to solve the system, with no further adaptations necessary due to the small size of the input matrices.

Finally, the total heat flux Q can be calculated as:

The error is quantified as: ERROR = |Q/q -1| (6.38) Where Q is the flux calculated in Full Direct Monte Carlo and q is the flux calculated with homogenized methods.