K Anoukou

R Brenner
email: renald.brenner@sorbonne-universite.fr

F Hong

M Pellerin

K Danas

Random distribution of polydisperse ellipsoidal inclusions and homogenization estimates for porous elastic materials

Keywords: RSA algorithm, ellipsoidal inclusions, FFT-based numerical homogenization, effective elastic properties, porous material

This work proposes an extension of the well-known random sequential adsorption (RSA) method in the context of non-overlapping random mono-and polydisperse ellipsoidal inclusions. The algorithm is general and can deal with inclusions of different size, shape and orientation with or without periodic geometrical constraints. Specifically, polydisperse inclusions, which can be in terms of different size, shape, orientation or even material properties, allow for larger volume fractions without the need of additional changes in the main algorithm. Unit-cell computations are performed by using either the fast Fourier transformed-based numerical scheme (FFT) or the finite element method (FEM) to estimate the effective elastic properties of voided particulate microstructures. We observe that an isotropic overall response is very difficult to obtain for random distributions of spheroidal inclusions with high aspect ratio. In particular, a substantial increase (or decrease) of the aspect ratio of the voids leads to a markedly anisotropic response of the porous material, which is intrinsic of the RSA construction. The numerical estimates are probed by analytical Hashin-Shtrikman-Willis (HSW) estimates and bounds.

Introduction

The effective physical properties of random heterogeneous materials are strongly dependent upon their microstructure. Therefore, the statistically accurate quantitative characterization of the microstructure is of great importance in their modeling. Specifically, knowledge of microstructural statistical information (e.g. volume fractions, shapes, connectivity, spatial distributions and orientations of constituents) is fundamental in understanding and interpreting the microstructure-property relationships. Statistical descriptors such as n-point probability functions, are widely used to quantify mathematically this information in a broad class of random microstructures (see comprehensive reviews in [START_REF] Milton | The Theory of Composites[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]). They are useful in the random generation or stochastic reconstruction of realistic virtual microstructures.

