Christelle Lusso

Xavier Chateau

Disordered monodisperse 3D open-cell foams: Elasticity, Thermal conductivity and Permeability

Keywords: homogenization, finite element, porous media, disordered foam, elasticity, thermal conductivity, permeability

This work is devoted to the modeling of the overall properties of 3D disordered monodisperse open cell foams manufactured by solidification of liquid foams. First, we develop a new and flexible meshing procedure to generate Representative Volume Elements (RVE) of disordered wet foam. We consider solid foam RVEs arising from solidification of the liquid Plateau border with wide range of volume fraction varying from 2.5% to 20%. Effective properties of foams such as elasticity, thermal conductivity and permeability are determined by averaging the solution of boundary value problems set on the RVE. The local problems are discretized and solved by the finite element method. We evaluate the accuracy of the resultant estimates, with respect to the size of the RVEs, by conducting a periodic study on a Kelvin's foam structure. In order to validate our foam numerical model, we compare results for the overall elasticity to the estimates of Gibson and Ashby, and results for the overall permeability to classical Kozeny-Carman equations, leading to good agreement. Finally, we propose new closed form estimates for both the elastic modulus and the thermal conductivity.

Introduction

Foam materials are ubiquitous, ranging from shaving cream to concrete. The main features of foam materials originate in their component's behavior, from liquid to solid, as well as their structure. Foams materials are involved in a broad range of applications such as firefighting, lightweight structures, cushions seats, foods, cosmetics, ... [START_REF] Bikerman | Foams[END_REF][START_REF] Prud | Foams: Theory, Measurements and Applications[END_REF]. Industrial interests together with academic researches motivate the development of foam sciences which draw chemists, physicists and mathematicians. Indeed, understanding the mechanical properties of foam materials is required in order to design new equipment and to improve performance of existing ones [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. In this context, the foam sciences span from fluid mechanics [START_REF] Koehler | A Generalized View of Foam Drainage: Experiment and Theory[END_REF][START_REF] Verbist | The foam drainage equation[END_REF] to solid mechanics and [START_REF] Areiasa | Analysis of experimentally assessed EVA foams with mixed solid-shell elements capable of very large strains[END_REF] and soft matter physics [START_REF] Kraynik | Foam flows[END_REF][START_REF] Kraynik | Foam micromechanics: Structure and rheology of foams, emulsions, and cellular solids[END_REF][START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF][START_REF] Weaire | Stress and strain in liquid and solid foams[END_REF].

The structure of foam materials involves several lengths ranging from microscopic to macroscopic scales. The microscopic scale exhibits heterogeneities corresponding to interfaces between the bubbles with thickness of a few nanometers, Plateau borders with lateral dimensions of a few micrometers and gas bubbles with diameters ranging from 10 µm to few millimeters. At the macroscopic scale (characterized by a typical length scale of the order of 10 bubbles size), the heterogeneous features tend to average out, conferring an homogeneous aspect to foamy materials. Thus foam can be modeled as an homogeneous material at the macroscopic scale. The macroscopic behavior describing the overall response to a macroscopic loading is derived from the behavior and spatial arrangement of material's components at the microscopic scale such that macroscopic state variables linked by the macroscopic state law must be equal to average of their microscopic counterpart computed over a Representative Volume Element (RVE) of the heterogeneous material. Predicting the macroscopic behavior of heterogeneous materials from that of their constituents properties and microstructure is the aim of homogenization methods [START_REF] Allaire | homogenization of the Stokes flow in a connected porous medium[END_REF][START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Dormieux | Microporomechanics[END_REF][START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF][START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF][START_REF] Francfort | Homogenization and optimal bounds in linear elasticity[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF][START_REF] Tartar | Incompressible fluid flow in a porous medium -convergence of the homogenization process[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF][START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF] consisting in three main steps: the representation or description of the material properties, the definition of a localization problem and the homogenization procedure [START_REF] Zaoui | Continuum micromechanics: survey[END_REF].

In this work, we investigate the macroscopic properties of open cell monodisperse disordered solid foams in the framework of homogenization techniques. We restrict ourselves to solid foams arising from the solidification of liquid contained in the Plateau borders of a liquid foam. Therefore, samples of disordered foams are generated with Surface Evolver [START_REF] Brakke | The Surface Evolver[END_REF][START_REF] Phelan | Computation of equilibrium foam structures using the Surface Evolver[END_REF], a program for modeling the shape of liquid interfaces. Once a stable equilibrium has been obtained, we assume that the liquid phase turns to a solid state without any geometrical changes. Therefore the liquid foam turns to a solid foam having open cells structure (we also assume that the films separating bubble one another are destroyed during the solidification process). Such a foam can be generated from porcine gelatin which is a thermoreversible gel: above 29 o C, it is a liquid and below this temperature it is a soft solid [START_REF] Khidas | Critical size effect of particles reinforcing foamed composite materials[END_REF]. Then, the homogenization problems are set on the solid foam samples assuming that applied loadings do not induce any change of the foam structure. From a practical point of view, this condition is fullfilled for the thermal conduction problem only if the prescribed temperature is lower than the solid phase melting temperature.

Similar homogenization problems have been solved for the elastic properties of random open-cells foam in [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF] for solid foams. Enhancing what have been done in [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF], we manage to handle foam samples complying with the Plateau border rules [START_REF] Weaire | The Physics of Foams[END_REF].

Before we present in detail the work we have performed, it is worth summarizing the three steps homogenization procedure we have employed.