Random generation or stochastic reconstruction of virtual microstructures needs to take into account an infinite amount of microstructural statistical properties. However, in practice, only limited (lower-order) statistical information is available either experimentally or theoretically. The accessibility of the high-level statistical information remains an important focus of research [START_REF] Gillman | Third-order thermo-mechanical properties for packs of platonic solids using statistical micromechanics[END_REF]. Statistical microstructural descriptors (e.g., n-point probability functions, lineal-path function, pore-size distribution functions, etc.), are widely used to quantify mathematically this information in a broad class of random microstructures. A review of several distribution functions can be found in [START_REF] Torquato | Random heterogeneous media: Microstructure and improved bounds on effective properties[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Lu | Linear-path function for random heterogeneous materials[END_REF].

In homogenization based statistical continuum theories [START_REF] Beran | Statistical Continuum Theories[END_REF][START_REF] Willis | Variational and related methods for overall properties of composites[END_REF][START_REF] Kröner | Statistical Continuum Mechanics[END_REF][START_REF] Sen | Effective conductivity of anisotropic two-phase composite media[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Mikdam | Effective conductivity in isotropic heterogeneous media using a strong-contrast statistical continuum theory[END_REF], the microstructural statistics of the representative volume element of random multiphase materials is assumed to be spatially uniform and ergodic. This means that probability functions are insensitive to translations and sample realizations. It follows that, the one-point probability function is simply the volume fraction of the phases. In turn, the phase distribution can be described by two-or higher order probability functions and can be statistically isotropic or statistically anisotropic. A number of two-point correlation function methods have been proposed in the literature. The ellipsoidal symmetry can be traced back to Willis [START_REF] Willis | Bounds and self-consistent estimates for the overall modui of anisotropic composites[END_REF][START_REF] Willis | Variational and related methods for overall properties of composites[END_REF] while Ponte Castañeda and Willis [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] proposed estimates that allow for different inclusion shapes and distribution functions. In the present work, the focus is on "almost isotropic" responses with randomly oriented ellipsoidal voids and use will be made of the Hashin-Shtrikman-Willis (HSW) estimates [START_REF] Te | The effect of inclusion shape on the elastic moduli of a twophase material[END_REF][START_REF] James | Theory of elastic properties of composite materials[END_REF][START_REF] Gatt | Elastic behavior of porous ceramics: application to nuclear fuel materials[END_REF], which are considered to be sufficient for the purposes of the present study.

Random process based on molecular dynamics

The basic idea of this process is due to Lubachevsky and Stillinger [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] and Lubachevsky et al. [START_REF] Lubachevsky | Disks vs. spheres : contrasting properties of random packings[END_REF] who proposed an algorithm based on molecular dynamics for random packings of disks and spheres in a square primitive cell. This algorithm was extended for a system of non-spherical particles in [START_REF] Donev | Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. i. algorithmic details[END_REF]. The algorithm treats all inclusions on an equal footing consisting in simultaneous generation of inclusions.

Recently, Ghossein and Lévesque [START_REF] Ghossein | Random generation of periodic hard ellipsoids based on molecular dynamics: A computationally-efficient algorithm[END_REF] were inspired by this approach to propose a fully automated numerical tool for a comprehensive validation of linear elastic homogenization models in the case of spherical or ellipsoidal particles reinforced composites. Vincent et al. [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: Ii. full-field simulations[END_REF] have used a similar random process based of the molecular dynamics for the study of the effective plastic flow surface of a fluid saturated bi-porous material. The main advantage of this algorithm is that high volume fractions (i.e. greater than 30% ) can be achieved. Polydisperse microstructures can be simply obtained by considering nonuniform growth rates for the inclusions [START_REF] Lavergne | Effects of grain size distribution and stress heterogeneity on yield stress of polycrystals: A numerical approach[END_REF].

General overview of the random sequential addition process

The algorithm proposed in the present study is based on a consecutive generation of inclusions in a periodic cell. It consists in placing randomly, irreversibly, and sequentially non-overlapping objects into a volume (or onto a surface) [START_REF] Talbot | Random sequential addition of hard spheres[END_REF][START_REF] Rintoul | Reconstruction of the structure of dispersions[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]. In this process, the acceptance of subsequent inclusion is constrained by some conditions reliant to the previously accepted inclusions and to the cell faces. Recent investigations have used the RSA algorithm in the context of linear and non-linear computational homogenization of composites or porous materials containing spherical particles or pores [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Böhm | Multi-inclusion unit cell models for metal matrix composites with randomly oriented discontinuous reinforcements[END_REF][START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF][START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF], monodisperse spheroidal inclusions [START_REF] Sherwood | Packing of spheroids in three-dimensional space by random sequential addition[END_REF] or even inclusions with octahedral, tetrahydral or cubic shapes [START_REF] Helmut | Effects of particle shape on the thermoelastoplastic behavior of particle reinforced composites[END_REF]. The drawback of this algorithm is the fact that in the case of monodisperse microstructure, it is difficult to achieve volume fractions greater than 30%. To overcome this limitation, Segurado and Llorca [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF] (see also [START_REF] Kari | Computational evaluation of effective material properties of composites reinforced by randomly distributed spherical particles[END_REF][START_REF] Sun | Finite element analysis of elastic property bounds of a composite with randomly distributed particles[END_REF]) developed a modified version of the RSA algorithm applied to spherical particles to reach higher volume fractions. In an alternative perspective and perhaps more closely related to analytical homogenization methods [START_REF] Willis | Bounds and self-consistent estimates for the overall modui of anisotropic composites[END_REF][START_REF] Willis | Variational and related methods for overall properties of composites[END_REF], Lopez-Pamies et al. [START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF] applied the RSA algorithm to polydisperse microstructures allowing for much higher volume fractions without any additional operations and without percolation thresholds (but also look at multiscale morphological models such as in [START_REF] Jean | Large-scale computations of effective elastic properties of rubber with carbon black fillers[END_REF] albeit for overlapping inclusions). Most of the studies in the literature pertain to two-phase composites. A recent extension of the RSA algorithms to N -phase composites that contain inclusions with different material constitutive response has been proposed in [START_REF] Papadioti | A methodology for the estimation of the effective yield function of isotropic composites[END_REF].

Scope of the study

In the present study, we propose an extension of the RSA algorithm to obtain random periodic (or not) distributions of mono-and polydisperse non-overlapping ellipsoidal inclusions of arbitrary shape and orientation. The centers of the ellipsoidal inclusions are randomly generated in a cuboidal cell following a uniform distribution probability. All ellipsoids that share the same shape and orientation shall belong to the same particle family. In the case of polydisperse distributions (i.e. inclusions of different sizes), one can reach with appropriate selection of the inter-distance parameters rather high volume fractions (> 60%) without need of additional operations contrary to previous studies [START_REF] Ghossein | Random generation of periodic hard ellipsoids based on molecular dynamics: A computationally-efficient algorithm[END_REF]. Then, we consider the numerical computation of the effective linear elastic properties of porous periodic microstructures by using a FFT-based numerical scheme and the FE method. We investigate particularly the deviation from isotropy and compare them with HSW estimates [START_REF] Willis | Bounds and self-consistent estimates for the overall modui of anisotropic composites[END_REF] which in some cases correspond to rigorous upper bounds. Voided microstructures are especially of interest in the context of geomaterials and rock physics. We show that the HSW estimates, which correspond to infinitely polydisperse matrix-inclusion systems, are extremely accurate for porous materials (but this is not the case for rigid particles [START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF]) even for monodisperse microstructures up to fairly high volume fractions (up to 35%). Beyond that volume fraction, the present algorithm needs to be modified (e.g. as in [START_REF] Ghossein | Random generation of periodic hard ellipsoids based on molecular dynamics: A computationally-efficient algorithm[END_REF]) to reach higher volume fraction of monodisperse inclusions but that is beyond the scope of this work. Instead, higher and more realistic high-volume-fraction microstructures can be easily attained by using polydisperse sizes of inclusions.

In Section 2, we discuss in detail the extension of the RSA algorithm in the context of random monoand polydisperse ellipsoids. Subsequently, in Section 3, boundary conditions are briefly presented and some details on the FFT and FEM methods, needed to carry out our numerical analysis of the effective elastic properties, are given. In Section 4, we show representative results for "almost" isotropic porous materials with random spherical, oblate and prolate voids and discuss their effect on the effective bulk and shear moduli, as well as their deviation from a purely isotropic response. The complete RSA algorithm is presented in appendix together with a convergence analysis of the FFT computations in terms of number of pores, voxel size and number of realizations.

Extension of the RSA algorithm for ellipsoidal heterogeneities

In this section, we extend previous approaches to include non-overlapping ellipsoids of arbitrary orientation and relative size. The extension consists in the generalization of the RSA algorithm proposed in [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF] for polydisperse inclusions. First, we define the microstructural parameters used in our process. Then, similar to most studies in the literature using the RSA algorithms, we introduce the geometric conditions imposed in order to allow for an adequate spacial discretization. Three types of geometric constraints are generally imposed:

1. For non-overlapping inclusions, the spacing between the inclusions must exceed a minimum value noted s 1 for proper discretization. Assuming periodicity of the inclusion distribution, this condition has to be checked 26 times for each pair of inclusions, i.e. between the newly added inclusion and any previously accepted inclusion as well as its periodic images near the opposite faces.

2. The inclusion surface must be sufficiently far from the cuboidal boundary faces to prevent, for example, the presence of distorted finite elements during meshing. This imposes that the distance from a point on the surface of the inclusion to any of the faces of the cuboidal cell must take a minimum value s 2 . Note, however, that this is less constraining in the context of FFT calculations mainly used in the present work.

3. In the case of periodic microstructures, any inclusion which intersects with any of the cuboidalcell boundary faces must be copied to the opposite face in order to impose periodicity of the microstructure.

Note that, by contrast with Pierard et al. [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF], the extension of the RSA algorithm proposed here takes into account an arbitrary orientation of inclusions. This algorithm can be used to generate random microstructures containing spherical, spheroidal or general ellipsoidal inclusions distributed in a monodisperse or polydisperse manner.

Microstructure description

Let us consider a cuboidal unit-cell (UC) with dimensions L 1 ,L 2 and L 3 and made of N p phases with volumes (r) and

V (r) (r = 1, • • • , N p ), such that Np r=1 V (r) = V . The UC comprises distinct families of ellipsoidal inclusions with volume fraction f (r) = V (r) /V , aspect ratios ω (r) 1 = c (r) /a
ω (r) 2 = c (r) /b (r) (r = 2, . . . , N p ) embedded in a matrix phase (r = 1). The vectors n (r) i (i = 1, • • • , 3)
form an orthonormal basis attached to the ellipsoid (Figure 1). A family of inclusions is a phase characterized by N (r) inclusions which have the same elastic properties, size, shape and orientation.

For simplicity in the following, attention is restricted to inclusions which have the same constitutive properties. The extension to different ones is straightforward and has been described for spherical inclusions in [START_REF] Papadioti | A methodology for the estimation of the effective yield function of isotropic composites[END_REF]. In this study, a monodisperse periodic microstructure is represented by a UC with randomly orientated inclusions having the same size and shape while a polydisperse microsctructure refers to inclusions of different (relative) size and/or shape. Following [START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF] for the generation of polydisperse microstructure, we define a reference size of inclusion which is used to calibrate every other inclusion phase/size. The reference size, for consistency with the monodisperse microstructures, is defined via a reference number of monodisperse inclusions, N ref , so that N ref ≤ N with N = Np r=2 N (r) being the total number of inclusions with volume fraction f = Np r=2 f (r) and aspect ratios ω ref ) allows to obtain uniquely the reference size of the inclusions in the unit-cell. In the following, we identify the different cases resulting from the above descriptions.

Polydisperse only in size

We consider inclusions of different size but same shape implying that the shape of all inclusions coincides with the reference one:

ω ref 1 = ω (r) 1 and ω ref 2 = ω (r) 2 (r = 2, • • • , N p )
. For instance, one can take the first family of inclusions to be the reference one and set one of its axes, say for instance c ref equal to 1. The lengths of the semi-axes of the reference inclusions are then given by

a ref = c ref ω ref 1 , b ref = c ref ω ref 2 , c ref = 3 V f ω ref 1 ω ref 2 4πN ref 1/3 . ( 1
)
The semi-axes lengths of the other inclusions are determined by using c ref and the size coefficient χ (r) , such that

c (r) = χ (r) c ref a (r) = c (r) ω (r) 1 , b (r) = c (r) ω (r) 2 , χ (r) ≤ 1. (2)
In practice, the microstructure is generated by starting with the phase (r) which has the largest size 130 χ (r) until its concentration f r is reached. The process is continued sequentially with the second largest phase and so on until Np r=2 f (r) f . Note that contrary to the monodisperse case, the targeted volume fraction is only approximately attained. The addition of inclusions of decreasing sizes allows to get a better accuracy on the volume fraction and to reach higher volume fractions.

Polydisperse only in shape 135

In this case, the inclusions are defined by different shapes (i.e. different ω (r) [START_REF] Milton | The Theory of Composites[END_REF] and ω (r) 2 for different r) but have the same volume. The reference inclusion is then for simplicity taken to be a sphere of radius R ref , and the calibration is done with respect to the volume of this sphere by setting

R ref = 3V f ref 4πN ref 1/3 . ( 3
)
This amounts to define the characteristic lengths of the semi-axes of each family of inclusion as

a (r) = c (r) ω (r) 1 , b (r) = c (r) ω (r) 2 , c (r) = ω (r) 1 ω (r) 2 1/3 R ref . ( 4
)
In this case, there is no preferential order in the random generation of inclusion/phases since they all have the same relative size. One can then begin with the phase r which has the smallest volume fraction f (r) and progress until Np r=2 f (r) f . It is worth noticing that instead of using a reference number of inclusions N ref to define a reference size, one can use directly a reference radius R ref .

Polydisperse in shape and size

This case is a simple combination of the two previous cases. In this general case, the relative size of each family of inclusions is determined using V (r) = 4π(χ (r) R ref ) 3 /3, where R ref is given by (3).

Again, the sequential addition process needs not to have a preferential order neither in terms of size nor in shape. Nonetheless, as mentioned previously, starting with the largest size of inclusion proves useful for attaining a better accuracy on the volume fraction as well as for reaching higher compacities.

Description of the extended RSA algorithm

In this section, we describe briefly the main steps of the extended RSA algorithm (see A for more details). The inputs of the algorithm are: the dimensions of the cuboidal cell L 1 , L 2 and L 3 , the volume fraction f , the number of phases N p , the number N ref of reference inclusions, the microstructural parameters ω

(r) 1 , ω (r)
2 , χ (r) , f r and two offsets distance ξ 1 and ξ 2 used in the calculation of the minimum distance parameters s 1 and s 2 , respectively. In this work, these parameters are fixed at ξ 1 = 0.02 and ξ 2 = 0.05, however they are user-defined and can be changed at will. In general, the highest the volume fraction the smallest should be ξ 1 and ξ 2 . When the desired volume fraction of the inclusions is moderate to low (< 20%), one can increase those values to obtain more uniform distributions. The outputs of the algorithm are: the position vector of the center of the ellipsoid v (r) i , the semi-axes lengths (a

(r) i , b (r) i , c (r) i
and the Euler angles (φ

(r) i , θ (r) i , ψ (r)
i ) of the ellipsoidal inclusion i belonging to the phase r. To determine periodic images of an inclusion, 26 vectors m = (m 1 , m 2 , m 3 ) are defined, where m 1 , m 2 and m 3 take the values (0, -L 1 , L 1 ), (0, -L 2 , L 2 ) and (0, -L 3 , L 3 ), respectively. Specifically, for each phase r = 2, • • • , N p , the RSA algorithm can be decomposed into four steps:

Step 1: Compute the semi-axis lengths of inclusion (i)

c (r) i = χ (r) ω (r) 1 ω (r) 2 1/3 R ref , b (r) i = c (r) i /ω (r) 2 and a (r) i = c (r) i /ω (r) 1
with R ref given by (1). For non-spherical inclusions, specify the orientation with Euler angles (φ

(r) i , θ (r) i , ψ (r) i ).
Step 2: In the sequential addition, generate a random center position vector v (r) i for inclusion i in phase r. Compute the minimum distance ∆ 1 between a new inclusion i and any previously accepted inclusion j = 1, . . . , i -1 including its 26 periodic images, and compare this distance to the minimum value s 1

s 1 = max a (r) i , b (r) i , c (r) i + max a (r) j , b (r) j , c (r) j * ξ 1
In this step, the algorithm for finding the minimum distance between two ellipsoids (see Algorithm Step 3: Use the algorithm for finding the minimum distance between ellipsoid and plane (see Algorithm 3 in A) to determine ∆ 2 and ∆ 2 . Compare the distance ∆ 2 with the minimum value s 2

2 in A) is called. If ∆ 1 < s 1 ,
s 2 = (max (a (r) i , b (r) i , c (r) i )) * ξ 2 If ∆ 2 < s 2 ,

then generate a new center position of inclusion i, otherwise accept the inclusion i.

Step 4: Ensure periodicity of the cuboidal cell by considering periodic images of the inclusion (see Section 2.5 below).

In the following, we describe in details the last three steps of the RSA procedure.

RSA Step 2: Minimum distance between two ellipsoids

In this section, we discuss the evaluation of the distance between two ellipsoidal inclusions. This step is extremely important since it allows to decide if a newly added inclusion is accepted or rejected. Specifically, the distance between two non-spherical inclusions is not simply determined by using the center-to-center distance, as is the case with spheres. There are many studies in literature devoted to the determination of the minimum distance between two implicit algebraic surfaces. For example, in the recent work of Chen et al. [START_REF] Chen | Computing minimum distance between two implicit algebraic surfaces[END_REF], a useful analytical method has been proposed to compute this distance. The drawback is that this method leads to a very complicated system of three non-linear equations of degrees 2, 3 and 6 to be solved, at least in the case of two ellipsoids. In the work of Pierard et al. [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF] on elasto-plastic composite materials reinforced with aligned elastic ellipsoidal particles, an algorithm proposed by Lin and Han [START_REF] Lin | On the distance between two ellipsoids[END_REF] is used to determine iteratively the minimum distance between two ellipsoids. This algorithm based on the local approximation of the ellipsoid by a spherical ball has excellent convergence properties. This approach is simple and easy to implement in the present study. Hence, following [START_REF] Lin | On the distance between two ellipsoids[END_REF], we implement and embed in the RSA algorithm a slightly modified version of the algorithm. The modification consists in a much simpler and more general way to check the non-overlapping condition of inclusions during the sequential addition. In Lin's algorithm, the non-overlapping is checked by solving two one-dimensional quadratic inequalities at each iteration.

With our modification, we are able to check this condition once by using a lemma on the intersection of ellipsoids proposed in the book of Kurzhanski and Vályi [START_REF] Kurzhanski | Ellipsoidal calculus for estimation and control[END_REF]. For clarity, in this section, the superscript (r) denoting the different ellipsoidal families is omitted.

Statement of the problem

Given two non-overlapping arbitrarily oriented ellipsoids E 1 and E 2 defined respectively by

E n = E(v n ; Z n ), n = 1, 2, such that E(v n ; Z n ) := {x : Q n (x) ≤ 0} . ( 5
)
Q n denotes a quadratic function which reads

Q n (x) = (x -v n ) T Z n (x -v n ) -1 (6)
with v n rhe vector position of the center of the ellipsoid. In turn, the matrix Z n describess the shape and the orientation of the ellipsoid E n and is a positive definite square 3 × 3 matrix. Its eigenvalues are the reciprocals of the squares of the semi-axes lengths of the ellipsoid (i.e., a, b, c and equivalently the aspect ratios ω 1 and ω 2 ) and its eigenvectors n i (i = 1, 2, 3) define the principal axes of the ellipsoid E n (Figure 1). For an arbitrarily oriented ellipsoid E(v; Z) with center position v, the matrix Z is defined as

Z = 1 c 2 ω 2 1 n 1 ⊗ n 1 + ω 2 2 n 2 ⊗ n 2 + n 3 ⊗ n 3 . ( 7
)
The orientation of the ellipsoid is specified by the rotation matrix R(φ, θ, ψ) such that n i = Re i (i = 1, 2, 3). (e 1 , e 2 , e 3 ) defines an orthonormal basis of the reference frame (Figure 1). When R = I, E is thus an axis-aligned ellipsoid (i.e. its principal axes coincide with the reference frame axes). The minimum distance between two ellipsoids reads

d(E 1 , E 2 ) = min x1∈Ω1, x2∈Ω2 x 1 -x 2 (8)
where Ω 1 and Ω 2 represent the boundaries of the ellipsoids E 1 and E 2 .

Algorithm on the local approximation of the ellipsoid by ball

Before we proceed further, it is recalled that one of the main conditions that need to be checked is that the distance between the ellipsoids shall not exceed a user-defined value s 1 . Therefore, the procedure of evaluating the distance between two different ellipsoids can be split in two main parts.

The first part consists in checking whether the ellipsoids overlap/intersect. Provided that the ellipsoids

195
do not overlap, one can then proceed to the second part and evaluate the distance between the two ellipsoids.

First part. In the present work, we propose a simple method to check if the ellipsoids intersect at one or more point by using the lemma of Kurzhanski and Vályi [START_REF] Kurzhanski | Ellipsoidal calculus for estimation and control[END_REF] which reads:1

Two ellipsoids E 1 and E 2 have a nonempty intersection if, and only if, there exists an ellipsoid

E α = E(v α ; Z α ), for any α ∈ [0, 1], such that det(Z α ) > 0 with Z α = 1 1 -h α 2 n=1 α n Z n , α 1 = α, α 2 = 1 -α, ( 9
)
where

h α = 2 n=1 α n v n •Z n •v n - 2 n=1 α n Z n •v n •v α and v α = 2 n=1 α n Z n -1 2 n=1 α n Z n •v n .
The use of this lemma has the advantage of checking only once the overlapping condition of any pair 200 of ellipsoids.

Second part. If the ellipsoids do not intersect, we proceed to solve the minimization problem [START_REF] Kröner | Statistical Continuum Mechanics[END_REF].

The basic idea for solving numerically this problem is due to Lin and Han [START_REF] Lin | On the distance between two ellipsoids[END_REF], who proposed an iterative resolution in which each ellipsoid E n is approximated locally in its interior by a spherical ball which marches along the internal surface of the ellipsoid, i.e. the spherical ball is tangent to the ellipsoid at a single point x n ∈ Ω n (Figure 2). The spherical ball

B n = B(c n ; r n ) (n = 1, 2) is defined by B(c n ; r n ) := {y : y -c n ≤ r n } , ( 10
)
with

c n = x n - γ n 2 N n and r n = γ n 2 N n .
Here, N n = ∇Q n is the normal to the ellipsoid E n at the point x n ∈ Ω n . The parameter γ n is a measure (matrix norm) of the matrix Z n related to its spectral radius ρ(Z n ), such that

0 < γ n ≤ 1 ρ(Z n ) and γ n = 1 Z n F ( 11
)
where • F indicates the Frobenius (Euclidean) norm.

The algorithmic process involving both parts is schematically explained in Figure 2 and can be described as follows:

Step 1: Check the overlapping at the beginning of iterations using [START_REF] Sen | Effective conductivity of anisotropic two-phase composite media[END_REF]. If det(Z α ) > 0, then set the Initialization: Choose the centers v 1 and v 2 of E 1 and E 2 as starting interior points c 1 and c 2 , respectively.

Step 2: Let us denote by c k 1 and c k 2 two arbitrary interior points at the k th iteration of the ellipsoids E 1 and E 2 , respectively. We generate two points x k+1

= c k 1 + t 1 (c k 2 -c k 1 ) and x k+1 2 = c k 1 + t 2 (c k 2 -c k 1 )
, where t 1 and t 2 are solutions of the two one-dimensional quadratic equations:

t n = t ∈ [0, 1] : A n t 2 + B n t + C n = 0 , n = 1 or 2, ( 12
)
with

A n = (c k 2 -c k 1 ) T Z n (c k 2 -c k 1 ), B n = (c k 1 -v n ) T Z n (c k 2 -c 1 ) and C n = (c k 1 -v n ) T Z n (c k 1 -v n ).
These two points are the closest points and hence the best candidate of the minimum distance if Step 3:

If the angles ∠(N k+1 1 , x k+1 2 -x k+1 1 ) and ∠(N k+1 2 , x k+1 1 -x k+1 2 ) are null, set d(E 1 , E 2 ) = x k+1 1 -x k+1
c k+1 n = x k+1 n - γ n 2 N k+1 n , n = 1, 2 (13)
and return to Step 2.

RSA Step 3: Minimum distance between an ellipsoid and a plane

In this section, we determine the minimum distance between an ellipsoid E and a plane P. It is worth noting that if the ellipsoid intersects a plane, the terminology minimum distance from a point lying on an ellipsoid to a plane is more appropriate. Then, this distance is obtained by solving a geometrical problem for finding the closest point on the boundary Ω of the ellipsoid E where the normal N = ∇Q(x) is collinear to the normal of the plane. Knowing this closest point, one can verify the second geometric condition and then impose the periodicity condition mentioned previously.

Any surface of the cuboidal cell is defined by a plane equation given as

P(n, η) := {y : y • n + η = 0} ( 14
)
where n is the normal unit vector of the plane P, and d is a constant obtained as a scalar product of n and a known point on P. For a cuboidal cell with dimensions L 1 , L 2 and L 3 , with a corner as origin, let us adopt the following convention: n = (-1, 0, 0) and η = 0 → face 1, n = (1, 0, 0) and η = -L 1 → face 2, n = (0, -1, 0) and η = 0 → face 3, n = (0, 1, 0) and η = -L 2 → face 4, n = (0, 0, -1) and η = 0 → face 5 and n = (0, 0, 1) and η = -L 3 → face 6. The problem to solve then reads

N × n = 0 and Q(x) = 0. ( 15
)
The solutions of (15) are two points x 1 and x 2 representing the closest and the farthest points, or vice versa depending on the face of the cuboidal cell. The related minimum or the maximum distance, denoted as d(E, P), is obtained by using the formula giving the distance of a point x to a plane

d(E, P) = |x • n + η| n (16)
Obviously, the farthest point from a plane with η = 0 is the closest one to the opposite plane with η = -L 1 . In practice, to distinguish the nearest and the farthest points, one can determine the two distances corresponding to the two solutions at fixed n and η and then take the minimum. We denote by

x min the solution (15) which minimizes the distance d(E, P). The corresponding algorithm is detailed in A. ∆ 2 is the minimum distance and ∆ 2 denotes the quantity n•x min + η.

RSA Step 4: Enforce periodicity of the cuboidal cell

To impose the periodicity of the microstructure in the cuboidal cell, any inclusion i that intersects any of the cuboidal cell faces j defined by the normal n and η is relocated at the opposite face by adding or subtracting the vector m to the center of the inclusion. To do this, the sign of ∆ 2 = n • x min + η must be checked. If the inclusion intersects the face j, ∀ j ∈ {1, 3, 5}, ∆ 2 < 0, and ∀ j ∈ {2, 4, 6}, to the three faces that share the same normals. Finally if the ellipsoid intersects only one face, it is duplicated to the opposite face of the cuboid. (L 1 = L 2 = L 3 = L = 1) with periodic boundary conditions. We analyse monodisperse and polydisperse microstructures consisting of random and uniform void distributions with three types of shapes; (i) spherical pores, (ii) prolate spheroidal pores (ω 1 = ω 2 = ω > 1) with aspect ratios ω = 2 and ω = 5 and (iii) oblate spheroidal pores with aspect ratios ω = 0.1 and ω = 0.5. The effective elastic properties of those microstructures are obtained through numerical periodic homogenization. In all subsequent calculations, the elastic properties of the matrix phase are set to: Young's modulus E = 1 GPa and Poisson's ratio ν = 0.3.

Illustrative examples

Numerical homogenization

The numerical computation of the linear effective elastic properties is carried out using the fast Fourier Transform (FFT) numerical scheme and the finite element (FE) method. The FEM calculations have been performed with the commercial software ABAQUS in the context of linear elasticity and small strains. The application of the periodicity conditions is discussed in the next section. The FFT-based numerical approach has been proposed by Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] to solve periodic boundary-valued heterogeneous unit-cell problems. The principle of this method is to solve iteratively the implicit integral equation for the strain field ε(x)

ε(x) = E + V Γ (0) (x -x ) : τ (x )dx , τ (x) = (C(x) -C (0) ) : ε(x) (17)
with E the macroscopic strain and Γ (0) denoting the strain Green operator corresponding to a reference homogeneous medium with elasticity C (0) . Following [START_REF] Kröner | Bounds for effective elastic moduli of disordered materials[END_REF], [START_REF] Lubachevsky | Disks vs. spheres : contrasting properties of random packings[END_REF] is usually called the Lippman-Schwinger equation by analogy with the quantum mechanics scattering theory. Developments of the FFT-based method have been proposed by several authors, in particular to improve its convergence in the case of a high mechanical contrast on the local properties [see, among others, [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement*[END_REF][START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF][START_REF] Brisard | Fft-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. In the present study, we have used the augmented Lagrangian scheme originally proposed by Michel et al. [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF] and later reinterpreted by Moulinec and Silva [START_REF] Moulinec | Comparison of three accelerated fft-based schemes for computing the mechanical response of composite materials[END_REF] as a special case of the polarization-based scheme of Monchiet and Bonnet [START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF].

At this point it is important to note that the FE approach may require tedious meshing that more than often does not converge, especially for penny-shaped (resp. elongated) inclusions presenting very low (resp. high) aspect ratios. This is the case in the present work for oblate spheroids with aspect ratio lower than 0.3 and for prolate spheroids with aspect ratio higher than 3. By contrast, in the FFT method, a regular grid is employed to discretize the unit-cell (i.e. voxel-based meshing). To describe accurately the highly elongated ellipsoids, a large number of voxels is required but finally yields to converged results. A convergence analysis with respect to the number of voxels, the number of pores and the number of realizations is presented in B. To assess the relative accuracy of the results obtained with the FFT-based method and the FE method, we carry out in Section 3.2 calculations for randomly oriented, monodisperse pore microstructures with aspect ratios equal to 0.5, 1 and 2.

Periodic Boundary conditions and effective stiffness tensor

For both numerical techniques, the unit-cell is subjected to periodic boundary conditions [START_REF] Suquet | Homogenization techniques for composite media[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF].

The displacement field u can be split into an affine part E•x and a periodic correction term u * (x) such that

u(x) = E•x + u * (x), ∀ x ∈ V, ( 18
)
where E is the overall strain which is equal to the average strain in the unit-cell (E = ε )2 and u * is a periodic field characterizing the fluctuation of the displacement due to the presence of the heterogeneities. The local strain field ε(x) deriving from the displacement u thus admits the following

decomposition ε(x) = E + ε * (x) with ε * (x) = ∇ s u * (x) = 0.
The components of u * take identical values at points on opposite faces of the unit-cell, in such a way that, considering a Cartesian frame of reference with origin placed at a corner of the cubic unit cell and axes e k (k = 1, 2, 3) are aligned with the principal axes of the cubic unit cell, we have [START_REF] Mbiakop | An analytical model for porous single crystals with ellipsoidal voids[END_REF]

u k (L, x 2 , x 3 ) -u k (0, x 2 , x 3 ) = E k1 L, u k (x 1 , L, x 3 ) -u k (x 1 , 0, x 3 ) = E k2 L, u k (x 1 , x 2 , L) -u k (x 1 , x 2 , 0) = E k3 L (19)
with k = (1, 2, 3). These relations constitute simple linear constraints that can be entered in any numerical code with a number of available techniques, such as elimination, Lagrange multipliers and penalty methods. In the context of our FE calculations in the Abaqus software, the linear constraint equations [START_REF] Ghossein | Random generation of periodic hard ellipsoids based on molecular dynamics: A computationally-efficient algorithm[END_REF] are implemented by use of the * Equation command, which uses the elimination technique.

As a practical remark for the finite element calculations, keeping the origin of the unit-cell fixed, loadings were done at its three corners having coordinates (L, 0, 0), (0, L, 0) and (0, 0, L) since the periodic boundary conditions ( 19) can be rewritten in terms of the displacements of these three corner nodes (see details in Appendix B of [START_REF] Danas | Effective response of classical, auxetic and chiral magnetoelastic materials by use of a new variational principle[END_REF]). The equations relating the displacements of the nodes on opposite faces of the unit-cell were coupled to those of the corner nodes. Note that the above procedure can be applied to any cuboidal shape. In the special case of voids, one should make sure that no void (nodeless point) is lying in any of the above four corners.

Once the local stress and strain fields are calculated, the effective stiffness tensor C is obtained with the relation between average stress and strain tensors

σ = C : ε (20)
The computation of the twenty-one independent coefficients of C is classically done by considering six independent macroscopic strain loadings.

Deviation from isotropy

Theoretically, a completely random distribution and orientation of ellipsoidal inclusions leads to an isotropic effective response. In practice, due to the finite number of inclusions and/or pores in the unit-cell, the resulting numerical effective stiffness tensor C is not exactly isotropic. To evaluate the deviation from isotropy of the tensor C, use is made of its isotropic projection C iso on the fourth-order deviatoric and hydrostatic isotropic tensors

C iso = 3 κ J + 2 µ K with κ = 1 3 C :: J = 1 9 C iijj and µ = 1 10 C :: K = 1 10 C ijij -9 κ . ( 21
)
Here, κ and µ are the isotropic bulk and shear moduli. The isotropic projectors are defined by

J = 1 3 i ⊗ i and K = I -J (22)
with i and I the identity tensors, respectively for symmetric second and fourth-order tensors. Note that J and K are orthogonal idempotent tensors, i.e J : J = J, K : K = K, and J : K = K : J = O.

The deviation from isotropy can be quantified in two ways: either geometrically (material symmetry)

or mechanically (elasticity symmetry). On the one hand, the geometrical approach3 makes use of morphological descriptors such as statistical correlation functions to characterize statistically the microstructure (e.g. two-point probability function is often used to quantitatively ascertain at least lower-order information on the geometrical arrangement of phases). On the other hand, the mechanical approach consists in evaluating the deviation from elastic isotropy by choosing an appropriate measure of the elasticity tensor. One of the first measure of this kind is the well-known Zener ratio [START_REF] Zener | Elasticity and anelasticity of Metals[END_REF], for cubic symmetry materials, which is simply the ratio of the shear moduli. Later on, several authors have adressed the more general question of the distance between two elasticity tensors of given symmetry [see, for instance, [START_REF] Gazis | The elastic tensor of given symmetry nearest to an anisotropic elastic tensor[END_REF][START_REF] François | Determination of the symmetries of an experimentally determined stiffness tensor: Application to acoustic measurements[END_REF][START_REF] Moakher | The closest elastic tensor of arbitrary symmetry to an elasticity tensor of lower symmetry[END_REF][START_REF] Bucataru | Invariant properties for finding distance in space of elasticity tensors[END_REF]. At this point, it is worth recalling that the symmetry of a physical property (elasticity, for instance) is necessarily equal or higher than the material symmetry. Consequently, the geometrical isotropy implies the mechanical isotropy but not the other way around. In addition, in the RVEs studied here one can never obtain an exact geometrical isotropy due to the finite number of the embedded inclusions and the periodicity of the unit cell. In this case then, it is extremely difficult to estimate quantitatively the deviation from mechanical isotropy (which is the measure of interest in this study) by using the deviation from geometrical isotropy.

The present study is, thus, focused on the analysis of the isotropy of the overall elastic response for different random distributions of ellipsoidal voids in the matrix. The deviation from isotropy δ iso of the effective elasticity tensor C is evaluated with the normalized common Euclidean distance

δ iso = C -C iso F C F ( 23
)
where A F = tr A : A T is the Frobenius norm of the tensor A. The value of δ iso = 0 thus corresponds to isotropy. This approach has been previously used, for example, in [START_REF] Ghossein | A comprehensive validation of analytical homogenization models: The case of ellipsoidal particles reinforced composites[END_REF] for the case of linear elastic response of ellipsoidal particles reinforced composites, as well as for nonlinear mechanical responses [START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF][START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: Ii. full-field simulations[END_REF].

Assessment of FFT and FE numerical computations

FE calculations of the effective elastic properties were carried out using the FE package [56], only on monodisperse microstructures with pores aspect ratios 0.5, 1 and 2 and volume fractions 5% and 15%.

The mesh generator NETGEN [START_REF] Schöberl | Netgen an advancing front 2d/3d-mesh generator based on abstract rules[END_REF] was used to create meshes of the unit-cell with ten-node quadratic tetrahedral elements (C3D10 in ABAQUS notation). For certain microstructures containing spheroidal pores, we faced severe convergence issues with the meshing procedure when using quadratic ten-node elements. To overcome these difficulties, we have resorted to linear four-node tetrahedral elements (C3D4 in ABAQUS notation) and substantial increase of the number of elements to compensate for the inaccuracy of the linear shape functions. For microstructures with spherical voids, a number of pores N = 30 is found to be sufficient for convergence. However, for non-spherical voids and large porosities, convergence could require a much larger number of pores (more than N = 500). Unfortunately, in that case the FEM calculations become substantially heavy. To keep the number and time of FEM calculations tractable, we have chosen to use N = 100 for microstructures with spheroidal pores to get only an idea of the FEM outcome. This issue has motivated the use of the FFT numerical scheme in the rest of this study.

The FFT calculations have been performed on a regular grid with 128 × 128 × 128 voxels. The choice of this spatial resolution is justified by the convergence analysis in terms of number of voxels presented in B (Figure 1). In particular, for monodisperse prolate and oblate spheroids with aspect ratios ω 1 = ω 2 = 2 and 0.5, respectively, a number of pores N = 500 was found to be sufficient for achieving representativity (Figure 2). In the case of spherical voids, N = 30 pores are sufficient to obtain accurate results as for the FE computations. Note that the material properties have been assigned to each voxel as follows: if the center of a voxel does not belong to any ellipsoid, then the entire voxel has the mechanical properties of the matrix, otherwise the voxel has null mechanical properties (void). It can be mentioned that improved methods allowing for the use of "filtered" elastic properties (smoother transitions) between neighboring phases have been proposed [START_REF] Brisard | Fft-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Gélébart | Filtering material properties to improve fft-based methods for numerical homogenization[END_REF].

In Table 1, the effective bulk and shear moduli ( κ, µ) computed with the FFT method are compared to those calculated with the FEM. The comparison is done by calculating the relative differences δ κ and δ µ . Note that in all calculations, the effective moduli obtained with the FFT scheme are smaller than the one obtained with the FE method. It can be observed in Table 1 that for all microstructures the differences are relatively small, keeping in mind that the FE results are obtained for a smaller number of pores than the FFT ones for ellipsoidal voids. Consistently, in the case of spherical voids, the FE and FFT results are in very good agreement. For the reasons evoked above, all results reported in the sequel have been obtained with the FFT numerical method. In closing this section, it is relevant to mention that the present RSA algorithm can deal with extremely elongated inclusion shapes, such as aspect ratios of 0.001 or 1000. Nonetheless, their discretization becomes prohibitive in the FEM context, while extremely large numbers of voxels would be required to achieve convergence with the FFT formulation. Furthermore, one should also be careful when analyzing very elongated ellipsoids since they tend to cracks and thus can lead to numerical singularities. In those cases, perhaps the analytical approach would be more appropriate, especially in the linear elasticity context.

f (%) Monodisperse microstructures Sphere Oblate (ω = 0.5) Prolate (ω = 2) (δ κ , δ µ ) (%) (δ κ , δ µ ) (%) (δ κ , δ µ ) (%)

Results

The FFT results have been obtained with N ref = 50 for the spherical voids and N ref = 500 for the spheroidal voids. To make the present section concise but representative, we have chosen to consider spheroidal prolate and oblate void (i.e. ω 1 = ω 2 = ω). In the sequel, ω = 1 refers to spherical, ω < 1 for oblate and ω > 1 for prolate voids.

Next, it is recalled that N ref = N = 500 for the case of monodisperse microstructures. For the case of polydisperse microstructures, we have used four different phases or equivalently four different relative sizes of voids (but same shapes) with size coefficient χ (r) = 1, 0.7, 0.4, 0.1 (i.e.

the volume of the larger inclusion is ten times the one of the smaller one) and volume fraction distribution f (r) = {0.5, 0.3, 0.15, 0.05} × f such that N r=2 pf (r) = f (with N p = 4 being the number of families/phases). A value of N ref = 500 leads on average to 800-1500 inclusions of different sizes.

The final number depends on the random process, the choice of the number of families/phases N p and the size coefficient. In this connection, it is perhaps important to mention that at high volume fractions, where lies the main difficulty in obtaining elongated inclusion shapes, the resulting polydisperse microstructures with different final number of voids provide almost identical responses, which for the case of voids are in close agreement with the corresponding monodisperse ones. This is a strong indication of convergence given the N ref used in our calculations. In addition to this observation, the deviation between various realizations both in mono-and polydisperse microstructures is extremely small and for that reason no error bars are shown in the results that follow.

The numerical FFT results are also compared with HSW analytical estimates which are obtained by employing the approach proposed by Gatt et al. [START_REF] Gatt | Elastic behavior of porous ceramics: application to nuclear fuel materials[END_REF], see also [START_REF] Walpole | On the overall elastic moduli of composite materials[END_REF], for isotropic composites comprising randomly oriented polydisperse ellipsoidal voids. In Figure 4, the deviation from isotropy δ iso (c.f. equation ( 23)) is shown for monodisperse and polydisperse microstructures. Hence, the spherical voids lead to the most isotropic effective response with δ iso < 0.5% for both monodisperse and polydisperse distributions. As one increases or decreases the aspect ratios from unity to get either prolate or oblate shapes, the deviation from isotropy increases. Even for a high number of inclusions, the transversely isotropic shape of the void affects the evaluation of the elastic properties. The deviation is markedly strong for oblate voids with aspect ratios ω = 0.1, where δ iso ∼ 24% for both monodisperse and polydisperse distributions at volume fractions f > 15%. This implies that even if theoretically one can consider very large number of orientations, in practice, with our present RSA method, a significant effect of the void shapes remains and lead to an anisotropic response.

In the rest of the study, we show only the isotropic parts corresponding to spheroidal voids and compare them with exactly isotropic (i.e. infinite number of orientations) analytical results. Nonetheless, their interpretation should be done with caution since the numerical estimates exhibit in some cases (such as ω = 0.1) strong deviation from isotropy. Moreover, it is noted that using a larger number of realizations does not reduce significantly the deviation from isotropy in these extreme cases since already a very large number of voids (∼ 1500) is considered in a single unit-cell.

Numerical versus analytical estimates for spherical, prolate and oblate pores

Figure 5 shows numerical (FFT) and analytical (HSW) results for monodisperse and polydisperse spherical voids. We observe that up to f = 30%, the numerical results for monodisperse and polydisperse microstructures are almost identical and only slightly lower than the HSW estimates, which in that case are also rigorous upper bounds. Of course, it is progressively more and more difficult to generate monodisperse distributions when the targeted volume fraction f gets closer to the theoretical saturation value f sat = 38.4% [START_REF] Zhang | Precise algorithm to generate random sequential addition of hard hyperspheres at saturation[END_REF]. By contrast, the polydisperse microstructures can easily reach higher volume fraction. In this work, results are shown up to f = 50% but it is indeed possible to reach even higher volume fraction [START_REF] Danas | Numerically-aided 3d printed random isotropic porous materials approaching the hashin-shtrikman bounds[END_REF]. In Figure 6, FFT and HSW estimates for prolate monodisperse and polydisperse voids with aspect ratios ω = 2 and 5 are reported. First, it should be noted that by referring to Figure 4, the maximum deviation from isotropy occurs for ω = 5 and is δ iso = 6%. Therefore, the curves shown correspond to only the isotropic part of the elastic properties as defined in [START_REF] Lavergne | Effects of grain size distribution and stress heterogeneity on yield stress of polycrystals: A numerical approach[END_REF]. Again a minor difference is observed between the mono-and polydisperse microstructures with the later being only slightly more compliant.

The numerical estimates are still lower than the HSW ones with their difference being larger for the shear modulus µ. It is interesting to observe that the HSW analytical estimates are very close to the numerical estimates but it is recalled that for ω = 0.1 the deviation from isotropy is very large δ iso = 24%. Moreover, one can observe that the effect of the aspect ratio is more pronounced on the bulk modulus κ than in the shear modulus µ. We also note that in such extreme cases of elongated ellipsoids, it is progressively more difficult to reach higher volume fractions even with polydisperse sizes. For instance, in the case of ω = 0.1, we were able to reach f = 20% for the monodisperse and f = 25% for the polydisperse.

Effect of aspect ratio in polydisperse microstructures

Figure 8 shows a comprehensive plot of the effect of the aspect ratio for polydisperse microstructures upon the effective bulk and shear moduli. We observe that spheroidal (prolate and oblate) voids tend to make the material more compliant as already seen in a number of studies [see, for instance, 15]. The oblate voids tend to lead to stronger drops of both κ and µ than prolate ones but one again should recall the strong deviation from isotropy for the oblate voids. We close by noting that the spherical voids lead to the less compliant response and hence are better candidates for materials that can approach the Hashin-Shtrikman bounds.

ω = 0.1 δ max iso = 24% ω = 0.5 δ max iso = 6% ω = 1 δ max iso = 0.1%% ω = 2 δ max iso = 4% ω = 5 δ max iso = 6% (

Concluding Remarks

In the present study, we extend the classical Random Sequential Adsorption (RSA) algorithm to the context of randomly oriented ellipsoidal inclusions of arbitrary size, shape and orientation (i.e.

425

polydisperse). The main algorithm is based upon a proper simplification of Lin and Han [START_REF] Lin | On the distance between two ellipsoids[END_REF] algorithm, which consists in employing a more general way to check the non-overlapping condition of inclusions during the sequential addition process. In that algorithm, the non-overlapping condition is checked by solving two one-dimensional quadratic inequalities at each iteration. Instead, with our modification, we are able to check this condition once by using a lemma on the intersection of ellipsoids proposed in the 430 book of Kurzhanski and Vályi [START_REF] Kurzhanski | Ellipsoidal calculus for estimation and control[END_REF]. The proposed algorithm is described in detail, is robust and can be easily implemented.

Secondly, the present study analyzes numerically the effective elastic properties of a large number of representative volume elements (RVEs) including monodisperse as well as polydisperse voids randomly distributed in the unit cell. Those numerical estimates are probed by the classical Hashin-Shtrikman-Willis (HSW) bounds and estimates. The deviation from isotropy of the numerically-obtained effective elastic tensors is critically assessed by using a simple measure based on a normalized Euclidean distance between two fourth order tensors and the Frobenius norm. In simple words, we project the numerically obtained effective elastic tensor (which may exhibit certain degree of anisotropy) to the isotropic fourth order tensor space to obtain an "isotropized" one. The Euclidean difference of the original tensor and the isotropized one, gives a simple to understand scalar measure of the deviation of isotropy.

In this regard, we find that as the aspect ratios of the voids increase (or decrease) substantially, the proposed RSA algorithm leads to substantially anisotropic porous materials which can deviate from isotropy by as much as 24% for oblate voids with aspect ratio ω 1 = ω 2 = 0.1. This deviation is present even for very large number of voids N > 500 both in the case of monodisperse and polydisperse microstructures. The polydisperse microstructures exhibit in principle lower anisotropy albeit an important one for ω 1 = ω 2 = 0.1 (reaching 20%). Interestingly, the isotropized shear and bulk moduli corresponding to the very elongated void aspect ratios remain in good agreement with the corresponding isotropic HSW estimates. As expected by previous works, the spherical voids exhibit the stiffest response out of all cases considered.

In closing, we remark that the present RSA microstructures can be used as a test bed for nonlinear analytical homogenization models in the context of finite strains and nonlinear mechanical and multiphysics responses (see for instance studies in nonlinear elasticity [START_REF] Lopez-Pamies | The nonlinearelasticresponseofsuspensionsofrigidinclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF], in elasto-plasticity [START_REF] Mbiakop | An analytical model for porous single crystals with ellipsoidal voids[END_REF][START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF], as well as in nonlinear magneto-elasticity [START_REF] Danas | Effective response of classical, auxetic and chiral magnetoelastic materials by use of a new variational principle[END_REF][START_REF] Lefèvre | A general result for the magnetoelastic response of isotropic suspensions of iron and ferrofluid particles in rubber, with applications to spherical and cylindrical specimens[END_REF]). Furthermore, with the current 3D printing technology such RVEs can be readily printed and experimentally tested in the lab. Such work is in progress [START_REF] Danas | Numerically-aided 3d printed random isotropic porous materials approaching the hashin-shtrikman bounds[END_REF].

Finally, the present algorithm can be potentially used to generate, 3D print and analyze microstructures where the aspect ratios of the inclusions are provided by external sources (for instance imaging or 3D tomography) via the use of sufficiently defined probability density functions (PDFs) (see for instance the very recent work of [START_REF] Lydżba | Equivalent microstructure problem: Mathematical formulation and numerical solution[END_REF] in this context).

Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant agreement no. 636903-MAGNETO).

A. Extended RSA algorithm for randomly oriented mono-and polydisperse inclusions Nomenclature:

• L 1 , L 2 , L 3 : dimensions of the cuboidal unit cell

Input: L 1 , L 2 , L 3 , f , N ref , N p , ω (r) 1 , ω (r) 2 , χ (r) , f (r) , ξ 1 and ξ 2 . 2: Ouput: v (r) i , a (r) i ,b (r) i , c (r) i , φ (r) i ,θ (r) i ,ψ (
if f - Np-r+2 s=1 f (s) ≤ f i ≤ f (r) then 8:
Compute the dimensions of the inclusion of each phase c for k = 1 to 27 do

(r) i = χ (r) ω (r) 1 ω (r) 2 1/3 R ref , b (r) i = c (r) i /ω (r) 2 and a (r) i = c (r) i /ω (r) 1 9: if ω (r) 1 =

17:

Compute

s 1 = max (a (r) i , b (r) i , c (r) i ) + max (a (r) j , b (r) j , c (r) j ) * ξ 1 18: if v (r) i -v (r) j + m k ≤ s 1 (1 + ξ -1 1 ) then 19:
Compute the minimum distance ∆ 1 using Algorithm 2 Compute s 2 = (max (a for face = 1 to 6 do

28:

Compute the minimum distance ∆ 2 and ∆ 2 using Algorithm 3

Update the volume fraction f

i = f i + (4/3)πa (r) i b (r) i c (r)
i /L x L y L z and set i = i + 1 34: end while 35: Create periodic images of inclusions using instructions given in Section 2.5 Algorithm 2: Minimum distance between two ellipsoids 1: for each pair of inclusions i and j (including its 26 periodic images which centers are defined using the vectors m) do

2:

Define

v 1 = v (r) i and v 2 = v (r) j -m k with k = 1 to 27 3: Define Z 1 = Z (r) i and Z 2 = Z (r) j 4:
Set α to any value in interval [0, 1]

5:

Compute Z α using (2.3.2)

6:

if det(Z α ) > 0 then the two ellipsoids overlap Assign the minimum distance ∆

1 = d(E i , E j ) = 0 12:
Define a maximum number of iterations (maxiter)

13:

while iter <= maxiter do 14:

Compute t 1 and t 2 by solving [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] 15:

Compute the angles between the normals and the vector

x 1 -x 2 16: if ∠(x 2 -x 1 , N 1 ) = 0 and ∠(x 1 -x 2 , N 2 ) = 0 then 17: ∆ 1 = x 1 -x 2 18:
exit jump out of the do loop Compute the two points x 1 and x 2 where N is collinear to n using (15) d(E i , P) using ( 16)

5:
Set the closest point x min to the solution x n which minimizes d(E i , P)

6:

Compute the algebraic measure of the segment defined by the distance ∆ 2 ∆ 2 = n • x min + η 7:

end for 8: end for

B. Convergence analysis

B.1. Convergence of FFT spatial resolution

The convergence analysis in terms of the number of voxels is conducted for three different microstructures made up of monodisperse pores with volume fraction 20% at a fixed realization. Cubic unit cells were generated comprising a number of pores N = 50 for microstructure with spherical pores, and N = 500 for microstructures with oblate (ω = 0.5) and prolate (ω = 2) spheroidal pores. Specifically, each microstructure was discretized with grid of 32 3 , 64 3 , 128 3 and 256 3 voxels. The error committed on the real volume fraction, denoted by δ f , was evaluated for each discretization (Figure 1(a))

δ f = |f -f voxel | f ( 24
)
where f voxel is the volume fraction of pore voxels.

The effective bulk and shear moduli were derived from C using 21. The deviation of C from isotropy, δ iso is evaluated for each discretization using (23)(Figure 1(b)). The choice of a sufficient number of voxels is determined by the convergence of δ f and δ iso to an asymptotic value. Following Ghossein and

Lévesque [START_REF] Ghossein | A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites[END_REF], we introduce a tolerance of 3% on δ f and δ iso as a criterion of convergence. We observe that both variables tend to converge for 128 3 and 256 3 voxels.

B.2. Convergence in terms of number of inclusions/pores

In this section, a convergence analysis in terms of the number of pores was performed on microstructures with prolate (ω = 2) and oblate (ω = 0.5) spheroidal pores with a volume fraction of 20%. At fixed realization, cubic unit cells containing monodisperse randomly distributed and oriented spheroidal pores were generated. Each microstructure was discretized with 128 3 voxels, as discussed in the previous We can observe in Figure 2(a,b) that for a number of spheroidal pores N = 500, the normalized elastic moduli κ and µ do not change significantly. In parallel, a convergence is also obtained on the result for δ iso (Figure 2(c)). Nonetheless, this converged value does not vanish and takes values which can be rather high for oblate voids with very low aspect ratio (Figure 4). This implies that further 510 increase of very elongated inclusions might not necessarily lead to a fully isotropic response. This point deserves further investigation.

B.3. Convergence in terms of realizations

Convergence of the effective moduli is studied in terms of realizations at fixed number of pores, volume fraction and discretization. Monodisperse microstructures containing 50 spherical pores (ω = 1) and 500 spheroidal pores with aspects ratios ω = 0.5, 2, 5 were investigated at volume fraction 20% and discretized with 128 3 voxels. On Figure 3, it can be observed that the normalized effective moduli do not significantly change for different realizations. All reported results in the main text have been averaged over five realizations.

Figure 1 :

 1 Figure 1: Cubic unit-cell and reference ellipsoid geometry, defined by the aspect ratios ω (r) 1 = c (r) /a (r) and ω (r) 2 = c (r) /b (r) and the orientation vectors n (r) i (with i = 1, 2, 3).

1 and ω ref 2 .

 2 The number N ref , the volume fraction f and the reference shape (ω ref 1 , ω ref 2

  generate a new inclusion center position and recheck the overlapping and minimum distance from existing ellipsoids. If ∆ 1 ≥ s 1 proceed to Step 3.

  190

205

  minimum distance d(E 1 , E 2 ) = 0 and terminate, otherwise continue to Initialization.

Figure 2 :

 2 Figure 2: Geometrical representation of iterative steps for finding the minimum distance between two ellipsoids.

1 and x k+1 2 as

 2 the intersections of the segment c k 1 , c k 2 with the ellipsoid surfaces Ω 1 and Ω 2 , so that x k+1 1

210 and only if the normals N k+1 1 and N k+1 2 and the vector x k+1 1 -x k+1 2

 1212 are colinear.

2 ,Step 4 : 1 and c k+1 2 of B 1

 24121 accept the new ellipsoid and move on to the addition of a new ellipsoid, otherwise continue to Step 4. Construct two new ball centers c k+1 and B 2 , respectively, as

Figure 3 Figure 3 :

 33 Figure 3 shows four representative microstructures obtained with the proposed RSA algorithm for a volume fraction f = 20%.

5 Figure 4 :

 54 Figure 4: Deviation from isotropy in monodisperse microstructures compared against that in polydisperse microstructures for different pores aspect ratios ω = 0.1, 0.5, 1, 2 and 5.

Figure 5 :

 5 Figure 5: Comparison between FFT and HSW estimates for monodisperse and polydisperse microstructures with spherical pores. Effective moduli as a function of porosity: (a) normalized bulk modulus κ/κ 1 and (b) shear modulus µ/µ 1

Figure 6 :

 6 Figure 6: Comparison between FFT and HSW estimates for monodisperse and polydisperse microstructures with prolate spheroidal pores. Effective moduli as a function of porosity: (a) normalized bulk modulus κ/κ 1 and (b) shear modulus µ/µ 1

Figure 7 :

 7 Figure 7: Comparison between FFT and HSW estimates for monodisperse and polydisperse microstructures with oblate spheroidal pores. Effective moduli as a function of porosity: (a) normalized bulk modulus κ/κ 1 and (b) shear modulus µ/µ 1 .

Figure 8 :

 8 Figure 8: Influence of the pore aspect ratio ω on the evolution of the effective moduli as a function of the porosity: (a) normalized bulk modulus κ/κ 1 and (b) shear modulus µ/µ 1 .

• 2 :iiAlgorithm 1 :

 21 f : total volume fraction of inclusions • N ref : reference number of inclusions in the unit-cell (equal to number of inclusions when monodisperse) • R ref : Reference radius of a spherical inclusion • N p : number of phases/families of inclusions (different shapes or sizes correspond to different families/phases) aspect ratios of the inclusions belonging to family r • χ (r) : size coefficient defining the allocation of sizes to the different families • f (r) : volume fraction of each family r • ξ 1 , ξ 2 : coefficients controlling the minimum distance between inclusions and between the unit-cell faces and the inclusions • v (r) i (i = 1, 2, 3): position vector of the center of an inclusion belonging to family r : lengths of semi-axes of an inclusion belonging to family r : the three Euler angles defining the orientation of the elliptical inclusion belonging to family r • m = (m 1 , m 2 , m 3 ): total of 27 vectors defining the translational periodicity of the cuboidal unit cell, where m 1 , m 2 and m 3 take the values in all combinations 0, -L 1 , L 1 , 0, -L 2 , L 2 and 0, -L 3 , L 3 . • ∆ 1 , ∆ 2 , ∆ 2 : Distances between inclusion-inclusion and inclusion-cuboidal face • Z (r) : matrix describing the shape and orientation of inclusion belonging to family r Main program 1:

r) i 3 : 4 : 6 :

 346 Compute the characteristic length R ref = 3V f /4πN ref 1/3 , and vectors m Initialize the number of inclusions i = 1, and the actual volume fraction f i = 0 5: while f i ≤ f do for r = 2 to N p do 7:

Initialize the two ball centers c 1 and c 2 10 :c 1 =

 101 v 1 and c 2 = v 2 at iter = 1 11:

19 : else 20 : 2 N n , n = 1 3 :

 1920213 update the center of the balls: c n = x n -γ n Minimum distance from an inclusion surface point to the cell faces 1: for each inclusion i do

4 :

 4 Compute the minimum distance ∆2 = min x∈(x1,x2)

Figure 1 :

 1 Figure 1: Convergence analysis in terms of the number of voxels performed on three different microstructures made up of monodisperse with volume fraction 20% at fixed realization: (a) error on the volume fraction δ f and (b) deviation from isotropy δiso as a function of the number of voxels along each side of the unit-cell (32, 64, 128 and 256).

Figure 2 :

 2 Figure 2: Convergence in terms of the number of spheroidal pores with aspect ratios ω = 2 and ω = 0.5 conducted on monodisperse microstructures. Effective moduli and the deviation from isotropy as a function of the number of pores (50, 100, 200, 300, 400, 500 and 1000): (a) normalized bulk modulus κ/κ1, (b) shear modulus µ/µ1 and (c) isotropy measure.

Figure 3 :

 3 Figure 3: Effective moduli for five different realizations of monodisperse microstructures with 50 spherical pores (ω = 1) and 500 spheroidal pores with aspects ratios ω = 0.5, 2, 5 and porosity f = 20%: (a) normalized bulk modulus κ/κ 1 and (b) shear modulus µ/µ 1 .

Table 1 :

 1 Relative difference between the FE and FFT computations of the effective bulk and shear moduli for three different microstructures

	5	(0.82, 0.86)	(2.2, 1.71)	(1.51, 0.86)
	15	(2.12, 1.75)	(6.32, 5.24)	(8.73, 7.43)

  1 and ω

		(r) 2 = 1 then
	10:	Give Euler angles φ	(r) i , θ	(r) i , and ψ	(r) i
	11:	end if
	12:	(r) i Assign a random center position vector v
	13:	end if
	14:	end for
	15:

for j = 1 to i -1 do 16:

Alternatively, one could check the intersection of two ellipsoids by solving two one-dimensional quadratic inequalities at each iteration, as in[START_REF] Lin | On the distance between two ellipsoids[END_REF], but this is more time consuming.

This algorithm is given in full detail in A.

3. Determination of effective elastic propertiesIn this section, we discuss the computation of the effective elastic properties of periodic porous random particulate materials with an elastic isotropic matrix by considering 3D cubic unit-cells V

In the sequel, the angular brackets • denote the average of a field over the unit-cell.

This approach has been used by many authors. For instance, Kanit et al.[START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] quantitatively characterized the geometry of 3D Voronoï tesselation based on the covariance function related to the covariograms of plane section images of the microstructure while Segurado and Llorca[START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF] have used the radial distribution function of sphere centroids to assess the randomness of the particle distribution, and hence the isotropy of the microsctructure.

Acknowledgments

K.D. would like to thank Prof. P. Suquet for useful comments on improving further the isotropy of the RVEs with very elongated inclusions, which will be included in a future work. The authors wish to acknowledge support from TOTAL S.A. K.D. also acknowledges support from the European Research