The first step is to describe the studied material. To do so, we model the microstructure of a RVE of an open cell disordered monodisperse foam. The formal definition of the RVE entails a condition of separation of scales (micromacro) [START_REF] Zaoui | Continuum micromechanics: survey[END_REF] assumed all along this work. As solid foams considered in this work are obtained from solidification of liquid foams, the starting point is to generate liquid foam samples. The microstructure of a liquid foam depends on the liquid volume fraction, yielding either a dry foam or a wet foam. In the dry limit, the liquid fraction is zero and the bubbles are compacted together forming polyhedra. Below the wet limit, the bubbles remain approximately spherical with some surface contacts. Thus, the wet foam structure is modeled as a continuous network of convex polyhedral surfaces. The meshing procedure adopted is based on a Voronoï tessellation, which in turn is based on random packings algorithm [START_REF] Krauth | Introduction to monte carlo algorithms[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. Once generated, the Voronoï tessellation is evolved towards a physical structure of dry foam [START_REF] Kraynik | Structure of random monodisperse foam[END_REF][START_REF] Kraynik | Foam microrheology: from honeycombs to random foams[END_REF] by a process called relaxation. The numerical procedure for relaxing the Voronoï tessellation is handle with the Surface Evolver software [START_REF] Brakke | The Surface Evolver[END_REF], which allows to evolve surfaces shaped by forces and constraints. The resultant stable structure models a disordered dry foam material. In this study, we focus on the monodisperse case so each cell is constrained to have an equal volume in the relaxation algorithm. Next, the conversion of dry foams to wet foams with arbitrary liquid fraction is carried out by first replacing the lines modeling the Plateau borders by liquid-carrying channels, and then computing equilibrium shape of the wet foam with Surface Evolver. The generated wet foam is transformed into an open cell solid foam by considering the Plateau borders as solid and removing the films located in between the bubbles. The solid foam is thus made of only the solid Plateau borders. Doing so, we numerically simulate elaboration of samples of open cell solid foams obtained through solidification of the liquid Plateau borders without change of shape. We also generate Kelvin periodic open cell monodisperse foam with the Surface Evolver to evaluate the accuracy of homogenized properties computed on finite-size RVE.

The localization step consists in solving a boundary value problem set on the RVE. The problem one has to solve depends on the property of interest (elasticity, thermal conductivity or permeability). For elasticity, the Plateau borders material behavior is modeled at the microscopic scale by a linearly elastic isotropic homogeneous law (Hooke's law) and the porous space is empty. The RVE is submitted to a no body force equilibrium condition and a prescribed displacement on the outer solid boundary accounting for the macroscopic strain. To determine the overall thermal conductivity of the solid foam, the heat transfer over the Plateau borders is modeled at the microscopic scale by the Fourier's law and it is assumed that the gas in the bubbles is a perfect insulating material at rest (no convection). A prescribed macroscopic temperature gradient is enforced on the outer solid surface and there is no heat source in the domain filled by the solid foam. Finally, we study the flow of a Newtonian fluid through the porous space of the solid foam to estimate its permeability. The motion of the fluid filling the porous space is described by the Stokes equation with no body force with a no-slip condition at the solid-liquid interface and a prescribed macroscopic pressure gradient on the outer porous surface. As it is well known, homogenization of the Stokes equations leads to the Darcy's law which describes at the macroscopic scale inertialess Newtonian fluid flows through a rigid porous medium [START_REF] Allaire | homogenization of the Stokes flow in a connected porous medium[END_REF][START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF][START_REF] Tartar | Incompressible fluid flow in a porous medium -convergence of the homogenization process[END_REF][START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF].

All these boundary value problems are discretized and solved by the Finite Element Method, with a MPI formulation for the resolution of the Stokes problem set on the porous space [START_REF] Bercovier | A finite element method for incompressible non-newtonian flows[END_REF][START_REF] Hecht | New development in FreeFem++[END_REF].

The last step in the homogenization approach is to average the solutions of the boundary value problems to compute the macroscopic state variables associated to microscopic fields. Macroscopic properties of the homogeneous material, equivalent to the studied heterogeneous material, are determined from these averaged quantities. More details concerning theses three steps' homegenization approach can be found in [START_REF] Zaoui | Continuum micromechanics: survey[END_REF].

The outline of the paper is as follows. The Section 2 is devoted to the representation step in the homogenization approach. The numerical procedure carried out to generate the microstructure meshing of foam material is detailed. In Section 3, we outline the localization step of the homogenization approach. The formulation of the local problems on the RVE are precised. In Section 4, we describe the homogenization step. In Section 5, the discretization of the problems is precised and a numerical study is conducted. More precisely, we perform a numerical study aiming to estimate the accuracy of the results as a function of the RVE size. This estimation is achieved through an implementation of the homogenization approach in a periodic configuration. Finally, Section 6 provides concluding remarks and observations.

Numerical modeling of the foam structure

The prerequisite to predicting overall properties of a material is to model its microstructure (geometry and behavior). In the framework of homogenization methods, overall behavior of heterogeneous material are predicted from the solution of a boundary value problem defined on a RVE of the material. This section outlines the procedure adopted to generate the discretized RVE, yielding the mesh of the microstructure. We present the numerical modeling of both disordered and periodic structures of foams.

Introduction

A foam material consists of gas cells (bubbles) enclosed in a liquid or solid. One of the salient feature of a foam is the volume fraction of the continuous phase (liquid or solid), i.e. the ratio of the continuous phase volume to the volume of the whole foam material, denoted by φ in the sequel. In this section we consider only liquid foam with gas bubbles separated by thin films made of two parallel interfaces between which a thin layer of liquid is enclosed. These films join together to form Plateau borders containing most of the liquid phase. The structure of liquid foams mainly depends on the value of the liquid volume fraction φ. On one hand, when the liquid volume fraction is low, we deal with dry foams.

Dry foams are made of polyhedral bubbles separated by thin films. Indeed, when the liquid volume fraction tends towards zero both the Plateau borders and the films thickness decrease to zero. On the other hand, when increasing liquid volume fraction, the Plateau borders thickness increases and thin films aera decreases, leading to more rounded bubbles. For volume fraction larger than the so-called wet limit, film separating bubbles does no more exist and bubbles shape is spherical: the gas-liquid mixture is no more a foam but rather a suspension of bubbles in a continuous liquid phase. For disordered monodisperse foam, the wet limit is reached for φ = 0.36, which corresponds to the random close packing density of the gas bubbles. In this work, we consider values of φ in [0.025, 0.2] which is between the dry and the wet limit. Foams can thus be modeled as a continuous network of convex polyhedral surfaces.

As homogenization problems we are interested in are set on solid foam, it is assumed in the sequel that liquid phase solidificates without modifying the shape of the liquid-gas interface. Furthemore, we consider that films separating bubbles are fully destroyed during the solidification of the Plateau borders. Thus, the solid foam RVE domain is perforated by pores formed by the network of bubbles of gas. The set of pores may be designated as the porous space. The porous space can be filled by a fluid that cannot penetrate the solid Plateau borders. In this work, we state problems on both the Plateau border and the porous space. The entire domain is denoted by Ω, including the solid Plateau borders network Ω s and the porous space denoted by Ω p . The Plateau borders and the porous space are complementary, that is Ω = Ω s ∪ Ω p and Ω s ∩ Ω p = ∅.

Disordered foam

First we model the microstructure of a disordered liquid foam. As an initial input, we place ourselves in the dry limit, then the structure is expanded to a wet one by thickening the Plateau borders. The structure of a dry foam is modeled as a set of convex polyhedra packed to fill space. The initial configuration of the foam RVE (i.e. before decreasing its energy with the Surface Evolver) occupies a spherical domain of diameter ℓ. The Voronoï tessellation algorithm is adopted to generate a partition of a sphere into convex polyhedral cells. The Voronoï tessellation is based on a spatial distribution of points S = x 1 , ..., x N x i ∈ R 3 with N the number of cells. The spatial distribution of points is based on random packing algorithms [START_REF] Krauth | Introduction to monte carlo algorithms[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. Herein, the adopted algorithm ensures disordered random Voronoï tessellation with approximately equal volume. At this stage, the meshing is performed with the open source program Neper [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF].

The next step is to evolve the surface of the Voronoï tessellation towards physically relevant geometry modeling the foam microstructure [START_REF] Kraynik | Structure of random monodisperse foam[END_REF][START_REF] Kraynik | Foam microrheology: from honeycombs to random foams[END_REF]. This step is called relaxation and the algorithm provides a stable structure of foam. It is reminded that structure of a foam at equilibrium is governed by three rules:

1. geometrical stability criterion, 2. topological equilibrium conditions, 3. minimization of the energy. The numerical procedure for relaxing the Voronoï tessellation is based on the Surface Evolver software developed by Ken [START_REF] Brakke | The Surface Evolver[END_REF] [START_REF] Brakke | The Surface Evolver[END_REF][START_REF] Phelan | Computation of equilibrium foam structures using the Surface Evolver[END_REF]. The Surface Evolver is specifically designed to find minimal energy configurations of surfaces under constraints. The meshed surface is associated with a constant surface energy per unit area, the surface tension, and is evolved using a gradient descent method of minimization. The vertices are shifted by reducing the energy, thus leading to approximately spherical cells in the absence of geometrical crowding. The algorithm adopted for the foam relaxation is inspired by Andrew Kraynik's algorithm [START_REF] Kraynik | Structure of random monodisperse foam[END_REF][START_REF] Kraynik | Foam microrheology: from honeycombs to random foams[END_REF]. We implement an enhanced version of the Kraynik's algorithm as follows. In the first step, the surface is evolved towards some local minimum in surface energy by iterating Conjugate Gradient descent method of optimization. Doing so, the polyhedra become surfaces with constant Gaussian curvature. Besides, the cells are constrained to have equal volume, thus modeling a monodisperse disordered foam. Then the surface is modified to comply with the geometrical stability condition. To this end, a criterion has been established by Andrew Kraynik based on the length of cell edges (see [START_REF] Kraynik | Structure of random monodisperse foam[END_REF]). In our algorithm, we ensure stability by eliminating short edges of length λ and tiny triangles of area β satisfying

λ < ǫ λ , β < ǫ β ,
where λ is the average length, β is the average area, and ǫ ∈ R is the threshold that restricts the edges length and the triangles surface with typically 0.01 ≤ ǫ ≤ 0.1 for stability [START_REF] Kraynik | Structure of random monodisperse foam[END_REF]. On the next step, the microstructure must satisfy topological equilibrium conditions, namely the Plateau's laws [START_REF] Weaire | The Physics of Foams[END_REF]. We remind that the topological requirements on edge and face connectivity in Plateau's laws are: (a) each film has constant mean curvature, in order to balance the pressure difference between adjacent cells separated by a film (Laplace's law), (b) 3 films meet at each cell edge and (c) 4 edges meet at each cell vertex. The topological equilibrium requirements are fullfilled by removing the nodes not complying with Plateau's law with a Surface Evolver's command. Furthermore, the mesh is regularized, the facets of the triangulation are shaped as equilateral as possible and each vertex is moved to the average position of its neighboring vertices. Indeed, during the relaxation procedure, some edges shrink to zero length and faces shrink to zero area, producing irregular skinny triangles. They are eliminated thanks to the criterion described above. We iterate this algorithm until the number of defaults (as regards to the Plateau's law) tends to 0. This algorithm allows the surface of the Voronoï tessellation to evolve towards a stable structure of foam by taking about O(10 4 ) iterations of minimization of the energy. The figure 1 depicts a Voronoï tessellation of a spherical domain with N = 64 cells before (left) and after (right) relaxation towards a physically relevant structure of foam. A wet foam structure can then be inferred starting from the relaxed dry foam structure. The extension to a wet foam structure with arbitrary volume fraction is achieved through the adjonction and thickening of the Plateau borders. Pipes that delineate the edges of the polyhedron are created whereas the films remain devoid of thickness. The surface tension of the films is allocated with twice the surface tension of the surfaces shared by cells and Plateau borders. An edge spread parameter determines the initial width of the additional Plateau borders, typically about 20% of the average edge length. Subsequently the Plateau border's volume is gradually increased and adjusted. The final stage of minimizing the energy is very sensitive to the parameters and requires to be closely controlled. The parameters are tuned by hand in each case (N, φ) with N the number of cells and φ the liquid volume fraction. The procedure involves minimizing the energy, regularizing the mesh and weeding out tiny elements, together with refining the mesh. The iteration stopping tolerance in the relative variation of the energy is 10 -4 . The figure 2 depicts the wet foam structure of a relaxed tessellation with the adjonction of the Plateau border network. On the left figure, the surface is generated from the relaxed Voronoï tessellation of N = 64 cells (right picture in figure 1) by adding a Plateau border (red pipe). On the right figure, the surface is optimized with no refinement. Surface Evolver scripts have been written to extract the surface mesh of the Plateau border network as well as the surface mesh of the porous space. Finally, the disjoint volume domains Ω s and Ω p are meshed: an unstructured 3d tetrahedral mesh is generated within the volume defined by the connectivity of the 2d triangular surface mesh. The meshing procedure is achieved using Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. For the sake of describing the boundaries of the domains, we introduce the following notations (the same notations hold for the periodic material)

Ω = Ω s ∪ Ω p , ∂Ω = Γ ext s ∪ Γ ext p , (2.1)
∂Ω s ∩ ∂Ω p = Γ int . (2.2)
where ∂Ω (resp. ∂Ω s , ∂Ω p ) denotes the boundary of the domain Ω (resp. Ω s , Ω p ), Γ ext s = ∂Ω s ∩ ∂Ω the outer boundary of the solid domain, Γ ext p = ∂Ω p ∩∂Ω the outer boundary of the porous domain and Γ int the inner solid-porous interface.

The figure 3 depicts the mesh of a disordered monodisperse foam with N = 64 cells and a solid volume fraction φ = 5%. We represent both the Plateau border network (left) and the porous space (right). We mention that the meshing procedure is substantially time consuming. For instance, a N = 512 cells tessellation needs around 24 hours for the relaxation algorithm leading to the stable structure of foam, followed by about few hours for the optimization algorithm. As a consequence, this configuration is the largest implemented in this work. The figure 4 depicts a monodisperse disordered foam structure with N = 512 cells and solid volume fraction φ = 5%.

Periodic foam

In the periodic case, the microstructure is a periodic network generated by a unit cell filling the space by translation. The unit cell of a periodic network is not uniquely defined and can be arbitrarily selected. The determination of the unit cell is motivated by geometrical symmetries which can be taken in advantage to simplify the boundary conditions of the localization problem. In the periodic context, the overall homogeneous properties are determined from the solution of a boundary value problem set on the unit cell rather than on the RVE. The effective properties, deduced from the unit cell, are defined independently of the choice of a unit cell among all those which make possible to generate the same medium by periodic repetition [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Sanchez-Palencia | Homogenization Techniques for Composites Media[END_REF]. The periodic media allow to determine the required size of the RVE in order to provide an accurate estimate for the macroscopic behavior, by comparing the solution obtained on a periodic structure made of numerous cells with the exact solution obtained on the unit cell (see the subsection 5.2.1 below for more details). The results obtained on a periodic configuration thus allow to assess the accuracy of the results obtained on a disordered configuration of a porous media.

We consider the Kelvin structure [START_REF] Kelvin | The Kelvin problem : foam structures of minimal surface area[END_REF][START_REF] Matzke | The Three-Dimensional Shape of Bubbles in Foam-An Analysis of the Role of Surface Forces in Three-Dimensional Cell Shape Determination[END_REF] which is a periodic network of Kelvin cells. The Kelvin cell has symmetry along each direction of space. Taking advantage of the three symmetries allows to reduce the problem on one eighth of the Kelvin cell. The mesh of one eighth of the Kelvin cell is build with the Surface Evolver. Both the Plateau border and the porous space are generated separately. The figure 5 depicts the mesh of the eight of a Kelvin cell with volume fraction φ = 5%. We duplicate the eigth of the Kelvin cell in order to generate the entire Kelvin cell, then we duplicate the entire Kelvin cell with N p cells per axe in order to generate a N = N p 3 cells mesh. In this framework, the size of the RVE relies on N p , which is chosen large enough to satisfy the condition d ≪ ℓ where d is the size of the heterogenities and ℓ is the size of the RVE (see paragraph 5.2.1 and figure 7). Moreover, we conventionally assume that the size of the basis cell d is of the same order of magnitude as the size of the heterogeneities (bubbles). The figure 6 depicts the mesh of N p = 2 Kelvin cells periodic foam sample with volume fraction φ = 5%. As in the case of disordered foams 2.2, the final optimization procedure involves iteratively minimizing the energy, regularizing the mesh and weeding out tiny elements, together with refining the mesh. Finally, the disjoint domains Ω s and Ω p are meshed with Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF].

Statement of the localization problems

The starting point of the homogenization approach is the description of material properties at the microscopic scale. We prescribe these local properties of the RVE by setting boundary value problems on the RVE. The overall responses of the RVE, in response to external sollicitations, are dictated by the loading conditions. In the context of homogenization approaches, only the average of the loading influences the overall response [START_REF] Zaoui | Continuum micromechanics: survey[END_REF]. Thereafter, as usual, we consider uniform loading in deformation (elasticity), temperature gradient (thermal conductivity) and pressure gradient (permeability).

Hereafter, we consider the behavior of solid foam structures. The Plateau borders network Ω s is elastic or rigid and the porous space Ω p is filled with a fluid flowing around the solid Plateau borders, or at rest, depending on the property under consideration. We consider two configurations: (a) disordered material, (b) periodic material, requiring two different boundary conditions. We first write down the constitutive equations governing the material behavior, and the conservation equations. Then, the particular boundary conditions suited to each case (disordered and periodic) are described. Finally, the associated variational formulations are written.

The elasticity problem

The elasticity model is adopted to describe the deformation of the solid phase. In what follows, we present the equations governing the equilibrium state of a compressible isotropic elastic medium, in the context of small deformation. The solid domain Ω s is filled by a solid and deformable material. Over the complementary domain Ω p , the stress is zero.

The solid phase obeys the linear Hooke's law, its constitutive equation reads

σ = 2µε(u) + λ trace (ε(u))δ in Ω s , (3.1)
where σ denotes the Cauchy stress tensor and ε(u) = 1 2 (∇u+∇ T u) the linearized strain tensor. Moreover u denotes the displacement field, x = (x, y, z) is the position vector at the microscopic scale, µ and λ stand for the Lamé coefficients and δ denotes the second order identity tensor.

As we consider the equilibrium state of the RVE, the momentum balance law with no body force reads div σ = 0 in Ω s .

(3.

2)

The equations are completed with the boundary conditions depending on the considered configuration wether disordered or periodic.

As mentionned in (2.1)-(2.2), the boundary ∂Ω s is divided into disjoint parts ∂Ω s = Γ ext s ∪Γ int denoting the external surface and the internal surface of the solid domain. The boundary Γ ext s corresponds to the part where the displacement u is prescribed and the boundary Γ int corresponds to the part where the traction σ • N is imposed (no-stress) where N stands for the outward unit normal to Ω s on ∂Ω s .

In order to determine the overall elastic response of the RVE, we apply an uniform loading E corresponding to uniform strain boundary conditions [START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF]. In the disordered framework, we enforce

u = E • x on Γ ext s , (3.3)
σ • N = 0 on Γ int , (3.4)
where E is a symmetric second order tensor.

In the periodic framework [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Sanchez-Palencia | Homogenization Techniques for Composites Media[END_REF], the problem is set on the unit cell of the periodic structure, herein the Kelvin cell. The displacement is seeked under the form of a sum of an affine field E • x and a periodic field ũ. Anew, the uniform tensor E stands for the macroscopic strain tensor in the homogenization approach. The corresponding boundary conditions read

u = E • x + ũ on Γ ext s , (3.5
) σ • N is anti-periodic on Γ ext s , (3.6)
We denote by σ and ε the local microscopic stress and strain fields associated to the prescribed macroscopic strain tensor E for both the disordered and periodic boundary value problems.

The problem is then written in a variational form. In the disordered framework, we set

V u = v ∈ H 1 (Ω s ) 3 v = E • x on Γ ext s . (3.7)
In the periodic framework, we set

V per u = v ∈ H 1 (Ω s ) 3 v = E • x + ṽ on Γ ext s , ṽ is periodic . (3.8)
Taking into account the boundary conditions, the variational formulation is written in both configuration:

find the displacement u ∈ V u (resp. u ∈ V per u ) such that Ω s 2µε(u) : ε(v) dV + Ω s λ div(u) div(v) dV = 0, for all v ∈ V u (resp. for all v ∈ V per u ).
(3.9)

The discretization of the problem is described in the Section 5.

The thermal conductivity problem

As for the elasticity problem set in subsection 3.1, the problem is set on the solid and conductive Plateau border Ω s , and the complementary domain Ω p is insulating. Heat transfer through the solid phase is described by the classical linear Fourier's law. The constitutive equation reads:

q(x) = -Λ∇θ(x) in Ω s , (3.10)
where θ(x) denotes the temperature at point x, q(x) is the heat flux and Λ is the conductivity coefficient. As we consider thermal equilibrium states of the material with no heat source, the thermal balance law is expressed as

div q(x) = 0 in Ω s . (3.11)
The equations are completed with the boundary conditions depending on the considered configuration wether disordered or periodic.

In order to implement the homogenization approach, we apply an homogeneous loading A standing for a macroscopic temperature gradient and corresponding to homogeneous boundary conditions.

In the disordered framework, we enforce

θ = A • x on Γ ext s , (3.12
)

q • N = 0 on Γ int , (3.13)
where A is a vector. These boundary conditions model a thermal conductivity behavior with insulation (no flux) at the pore-solid interface. In the periodic configuration, the temperature is seeked under the form of a sum of an affine function and a periodic function θ. The boundary conditions read

θ = A • x + θ on Γ ext s , (3.14)
q • N = 0 on Γ int , (3.15)
The solution represents the local flux q and the local temperature θ at the microscopic scale, induced by the macroscopic temperature gradient A.

The problem is then written in a variational form. In the disordered framework, we set

V θ = τ ∈ H 1 (Ω s ) τ = A • x on Γ ext s . (3.16)
In the periodic framework, we set

V per θ = τ ∈ H 1 (Ω s ) τ = A • x + τ on Γ ext s , τ is periodic . (3.17)
Taking into account the boundary conditions, the variational formulation is written for both problems: find the temperature θ ∈ V θ (resp. θ ∈ V per θ ) such that

Ω s Λ∇θ • ∇τ dV = 0, for all τ ∈ V (resp. for all τ ∈ V per θ ).
The discretization of the problem is described in the Section 5.

The Stokes problem

At the microscopic scale, the creeping flow of a viscous fluid is governed by the Stokes equations. The domain Ω p ⊂ R 3 standing for the porous part is saturated by an incompressible Newtonian fluid.

The constitutive equation reads

σ = 2µd(u) -pδ in Ω p , (3.18)
where u denotes the velocity, p is the pressure field, d is the strain rate tensor and µ > 0 is the dynamic viscosity. We consider the inertialess flow with no body force of the fluid through the porous space. Putting the constitutive equation (3.18) into the momentum balance equation yields

-µ div d(u) + ∇p = 0 in Ω p . (3.19)
Moreover, the fluid is assumed to be incompressible, thus the velocity satisfies the incompressibility condition

div u = 0 in Ω p . (3.20)
The equations are completed with the boundary conditions. For fluid flow through the solid foam, we restrict to disordered material because we do not think it is worth studying again the relationship between RVE size and precision of the overall property estimate.

As mentionned in (2.1)-(2.2), the boundary of the domain is disjoint into ∂Ω p = Γ int ∪ Γ ext p , where Γ int stands for the solid-fluid interface and Γ ext p is the external surface of the porous domain. We assume a no-slip boundary condition at the solid-fluid interface u = 0 on Γ int .

(3.21)

The homogenization approach is implemented by applying a macroscopic pressure gradient α α α on the RVE. We enforce homogeneous macroscopic pressure gradient loading α α α by imposing

σ • N = -(α α α • x)N on Γ ext p , (3.22)
where α α α is a vector. Finally, we write the variational formulation of the problem. The Stokes problem is formulated in a mixed variational form, the velocity u ∈ V s and the pressure p ∈ M s are approximated simultaneously. We set

V s = v ∈ H 1 (Ω p ) 3 v = 0 on Γ int .
Taking into account the boundary conditions, the variational formulation reads: find (u, p) ∈ V s × M s such that

Ω p (2µd(u) : d(v) -p div v) dV + Γ ext p (α α α • x)(N • v) dS = 0, (3.23
)

Ω p q div u dV = 0, for all (v, q) ∈ V s × M s .

Homogenization of the local fields

In this section, we describe the homogenization of the local fields for the linear elasticity problem and the thermal conductivity (on the solid domain), and for the permeability (on the porous domain). In the homogenization step, the macroscopic properties are derived from the volume average of the solutions of the localization problems set on the RVE. In the framework of homogenization of porous media, the solutions for the localization problems are not defined over the whole RVE: displacement and temperature fields are defined over the solid phase only, while velocity and pressure fields are defined over the porous space. It is therefore necessary to define suitable extension of these fields over the whole RVE domain in order to define averaged quantities. In this context, the derived macroscopic properties stand for the effective properties of the material.

Homogenization of the Hooke's law

The homogenization of the linear elasticity problem [START_REF] Hashin | On some variational principles in anisotropic and nonhomogeneous elasticity[END_REF] is addressed in this subsection. Let (σ, ε) defined over Ω s , (σ, ε) denotes the solution of the elasticity problem (3.1)-(3.4) with loading parameter E, standing for the macroscopic strain on the RVE. The displacement u of the elasticity problem (3.1)-(3.4) is extended over the pore space by any continuous field complying with the boundary condition u = E • x over Γ ext p . In this framework, the local stress field ε

satisfies ε = 1 |Ω| Ω ε dV = E,
whatever the extension of the displacement field. By extending the Cauchy stress tensor σ by zero over the porous space, the macroscopic stress field is defined as the average of the microscopic stress field over the RVE [START_REF] Dormieux | Microporomechanics[END_REF][START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF]]

Σ = σ = 1 |Ω| Ω σ dV. (4.1)
Thanks to the fact that the extended Cauchy stress tensor complies with a nobody force balance equation and a no force condition on internal surface Γ int , Σ reads

Σ = 1 |Ω| ∂Ω (σ • N) ⊗ x dS = 1 |Ω| Γ ext s (σ • N) ⊗ x dS . (4.2)
The equivalent homogeneous law states that the macroscopic stress Σ is related to the macroscopic strain E by

Σ = C hom : E, ( 4.3)
where C hom is the fourth order tensor of macroscopic elastic moduli. The linearity of the macroscopic state law (4.3) arises from the linearity of the elastic problem (3.1)-(3.4): σ and ε linearly depend on E. The elastic macroscopic tensor stands for the effective elastic properties of the material. Taking into account symmetry, the macroscopic tensor C hom is defined by 21 independent homogenized coefficients. For each loading E, the solution of the problem (3.9) together with the equations (4.1)-( 4.3) allows to compute six independent coefficients of C hom . Doing so, the whole 21 coefficients of C hom are obtained from the solutions of the elastic problem solved for six independent values of E.

In case of a macroscopic isotropic material, the tensor C hom depends only on 2 independent coefficients (λ hom , µ hom ). Therefore, computing the homogenized coefficients of a disordered RVE allows to check if the equivalent homogeneous material is isotropic.

A classical question of homogenization is whether the macroscopic property (i.e. the elastic tensor C hom ) depends on the particular boundary conditions (here the uniform strain boundary condition 3.4) used to define a well-posed localization problem. For non porous composite materials the classical Hill-Mandel theorem [START_REF] Hill | The essential structure of constitutive laws for metal composites and polycristals[END_REF][START_REF] Mandel | Plasticit classqiue et viscoplasticit[END_REF] states that the overall stiffness tensor C hom computed from the solution of a uniform strain boundary value problem and the overall compliance tensor S hom computed from the solution of a uniform stress boundary problem set on the same REV are the inverse of one another up to the third power of the ratio of the typical size of heterogeneity to the RVE's length. Then when the scale separation condition is fullfilled, the overall stiffness tensor does not depend on the particular conditions used to set the boundary value localization problem. When dealing with porous material it is no more possible to apply a uniform stress boundary condition over the RVE because the porous space intersects the RVE's boundary (see Figures 3 and4). Two solutions allow to show that Hill-Mandel theorem is still valid in this situation.

First it is possible to consider the whole RVE (solid + porous domain) as a heterogeneous linear elastic structure in which the elastic moduli µ micro (x) and λ micro (x) depend on the position vector with, µ micro (x) = µ, λ micro (x) = λ in Ω s and µ micro (x) = µ p , λ micro (x) = λ p in Ω p . Of course Hill-Mandel lemma applies to this classical heterogeneous elastic material whatever the values of the constituents' elastic moduli. As the porous material case is recovered by letting the elastic moduli µ p and λ p tend towards zero [START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF], one can readily conclude to the validity of Hill-Mandel lemma for elastic porous material.

Another way to demonstrate that Hill-Mandel lemma applies to linear elastic porous material is to compare the overall elastic moduli computed from the solution of uniform strain boundary value problem to the same quantity computed from another boundary value problem set on the same REV. For overcoming the impossibility of applying a uniform stress loading on a porous REV we turn to a periodic porous material to address this question. In this situation it is possible to define a boundary value problem set on the unit cell in which the boundary conditions are periodic. Furthermore, solving the uniform strain boundary condition problems for RVE of increasing size allows to evaluate the influence of the RVE size on the accuracy of the computed overall properties. Details of this study are described in section 5.2.1 for both the elasticity problem and the thermal conductivity problem.

Homogenization of the Fourier's law

Herein, we proceed to the homogenization step to complete homogenization of the Fourier's law over the Plateau borders network. We consider the energy flux field q solution of the thermal conductivity problem (3.10)-(3.11)-(3.12)-(3.13) set on Ω s , induced by the loading parameter A. The flux q is extended by zero to the pore space. The macroscopic flux Q is defined as the average of its microscopic counterpart q:

Q =< q >= 1 |Ω| Γ ext s (q • N) x dS = 1 |Ω| Ω q dV.
At the macroscopic scale, the thermal conductivity is described by the equivalent Fourier's law

Q = -Λ hom • A, (4.4)
where Λ hom is the macroscopic conductivity tensor of order two and A represents the macroscopic temperature gradient. We recall that Λ hom is symmetric definite positive tensor, which can be computed from the three solutions of thermal conductivity problem [START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF]. The homogenization of the Stokes equations yields the Darcy's law at the macroscopic scale, inferred as an homogenized limit model [START_REF] Allaire | homogenization of the Stokes flow in a connected porous medium[END_REF][START_REF] Tartar | Incompressible fluid flow in a porous medium -convergence of the homogenization process[END_REF].

We consider the solution (u, p) of the Stokes problem (3.19)-(3.21) induced by the uniform loading condition (3.22) defined by a given vector α α α. The velocity u is extended by zero over the solid space and the continuity is ensured by (3.21). In this framework, the uniform loading α α α represents the macroscopic pressure gradient equal to the average of the microscopic pressure gradient over Ω. As a consequence, at the macroscopic scale Darcy's law linearly relates the pressure gradient α α α to the fluid flux U by

U = -κ hom α α α, (4.5)
where κ hom denotes the macroscopic permeability tensor of order two. Furthermore, the macroscopic flux U is related to the microscopic velocity u through

U = 1 |Ω| Ω u dV = u . (4.6)
Using (4.5) and (4.6), the macroscopic permeability tensor κ hom is entirely determined from the solutions of the Stokes problem defined on Ω p induced by three independant loading vectors α α α. In the isotropic case, κ hom = k hom δ.

For isotropic media, we introduce the intrinsic permeability k ′hom = k hom µ which depends only on the geometrical features of the porous space Ω p . Without loss of generality, and for the sake of clarity, we solve the Stokes problem (3.19) to (3.22) for µ = 1Pa.s which allows to directly compute k ′hom = k hom using (4.5) and (4.6). Thus, the effective isotropic macrocopic permeability κ hom = k hom δ = µk ′hom δ of the heterogeneous solid foam is completely determined.

Numerical study

This section is devoted to the numerical resolution of the homogenization problems. We first present the discretization of the localization problems for elasticity, thermal conductivity and permeability as stated in Section 3. Then, homogenized coefficients are simply determined through averages of the solution of localization problems as explained in Section 4.

We perform a numerical study and present our result for elasticity, thermal conductivity and permeability. We begin by investigating a periodic configuration, which helps shed light on the influence of the size of the RVE on the homogenized coefficients. Afterwards, we investigate the disordered configuration. The accuracy of our results in the disordered case is assessed through the numerical results in the periodic configuration.

Discretization of the problem

In this section, we precise the discretization of the problems presented in subsections 3.1, 3.2 and 3.3. The domain Ω denotes the Plateau border Ω s for both the elasticity and thermal problems, and the porous part Ω p for the Stokes problem. We introduce a regular family of triangulations {T h } h>0 of the domain Ω. More precisely, the domain Ω is discretized by tetrahedral finite elements T h where h > 0 is the discretization parameter standing for the mean edge length. The approximation of the domain, designated by Ωh , is thus the union of T h . The space discretization is based on finite element.

First, we consider the elasticity problem. To approximate the displacement, we consider the finite element space

V uh ⊂ V u defined as V uh (v d ) = v h ∈ (C 0 ( Ω)) 3 ∀T ∈ Ωh , v h | T ∈ P 1 , v h | Γ ext s,h = v d
, where C 0 ( Ω) stands for the space of continuous functions on Ω and P 1 represents the space of polynomial functions of degree 1. Thus, this P 1 setting represents continuous and piecewise linear finite elements and we strongly impose the velocity boundary condition on the boundary Γ ext s,h with a given velocity v d . The approximate velocity field is denoted u h ∈ V uh (E • x). The approached problem associated with the variational formulation (3.9) is derived from the Galerkin method and reads: find

u h ∈ V uh (E • x) such that Ωh 2µε(u h ) : ε(v h ) dV + Ωh λ div u h div v h dV = 0, for all v h ∈ V uh (E • x).
(5.1)

Different boundary conditions must be taken into account for periodic porous media (Eq. (3.5) instead of Eq. (3.3)). In the general case, the definition of the finite element space should be modified to account for these differences. However, as we consider only unit cells that possess geometrical symmetries allowing to replace boundary condition (3.5) by classical Dirichlet conditions (3.3), we do not need to address this problem in this paper. This remark also applies to the thermal conductivity problems.

Then, we consider the thermal conductivity problem. The problem is discretized using the Galerkin method together with the finite element method, as for the elasticity problem, leading to the discrete variational formulation:

find θ h ∈ V θh (A • x) such that Ωh Λ∇θ h • ∇τ h dV = 0, for all τ h ∈ V θh (A • x). (5.2)
where V θh is an approximation of V θ defined in (3.16). Finally, we consider the Stokes problem. To approximate the velocity and the pressure fields, we adopt mixed finite element spaces V sh ⊂ V s and M sh ⊂ M s . As explained above for the elasticity problem, the finite element spaces are spanned by continuous piecewise functions defined on Ωh . More precisely, we consider the P 2 /P 1 setting with tetrahedral Taylor-Hood finite element [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms[END_REF][START_REF] Hood | A numerical solution of the Navier-Stokes equations using the finite element technique[END_REF], that is, continuous piecewise quadratic finite elements for the velocity and continuous piecewise linear finite elements for the pressure. The approximate velocity and pressure fields are denoted by u h and p h . The approached problem associated with the variational formulation (3.23) is written: find

(u h , p h ) ∈ V sh × M sh such that Ωh (2µd(u h ) : d(v h ) -p h div v h ) dV + Γext p ,h (α α α • x)(N • v h ) dS = 0, for all (v h , q h ) ∈ V sh × M sh ,
Ωh q h div u h dV = 0, for all q h ∈ M sh .

(5.3)

Note that it would have been enough to add the periodicity condition respectively for the displacement, the temperature, the velocity and the pressure to the sets V uh , V θh , V sh and M sh to define the discretization of periodic problems.

Numerical results

Thereafter, we precise the numerical parameters of the simulation, then we present and discuss our numerical results. As we provide only dimensionless results in the sequel, the adopted numerical values for lengths and material properties are of no importance. The only relevant parameters are the porosity for the three overall properties (elasticity, thermal conductivity and intrinsic permeability), and the Poisson ratio for the elasticity. For the elasticity problem, the material parameters modeling the elastic behavior are the Young's modulus E = 1000 Pa and the Poisson's ratio ν = 0.3. The corresponding Lamé coefficients are commonly expressed as

λ = Eν (1 + ν)(1 -2ν) , µ = E 2(1 + ν) .
We implement the homogenization approach by considering six loadings: an uniaxial deformation in each direction of space and a pure shear in each couple of directions of space. More specifically, in case of uniaxial deformation the macroscopic loading is set to E ii = e i ⊗ e i with i ∈ {1, 2, 3} (without summation over the repeated indice i) where e i stands for the unit normal vector of the canonical Euclidean basis of R 3 . In case of shearing, the macroscopic loading is set to E i j = e i ⊗ e j + e j ⊗ e i with i, j ∈ {1, 2, 3} and i j. For the thermal conductivity problem, the conductivity coefficient is set to Λ = 1000 W.m -1 .K -1 . In the context of the homogenization approach, we apply a loading in each direction of space A = e i with i ∈ {1, 2, 3}.

For the Stokes problem modeling the fluid flow through the porous space, the dynamic viscosity µ is set to 1 Pa.s for convenience, and the macroscopic loading is set to α α α = e i with i ∈ {1, 2, 3}.

Periodic configuration

In the following, we place ourselves in the periodic configuration of a Kelvin structure [START_REF] Kelvin | The Kelvin problem : foam structures of minimal surface area[END_REF]. For an overall overview of the mechanical behavior of a porous media in the periodic framework, the reader is referred to [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Sanchez-Palencia | Homogenization Techniques for Composites Media[END_REF]. As already mentioned above, our purpose is to evaluate the influence of the RVE size on the accuracy of the homogenized coefficients computed from the uniform boundary loading (Eq. (3.3), (3.12) and (3.22)). Indeed, in the context of the homogenization approach, a condition upon the RVE size is necessary to ensure the convergence of the homogenized coefficients towards the effective properties of the material. Denoting by d the size of the unit cell, by ℓ the size of the RVE and by L the size of the global structure, a required condition (for the validity of the homogenization approach) is d ≪ ℓ ≪ L. More precisely, the ratio δ of the size of the unit cell d (which is close to the bubble radius, see figure 5) to the size of the RVE ℓ has to tend towards 0, that is the size of the unit cell must be negligible compared to the size of the RVE. The condition of separation of scale d ≪ ℓ ≪ L is assumed all along this work, but the ratio δ = d/ℓ can be significantly larger than 0 in the numerical simulations. Therefore, in our numerical study we slightly increase the RVE size for the sake of decreasing the ratio and in turn improve the accuracy of the results.

Furthermore, we emphasize that from the numerical point of view, there are two differents convergences. The first convergence relies on the homogenization approach. Let us slightly clarify: the increase of the size of the RVE leads to δ → 0 + and lim

δ→0 + C hom δ = C hom = C per
where C hom δ is the macroscopic tensor determined on a RVE of size ℓ, made of unit cells of size d = δℓ, subjected to uniform loading (3.3), and the macroscopic tensor C hom characterizes the overall elastic behavior of the structure's constitutive material. In this context, C hom can be determined from the solution of the homogeneous boundary value problem set on the RVE, as δ → 0 + (depicted in the middle of figure 7). Besides, the homogenized macroscopic tensor rising from the solution of a periodic boundary value problem set on the unit cell, is denoted by C per . It has been proved that C per provides the homogenized macroscopic tensor C hom of a periodic RVE subjected to uniform loading when scale separation is achieved, (see [START_REF] Sanchez-Palencia | Homogenization Techniques for Composites Media[END_REF][START_REF] Zaoui | Continuum micromechanics: survey[END_REF]). Therefore, C hom can be defined from the solution of the periodic boundary value problem set on the unit cell, and we infer C hom = C per . Similarly, we have lim δ→0 + Λ hom δ = Λ hom = Λ per for the thermal conductivity problem, and lim δ→0 + κ hom δ = κ hom = κ per for the permeability problem. The framework is depicted in figure 7. The second convergence relies on the finite element method. The refinement of the mesh leads to lim h→0 u h = u where u stands for the exact solution of the problem set in the variational form (Section 3) and h is the mesh size for the discretization of the problem (Subsection 5.1). Besides, considering a Kelvin periodic structure, the determination of the exact solution of the boundary value problem set on the Kelvin unit cell is eased by the geometry. Indeed, the Kelvin unit cell has three planes of symmetry, and each symmetry allows to halve the unit cell. Thus, the local problem is reduced to a local problem on one eighth of the Kelvin cell, designated by Ω h K 8 . Considering a RVE made up of (2N p ) × (2N p ) × (2N p ) unit cell (see figure 8) also involves three planes of symmetry. Therefore, the solution of a boundary value problem with uniform loading can be determined from a problem set on one eight of the system. the deviation from the scale separation condition δ ≪ 1.

C hom δ C per d ℓ = d δ L σ = C hom : ε u = E • x + ũ u = E • x Structure unit cell RVE RVE

Elasticity.

We first consider the elasticity problem. We determine the homogenized macroscopic tensor C hom on a unit Kelvin cell, as explained in Section 3.1 with periodic boundary conditions (3.5)- (3.6). We implement the discrete formulation (5.1) and strongly enforced the Dirichlet boundary conditions in the definition of V uh . The resolution is achieved through the CAST3M software [START_REF]Cast3m[END_REF]. We determine the macroscopic elastic tensor C hom of Eq. ( 4.3) by averaging the stress field on the whole domain as explained in the subsection 4.1. In the periodic configuration, the structure doesn't exhibit isotropy, thus the overall macroscopic tensor C hom has six independent components. For the sake of simplicity, we focus on the coefficient C hom xyxy . We plot in figure 9 the ratio

C hom xyxy (Ω h N p )/C hom xyxy (Ω h K 8
) as a function of the number of cells per axe N p . We compare the results obtained on the r 0 mesh (left) to those obtained on the r 1 mesh (right). The plotted red line enables to evaluate the gap between the result on Ω h N p and the exact solution determined on Ω h K 8 (and reached for N p → +∞ at fixed h). The black line indicates 10% of deviation from the exact solution. Moreover, the second axe on top provides the corresponding number of cells in the disordered configuration (it is recalled that with symmetric boundary conditions the problem solved on N p cells per axe provides the solution of a (2N p ) 3 cells RVE). This precision sheds light on the estimation of the influence of the size of the RVE on the homogenized coefficient in the disordered configuration. We do not notice any effect of the foam porosity on the rate of convergence of the overall properties of the foam as the RVE size increases. It does not means that such an effect does not exist but, as we study foams with solid volume fraction ranging from 0.025 to 0.2, it is likely that this effect should be rather small.

The smooth convergence of the estimate for the overall elastic modulus as a function the number of cells is a classical result [START_REF] Michel | Composites à microstructure périodique[END_REF]. Indeed the rate of convergence of the method is linear in RVE size. More precisely, the data depicted in Figure 9 are well fitted by the equation

C hom xyxy (Ω h N p )/C hom xyxy (Ω h K 8 ) = 1 + αN -1 = 1 + α(d/ℓ) -1
with α ≃ 0.7. Then, for N p ≥ 6 the results are about 10% of difference from the converged homogenized solution. This observation provides an estimation of the accuracy of the result in the disordered configuration with N cells. Since the result on N p = 4 cells corresponds to the result on N = (2N p ) 3 = 512 disordered cells, we can consider that the homogenized coefficients obtained on a disordered RVE of N = 512 cells are about 20% different from the effective coefficients (at fixed h). In order to reach about 10% of deviation from the exact solution, one has to consider N = 1728 cells. The same overall trend is observed on the r 1 mesh. These results are in good agreement with the study conducted in [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF] which figures out 20% difference between numerical and theorical shear modulus.

Besides, let us further comment that, the convergence of the homogenization approach is faster on the reference mesh Ω h K 8 than on Ω h N p . This numerical effect is observed in figure 9 since more important variations of C hom xyxy are observed for the r 1 mesh figure (right), whereas in regards to the finite element convergence C xyxy (Ω h,r 0 N p ) > C xyxy (Ω h,r 1 N p ) obviously holds. On the r 1 mesh figure, the result for N p = 10 isn't available due to lack of memory in computing resource.

Thermal conductivity.

As a second step, we consider the thermal conductivity problem. As explained for the elasticity problem, we consider a mesh of the unit Kelvin cell, taking into account one eight of the Kelvin cell Ω h K 8 , and a N p 3 unit cells RVE Ω h N p , with N p ∈ {2, 3, 4, 5, 6, 7, 8, 9, 10}. The solution of reference, related to the mesh Ω h K 8 , is compared to the solution on Ω h N p . We recall that we take into account a foam structure with a solid volume fraction φ = 5%. We determine the homogenized macroscopic tensor Λ hom on a periodic Kelvin structure, as modeled in Section 3.2 with periodic boundary conditions (3.14)- (3.15). We implement the discrete formulation (5.2) and strongly enforced the Dirichlet boundary conditions in the definition of V θh . The resolution is achieved through the CAST3M software. We determine the macroscopic thermal conductivity tensor Λ hom of Eq. (4.4) by averaging the heat flux field on the whole domain as explained in subsection 4.2. The periodic RVE does not exhibit isotropy, thus the discrete macroscopic tensor Λ hom has two independent components. For the sake of simplicity, we focus on the coefficient Λ hom

x . We plot in figure 10 the ratio Λ hom

x (Ω h N p )/Λ hom x (Ω h K 8 ) with respect to the number of cells per axe N p . We compare the result obtained for both the r 0 mesh (left) to those obtained for the r 1 mesh (right). As for the elasticity problem, the plotted red line represents the result when convergence is achieved, and the black line indicates 2% of deviation from the exact solution. The second axe on top provides the corresponding number of cells of the RVE.

: Ratio Λ hom x (Ω h Np )/Λ hom x (Ω h K 8
) as a function of the number of cells N p , for the r 0 mesh (left) and for the r 1 mesh (right). Red: exact solution. Black: 2% deviation from the exact solution.

We observe that from N p = 2 to N p = 10, the difference with the converged homogenized solution (red line) decreases approximately from up to 6% to 1%.

Similarly to what is observed for the linear elasticity, the rate of convergence for thermal conductivity is linear in RVE size. The method error depicted in Figure 10 is accurately predicted by the equation

Λ hom x (Ω h N p )/Λ hom x (Ω h K 8 ) = 1 + βN -1 = 1 + β(d/ℓ) -1 with β ≃ 0.1.
Comparing to the elasticity problem, this is the reflect that the convergence is reached more rapidly with the conductivity problem. Since the results computed using symmetries and N p = 4 correspond to the results for a RVE containing N = 512 cells, it may be concluded that thermal conductivity computed with a disordered RVE of N = 512 cells may be regarded as an estimate of the homogenized thermal conductivity with an error of nearly 3%.

In regards to the finite element convergence, as we observe the same overall trend for both the r 0 and the r 1 mesh, we infer that the convergence on the reference mesh Ω h K 8 is similar to the convergence on the N p 3 cells mesh Ω h N p . Remark. We do not provide the results for the Stokes problem due to a lack of memory of the computing resource. However we solve the Stokes problem for small values of N p and periodic unit cell. We observe exactly the same trends as for the elastic and the thermic problems (estimates for the overall permeability of periodic material are decreasing functions of the number of cells tending towards the homogeneous permeability computed on the periodic unit cell). Consequently, we consider in the sequel that the size of disordered solid foam RVE is large enough so that our estimates are accurate to account for the effect of solid volume fraction.

If these trends are extrapolated to disordered materials, they suggest that a RVE should contain at least 1728 bubbles so that its overall elastic properties can be estimated with an error less than 10%, while the same RVE allows to compute the overall thermal conductivity with less than 2% error (see figure 10). In the sequel, the larger RVE we consider contains 512 bubbles (it was not possible to solve larger problems due to lack of power and memory of computer resource). Therefore, our computations provide estimates with a minimum percent error of 15% for the overall elastic modulus and 3% for the overall thermal conductivity. This must be borne in mind when comparing the overall properties of two samples with different morphological properties (mainly porosities in this work): differences between two macroscopic properties is significative only if they are greater than error margin which is the case for our results.

Disordered configuration

The RVE occupies a sphere of diameter ℓ, with ℓ ∈ [START_REF] Allaire | homogenization of the Stokes flow in a connected porous medium[END_REF][START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF] depending on the number of bubbles in the RVE. We consider an unstructured mesh of the spherical domain, containing N bubbles with N ∈ {32, 64, 128, 256, 512}. Each bubble is constrained to have equal volume v B . In the sequel, we considered RVE with solid volume fraction φ ∈ {2.5, 5, 10, 15, 20}, thus the RVE's volume is Nv B /(1 -φ).

We set the value of the mesh size h with respect to the computing resource. Then the choice of RVE size represents a compromise between computational practicability and precision of the computed quantities. Indeed, the mesh size h directly influences the number of elements N tet (tetrahedra) of the volume mesh. The mesh size is approximately set to h = 10 -1 by averaging the edge length of the surface mesh. We designate by r 0 the basis mesh with zeroth level of refinement and by r 2 the twice refined mesh.

Remark. On the domain Ω p representing the porous part, the resolution of the problem is too demanding in memory computing ressource to allow a refinement. As a consequence, we study the Stokes problem on the r 0 mesh only.

The number of tetrahedra with h = 10 -1 , is reported in Table 1 for φ = 2, 5%.

φ = 2.5% N = 32 N = 64 N = 128 N = 256 N = 512 Ω s :

Elasticity.

We implement the discret formulation (5.1). The Dirichlet boundary conditions are strongly enforced in the definition of V uh . The resolution is achieved through the FreeFem++ software [START_REF] Hecht | New development in FreeFem++[END_REF], the Conjugate Gradient method is employed to solve the system of linear equations. However, for larger systems such as those related to the r 2 mesh or the N = 512 cells mesh, the method doesn't converge since the bad matrix conditioning leads to a divergence of the method when the number of degrees of freedom increases. In this case, we adopt the GMRES (Generalized Minimum RESidual) method [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], an iterative Krylov subspace method [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF][START_REF] Freund | Iterative solution of linear systems[END_REF] for nonsymmetric and large systems, which is more general and more robust than the Conjugate Gradient algorithm, whereas it is more expensive in computing time due to the iterativity.

We compute the macroscopic elastic tensor C hom of equation ( 4.3) by averaging the stress field on the whole domain occupied by the RVE as explained in subsection 4.1. We have found that foam samples we have designed exhibit isotropy (taking into account numerical errors). Therefore, all the components of the elastic tensor C hom are function of only two coefficients: the bulk modulus K hom and the shear modulus µ hom . Dimensionless macroscopic shear modulus µ hom /µ and bulk modulus K hom /K with µ (resp. K) the solid phase shear modulus (resp. bulk modulus) as a function of the solid volume fraction are depicted in figures 11 and 12. The results for several RVE sizes (i.e. numbers of bubbles) computed with r 0 and r 2 meshes are presented. The macroscopic elastic properties are increasing functions of the volume fraction. The size of the RVE affects the results mainly for the shear modulus. Regarding the bulk modulus computed with the r 2 mesh, the converged homogenized solution is reached even with N = 32 cells, whereas for the shear modulus, the converged homogenized solution is not perfectly reached even with N = 512 cells. In addition, it is worth noting that mesh refinement leads to a decrease in overall shear modulus estimate. This trend is consistent with the minimizing principle associated with the variational formulation (3.9) of the homogenization problem [START_REF] Dormieux | Microporomechanics[END_REF]: for all u ∈ V u (E • x)

µ hom E : E - 1 3 (trE) 2 + K hom 2 (trE) 2 ≤ 1 |Ω| Ω s µε : ε + λ 2 (trε) 2 dV.
At first sight, V uh can be considered as subset of V u , which means that estimates of µ hom and K hom are upper bounds of exact solutions. As the velocity field V uh (r 0 ) is a subset of V uh (r 2 ) the r 0 -estimate for µ hom is larger than the r 2 -estimate.

Besides, we manage to improve the accuracy of the results by increasing the order of interpolation, using quadratic elements. In figures 13 and 14 we represent the results of a resolution with P 2 element on the r 0 mesh (left) and on the r 2 mesh (right). A closer look of results obtained for the overall shear modulus using different mesh refinements and degrees of interpolation reveals that using r 0 -mesh with P 2 element (figure 14 left) allows to obtain accurate estimates since refining the mesh two times (i.e. multiplying by about 40 the number of elements) does not substantially change the result (compare the N = 32, N = 64 and N = 128 shear modulus estimates depicted in figure 14 for r 0 and r 2 -mesh with same degree of interpolation). The same comment applies to the overall bulk modulus estimates (compare figure 13 left to figure 13 right). Therefore, we consider r 0 -meshes with quadratic interpolation yield sufficiently accurate estimates, especially given that estimates for the homogenized coefficients computed using a 512-bubbles RVE can be underestimated by about 20% (subsection 5.2.1).

In order to assess our results, we compare them to literature data. More precisely, we consider the Hashin-Shtrikman (HS) [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF] upper bounds for the Young modulus

E hom HS E = 2φ(7 -5ν) φ(15ν 2 + 2ν -13) -15ν 2 -12ν + 27 , (5.4)
and the shear modulus

µ hom HS E = φ(7 -5ν) 2(1 + ν)(10νφ -15ν -8φ + 15) . (5.5)
It is reminded that for porous media i) Hashin-Shtrikman lower bounds for E hom and µ hom are equal to zero and ii) Hashin-Shtrikman upper bounds and Mori-Takana estimates are the same [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF].

We also consider the Gibson-Ashby estimates for the Young modulus and the shear modulus [20]

E hom As E = C 1 φ 2 with C 1 ≈ 1, (5.6)
and

µ hom As E = C 2 φ 2 with C 2 ≈ 3 8 .
(5.7)

Comparing our estimates for the overall elastic moduli of solid foam to Gibson-Ashby estimates is a suitable way to check their ability to accurately fit experimental data since Gibson and Ashby have shown in [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF] (figure 5.10 page 192) that equations (5.6) and (5.7) reasonbly fit experimental data even if they are widely dispersed (for a given value of the porosity, experimental data collected by Gibson and Ashby can be distributed over a range of ±50% of the average value).

We compare the Hashin-Shtrikman and Gibson-Ashby estimates to the estimates of the dimensionless elastic modulus computed for a 512-bubbles RVE of foam on a r 0 -mesh with P 2 element. In figure 15 (left), we plot the dimensionless elastic modulus E hom HS /E and E hom As /E as a function of the solid volume fraction. In figure 15 (right), we plot the ratio µ hom HS /E and µ hom As /E as a function of the solid volume fraction. It is worth noting that change of solid volume fraction induces variation of the macroscopic elastic properties large enough to be significant with respect to the estimation error we discuss above. As a consequence, we are confident that our results allow to accurately estimate the overall properties of real monodisperse isotropic open cells solid foams. First, it can be seen in figure 15 that our estimates for the overall Young modulus and the overall shear modulus are smaller than the Hashin-Shtrikman upper bounds (5.4) and (5.5) and larger than the Gibson-Ashby estimates (5.6) and (5.7). Theorically, overall properties of the foam depend on the value of the Poisson ratio of the solid phase, as shown by Eq. (5.4) and (5.5) for the Hashin-Shtrikman upper bounds. In figure 15 we plot the Hashin-Shtrikman bounds for ν = 0, ν = 0.3 and ν = 0.5. The three curves overlap in the (φ -E hom ) diagram (figure 15 left) which means that the solid phase Poisson ratio has a weak influence on the Young modulus value. For the overall shear modulus, it can be seen in figure 15 right that Poisson ratio of the solid phase has a more pronounced effect on the overall shear modulus even if it remains small enough in comparison with the scatter of experimental data [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF].

We have fitted our numerical data to a power law as done by Gibson and Ashby (Eq. (5.6) and (5.7)). We have inferred that numerical data depicted in figure 15 are well described by 1.45 (5.8) and µ num E = 0.3φ 1.5 .

E num E = 0.7φ
(5.9)

When plotted in a solid volume fraction vs dimensionless overall elastic modulus diagram, estimates (5.8) and (5.9) fit accurately enough our data so that one does not see any gap between the points representing the data and the lines representing the power law in figures 15. The power laws (5.8) and (5.9) are in good agreement with results obtained in [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF] which state power laws with exponent from n = 1 to n = 3.15 depending on the samples structure. The closest structures of [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF] to our foam structure lead to exponent between 1 and 2 as in (5.8) and (5.9). To conclude on overall elastic modulus, it must be noticed that our numerical results do not confirm the Gibson and Ashby prediction ν AS = 0.3 for the overall Poisson ratio of a foam [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. We find values of ν num between 0.26 and 0.35 for monodisperse solid foam computed with ν = 0.3. Regarding the scatter of experimental data, this result is consistent with the Gibson and Ashby estimate. Lastly, let us mention that the estimate ν = 0.3 is not consistent with the relationship (5.6) and (5.7) since putting (5.6) and (5.7) in the equation ν AS = E AS /(2µ AS ) -1 yields ν AS = 1/3 = 0.33.

Thermal conductivity.

Rigid foams are mainly used for thermal insulation since they have low thermal conductivity. Heat flow through a rigid foam have four contributions: conduction through the solid, conduction through the porous space (filled by a gas), convection with the connected pores and radiation through the solid faces of the cells. Relative size of these contributions varies as a function of the properties of the foam's components (solid and gas) and of the cell size. Modeling the thermal properties of foam from a description of the four contributions is out of the scope of this paper. Nevertheless, it is interesting to study thermal conduction through the solid since this contribution becomes important for larger values of the solid volume fraction (i.e. φ > 5%) [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF].

We implement the discrete formulation (5.2) with the Dirichlet boundary conditions strongly enforced in the definition of V θh . The resolution is achieved through the FreeFem++ langage and the numerical method is wether the Conjugate Gradient or the GMRES method for larger systems such as on the r 2 mesh or on the N = 512 cells mesh.

We compute the macroscopic thermal conductivity as explicited in subsection 4.2. As the structure exhibits isotropy, the discrete macroscopic thermal conductivity tensor Λ Λ Λ h is reduced to Λ Λ Λ h = Λ h δ. We plot in the figure 16 the dimensionless macroscopic thermal conductivity as a function of the solid volume fraction computed with both the r 0 mesh and the r 2 mesh. As expected, macroscopic thermal conductivity increases with the solid volume fraction. The size of the RVE slightly influences the results since we obtain almost superimposed results with respect to the number of cells N. Therefore, it may be concluded that the size of foam samples we adopt in this work is large enough to be considered as representative of the overall material for solid thermal conductivity. This result confirms what we claim in the section 5.2.1: estimates computed with a N = 512 RVE should be very accurate (difference with the homogenized value less than 3%). On the other hand, it is clear in figure 16 that finite element convergence is not reach with the r 0 meshes since smaller estimates for the overall thermal conductivity are obtained with r 2 meshes.

In order to enhance the result, we increase the order of interpolation and perform a resolution with P 2 elements. We represent in figure 17 the results obtained with P 2 element on both the r 0 and the r 2 mesh. We observe that the resolution with P 2 elements on the r 0 mesh provides the same accuracy as on the r 2 mesh with P 1 elements (figure 16). Furthermore, the resolution with P 2 elements on the r 2 mesh provides the same results as the resolution with P 2 element on the r 0 mesh. In other words, for the thermal conductivity problem, with a r 2 refinement and a P 1 resolution the convergence is reached. By combining the results obtained in subsection 5.2.1 for the convergence of the estimate for the overall thermal conductivity of a RVE made of a periodic material, and the results depicted in figures 16 and 17, we can conclude that our computations allow to estimate macroscopic thermal conductivity of a disordered monodisperse foam with an error less than 3%.

In figure 18, we compare the estimate we obtain to the Hashin-Shtrikman upper bound for the overall thermal conductivity of a porous material when only conduction through the solid phase is taken into account. It is reminded that the Hashin-Shtrikman upper bounds reads It is readily seen that our estimates are always below the Hashin-Shtrikman upper bound. As we did not find an analytical estimate for the solid thermal conductivity of a monodisperse isotropic foam in the literature, we fit our numerical data to a power law. In the range of solid volume fraction we studied, solid thermal conductivity can be estimated by

Λ hom HS Λ = 2φ 3 -φ .
Λ hom Λ = φ 1.2 .
As for elasticity, when plotted in figure 18, this closed form equation almost perfectly fits the numerical results in the range φ ∈ [0.025, 0.2].

Permeability.

In this last section, we present the results we obtain for the intrinsic permeability of a monodisperse disordered solid foam. We implement the discret formulation (5.3) with the dirichlet boudary conditions strongly enforced in the definition of V p h . The resolution is achieved through the FreeFem++ langage, the formulation is a MPI (Message Passing Interface) formulation with MPI resolution on 4 processors. The numerical method implemented to solve the linear system is the Mumps (MUltifrontal Massively Parallel Solver) [START_REF] Amestoy | MUMPS MUltifrontal Massively Parallel Solver Version 2[END_REF] which is a sparse direct solver. The resolution of the discrete problem is achieved separately for u xh , u y h and u z h , since the cost of solving the linear system is prohibitive due to storage requirements. In addition, we have observed that the method doesn't converge for large number of unknowns. We have managed to proceed at resolution up to 1 220 184 degrees of freedom, corresponding to a 512-bubbles RVE with φ = 20% of solid volume fraction.

We determine the macroscopic permeability induced in (4.5) by averaging the velocity field on the whole domain as explained in subsection 4.3. As the structure exhibits isotropy, the discrete macroscopic permeability tensor κ h is reduced to κ ′hom = k ′hom δ. Moreover, in the monodisperse configuration the average cell volume v B stands for the only characteristic length scale which is adopted herein. In order to assess our results, we consider the Kozeny-Carman estimate of the intrinsic permeability of a porous media [START_REF] Carman | Flow of Gases Through Porous Media[END_REF]

k ′ KC = C ρ 2 (1 -φ) 3 φ (5.10)
where ρ is a characteristic length and C is a constant. In the left figure 19, we represent the ratio k ′hom /v B In the left figure 19, we observe that the macroscopic permeability decreases with the volume fraction. The size of the RVE significantly influences the results which is not surprising in view of the relative small size of the numerical problems we solve. More specifically, we observe that the number of cells has more influence on low solid volume fraction φ = 2, 5% than on high solid volume fraction φ = 20%. This fact relies on the slightly higher number of elements for meshes with φ = 20%. Nevertheless, as the permeability computed for N = 256 and N = 512 are close one to the other, we can consider that a N = 512-bubbles RVE is large enough with respect to the bubble size to compute accurately the overall permeability of the solid foam as a function of solid volume fraction. As ever said, we were not able to increase the order of interpolation of the finite elements we use due to lack of computing power; it is very likely that numerical convergence is not reach for the results depicted in figure 19. In the abscence of more reliable results we use these data for comparison of available results.

On the right figure 19, we compare our numerical data to the Kozeny estimate (5.10) computed for C = 1/350 and C = 1/1500. Even if same trends are predicted by numerical simulation and Eq. (5.10), it is obvious that the Kozeny-Carman estimate does not perfectly fits the numerical data. This is not surprising given the fact that Kozeny-Carman estimate is obtained by modeling the porous space by a network of ducts, which does not correspond to the morphology of disordered solid foam. For now, we do not have better estimate to propose.

Overview and concluding remarks

First of all, the starting point is to generate digital model of monodisperse open cell solid foams. We point out that this is not trivial to carry out since the microstructure of digital foams must reach a stable state of equilibrium in order to model physically suitable microstructure of real foams. Furthermore, physically relevant stable structure minimizes the surface energy of the foam, which entails geometrical as well as topological conditions. These conditions present some difficulties to tackle numerically, mainly owing to the geometrical irregularity of the surface such as corners where Plateau borders join membrane. More precisely, these singularities can lead to onset of surface intersections during the foam-making process, which prevents volumic meshing of the digital foam samples. Hence the optimization algorithm, minimizing the energy of the surface, requires a high regularity of the surface model to produce a physically realistic realization of a foam sample. Moreover, the minimization of surface energy is a constrained minimization problem which is ill-posed and thus can lead to unstabilities. In this context, we manage to entirely automate the foam meshing procedure. This enables us to produce digital samples of monodisperse open cell solid foam with solid fraction ranging from 2.5% to 20%. It is worth noting that our procedure is able to produce polydisperse solid or liquid foams with open or closed cells.

In a second step, we employ these digital foam RVEs to estimate overall properties (namely elasticity, thermal conductivity through the solid phase and permeability) of monodisperse open cell solid foam as a function of solid volume fraction in the framework of homogenization methods for disordered materials. Overall properties of the materials are estimated from the finite element solution of boundary value problems set on the foam RVEs. From the physical point of view, we have found that the overall properties of solid open cell solid foams depend strongly on the solid volume fraction for values between 2, 5% and 20% which is consistent with the literature. Besides, an isotropic behavior is predicted for all the digital samples we have produced, which proves that our numerical procedure really generates disordered material. Our results compare well to available results of the literature. We also propose new closed form estimates for both the elasticity and the thermal conductivity of monodisperse open cell solid foam that allow to accurately predict the overall properties of the material as a function of the solid volume fraction and components properties.

In regards to the homogenization techniques, we have studied the influence of RVEs size and mesh fineness on the estimates accuracy. We have shown that our numerical samples are large enough and the mesh fineness is small enough to evaluate the thermal conductivity and the elasticity with an accuracy of respectively ≃ 3% and ≃ 20%. Due to lack of computing power, we were not able to rigorously evaluate the accuracy of permeability estimates. Even if the imprecision associated with elasticity and permeability estimates may seem quite large, it should be kept in mind that experimental data for foams are rather dispersed. Therefore, we are confident that our estimates allow to account for the effects of solid volume fraction and components properties onto the overall properties of monodisperse open cell solid foams. However, for the sake of improving the results precision, it would be interesting to first improve the accuracy of the numerical approximation and then to increase the size of the RVE, which can be done using a more powerful computer.

As a further perspective, it is first possible to use the tools we have developed in this study to generate digital samples of polydisperse foams with open-cells, closed-cells and varying solid volume fractions. As regards the mechanical properties, it should be interesting to address nonlinear behavior originating from large deformation or inelastic behavior of the solid phase. For thermal conductivity, we believe that the key work task to be undertake is modeling the other modes of heat transfer through the porous foam: conduction in the porous space, convection and radiation. Finally, it should be interesting to study foam permeability as a function of relative proportion of closed cells versus open cells, as it is well known that cells membranes are not all destroyed or all preserved during the solidification of a foam.

Figure 1 :

 1 Figure 1: Voronoï tessellation of a spherical domain with N = 64 cells. Left: before relaxation algorithm. Right: after relaxation algorithm.

Figure 2 :

 2 Figure 2: Wet foam structure with a Plateau border (red pipe) of φ = 5% volume fraction. Left: before optimization. Right: after optimization (no refinement).

Figure 3 :

 3 Figure 3: Meshed microstructure of a disordered monodisperse foam with N = 64 cells and a solid volume fraction φ = 5%. Plateau border network Ω s (left) and porous space Ω p (right). Left figure: the green surface is the external surface Γ ext s and the yellow surface is the internal surface Γ int . Right figure: the green surface is the outer surface Γ ext p and the yellow surface is the internal surface Γ int . Both the solid domain Ω s and the porous domain Ω p are connected spaces since the interfaces between bubbles have been removed.

Figure 4 :

 4 Figure 4: Wet foam structure with φ = 5% volume fraction, and N = 512 cells.

Figure 5 :

 5 Figure 5: Meshed microstructure of the eight of a Kelvin unit cell with solid volume fraction φ = 5% and r 2 refinement (see subsection 5.2.1). Solid space (left) and porous space (right).

Figure 6 :

 6 Figure 6: Meshed microstructure of a Kelvin foam (N p = 2) with volume fraction φ = 5% and r 0 refinement (see subsection 5.2.1).Solid Plateau border network (left) and porous space (right).

  (3.10)-(3.11)-(3.12)-(3.13) induced by three linearly independant vectors A. 4.3. Homogenization of the Stokes equations The homogenization of the Stokes equations is briefly schematized in what follows, the reader is referred to [14] for more details. At the microscopic scale, the motion of a viscous fluid through the porous space is governed by the Stokes equations (3.19)-(3.20)-(3.21)-(3.22) set on Ω p . At the macroscopic scale, the Darcy's law describes the transport of fluid through a porous medium

Figure 7 :

 7 Figure 7: Different scales to consider in homogenization of periodic structures: unit cell, RVE and structure. The macroscopic tensor C per is determined on the unit cell and the macroscopic tensor C hom δ is determined on the RVE.

Figure 9 :

 9 Figure 9: Blue: ratio C hom xyxy(Ω h Np )/C hom xyxy (Ω h K 8) as a function of the number of cells N p , for the r 0 mesh (left) and for the r 1 mesh (right). Red: exact solution. Black: 10% deviation from the exact solution.

Figure 10 :

 10 Figure 10: Blue: Ratio Λ hom x (Ω h Np )/Λ hom x (Ω h K 8 )as a function of the number of cells N p , for the r 0 mesh (left) and for the r 1 mesh (right). Red: exact solution. Black: 2% deviation from the exact solution.

Figure 11 :Figure 12 :

 1112 Figure 11: Macroscopic dimensionless shear modulus µ hom /µ of the solid foam as a function of the solid volume fraction φ. Resultsobtained with the r 0 mesh (no refinement -left) and with the r 2 mesh (2 refinements -right).

Figure 13 :Figure 14 :

 1314 Figure 13: Macroscopic dimensionless µ hom /µ shear modulus of the solid foam as a function of the solid volume fraction φ. With a P 2 resolution on the r 0 -mesh (no refinement -left) and on the r 2 -mesh (2 refinements -right).

Figure 15 :

 15 Figure 15: Dimensionless Young modulus E hom (left)/shear modulus µ hom (right) as a function of the solid volume fraction φ.Numerical estimates compared to Hashin-Shtrikman and Gibson-Ashby estimates. HS bounds and numerical estimates are computed for ν = 0, 0.3 and 0.5.

Figure 16 :

 16 Figure 16: Dimensionless isotropic overall thermal conductivity as a function of the solid volume fraction φ. With the r 0 mesh (no refinement -left) and with the r 2 mesh (2 refinements -right).

Figure 17 :

 17 Figure 17: Dimensionless isotropic thermal conductivity as a function of the solid volume fraction φ. With a P 2 resolution on the r 0 mesh (no refinement -left) and on the r 2 mesh (2 refinements -right).

Figure 18 :

 18 Figure 18: Dimensionless isotropic thermal conductivity, as a function of the solid volume fraction φ, for 512-bubbles RVE of a monodisperse solid foam. Comparison with Hashin-Shtrikman upper bound.

Figure 19 :

 19 Figure 19: left: Discrete permeability k ′hom /v B 2 3 as a function of the volume fraction φ. right: Comparison of the discrete permeability k ′hom /v B 2 3 to the Kozeny-Carman estimate C(1 -φ) 3 /φ with C = 1/350 and C = 1/1500.

Table 1 :

 1 Number of tetrahedra N tet as a function of the number of cells N. Ω s : r 0 -solid phase mesh Ω h s with no refinement. Ω s : r 2 -solid phase mesh with two refinements and Ω p -porous space mesh Ω h p . Volume fraction of the RVE: φ = 2.5%.

	r 0	2 601	4 893	9 669	20 726	43 658
	Ω s : r 2	93 083 189 524 386 862 785 457 1 598 656
	Ω p	12516	27 908	57 739 100 665	199 176

Preprint submitted to ElsevierDecember 11, 2017

Acknowledgment. This work was supported by the ANR ProMAP project, grant ANR-13-RMNP-0003 of the French Agence Nationale de la Recherche.

We are also grateful to Chavance Pierre for assistance to the periodic foams numerical simuulations.

In this framework, the solution of the problem set on Ω h N p with N p = 1 (i.e. one Kelvin unit cell) is entirely determined by solving the problem on Ω h K 8 . We consider a mesh of

, with d = 100 (it is reminded that physical size of the RVE or of the unit cell does not impact the homogenized coefficient, thus it is not necessary to specify the true size of the unit cell here). The RVE size is increased by considering a periodic network of Kelvin cell with N p cells per axe, that is to say a mesh of N p 3 cells, designated by Ω h N p . Due to symmetry along each axe, the result on a mesh made of N p 3 cells provides the effective properties of a structure made of (2N p ) 3 In what follows, we focus on a foam structure with a volume fraction φ = 5%. The others values of φ lead to the same overall behavior. We consider a r 0 mesh referring to the zeroth level of refinement, and a r 1 mesh obtained by halving the r 0 mesh size. The r 0 mesh size is h ≃ 10 -1 .

Of course, we use the same mesh fineness to compute the finite element solution of homogeneous localization problem and periodic localization problem, so that the measured difference between C hom δ and C per only comes from