Foucault De Francqueville
email: foucault.de-francqueville@polytechnique.edu

Pierre Gilormini

Julie Diani

Representative volume elements for the simulation of isotropic composites highly filled with monosized spheres

Keywords: Micro-mechanics, Particulate media, Finite elements, Representative volume element

A method is proposed for generating reliable representative volume elements (RVEs) that allows reducing the statistical analysis required for the simulation of the mechanical behavior of isotropic composites highly filled with monosized spheres. The method combines (i) an algorithm inspired from molecular dynamics and associated with an analytical equation of state, and (ii) a geometrical analysis using the two-point correlation function and a nearest-neighbor distribution function. A restrictive selection process is defined, which leads to microstructures reasonably close to randomness and isotropy. The pertinence of the proposed generation and selection of RVEs is confirmed by the simulation of their elastic behavior with the finite element method. In particular, it is shown how the selection procedure allows reducing the computational effort required to reach reliable elastic moduli by operating on a limited number of suitable RVEs. The results are in good agreement with the generalized self-consistent model and with original experimental data obtained on a composite where an acrylate matrix was reinforced by sifted glass beads.

Introduction

In order to run numerical simulations, accounting for the local properties and microstructure, of materials highly filled with particles (above 50%) such as propellants, a first step lies in generating representative isotropic random microstructures. In the case of monosized spherical particles, several methods have been designed to generate random microstructures without any one winning over [START_REF] Bargmann | Generation of 3D representative volume elements for heterogeneous materials: a review[END_REF]. Microstructures extractions from tomographies or other imaging techniques may be used for their actual representation of the composites [START_REF] Louis | Application of image analysis for characterization of spatial arrangements of features in microstructure[END_REF][START_REF] Buffière | In situ X-ray tomography studies of microstructural evolution combined with 3D modeling[END_REF][START_REF] Maire | X-ray tomography applied to the characterization of cellular materials. Related finite element modeling problems[END_REF][START_REF] Tariel | 3D Multi-scale segmentation of granular materials[END_REF]. However they might induce some segmentation complications or restrain the applicable numerical boundary conditions. The random sequential addition (RSA) procedure [START_REF] Widom | Random sequential addition of hard spheres to a volume[END_REF], where particles are randomly and sequentially dispersed in a cell without overlapping, is often adopted [START_REF] Böhm | Multi-inclusion unit cell models for metal matrix composites with randomly oriented discontinuous reinforcements[END_REF]; [START_REF] Williams | Random packing of spheres and spherocylinders simulated by mechanical contraction[END_REF], for instance). This scheme is limited to spheres volume fractions below 38% due to the frozen position of each particle [START_REF] Cooper | Random-sequential-packing simulations in three dimensions for spheres[END_REF]. Therefore, alternative approaches have also been developed such as rate-dependent densification algorithm [START_REF] Jodrey | Computer simulation of close random packing of equal spheres[END_REF] or "drop and roll" algorithm [START_REF] Visscher | Random Packing of Equal and Unequal Spheres in Two and Three Dimensions[END_REF] which are rather dedicated to analysis of spheres close packing [START_REF] Torquato | Is random close packing of spheres well defined?[END_REF]. A more recent approach suggests to start from a Poisson distribution of spheres and exploit an optimization formulation to improve particle localization while enforcing a non-overlapping constraint [START_REF] Pathan | A new algorithm to generate representative volume elements of composites with cylindrical or spherical fillers[END_REF]. However the final volume fraction reached hardly exceeds 40%. A Monte Carlo methods allows reaching higher volume fractions [START_REF] Tobochnik | Monte Carlo simulation of hard spheres near random closest packing using spherical boundary conditions[END_REF]. It is based on the random perturbation of the spheres positions, from an initial microstructure that is periodic usually, enforcing the impenetrability of particles. Optimizing and restraining the interval of perturbations ensure the equilibrium of the final microstructure. Such an approach is known as the Metropolis algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. While up to 72% volume fractions of particles may be reached this way, the randomness of highly filled microstructures is questionable [START_REF] Rintoul | Computer simulations of dense hardsphere systems[END_REF], and combined approaches have been developed by mixing the previous methods. The modified RSA scheme, for instance, starts with a low volume fraction RSA microstructure in a cell where particles are submitted to random motions. The particles volume fraction is then increased by incrementally reducing the cell volume [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Bailakanavar | Automated modeling of random inclusion composites[END_REF][START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF]. Alternatively, the cell size can be kept fixed but the particles radius is increased while they undergo a Brownian motion, modeling molecular dynamics [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF][START_REF] Skoge | Packing hyperspheres in high-dimensional euclidean spaces[END_REF][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]. Combination of this latter method with Montecarlo steps (shaking process) has also been envisioned [START_REF] Buryachenko | Modeling of random bimodal structures of composites (application to solid propellants): I. Simulation of random packs[END_REF] but the molecular dynamics approach is retained here because an implementation is directly available where the pressure induced in the cell during the process is computed, which provides very useful information.

Once a microstructure has been generated, its linear elasticity can be estimated to assess its potential as representative volume elements (RVEs) of the composite material [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF][START_REF] Hazanov | Order relationships for boundary condition effects in heterogeneous bodies smaller than the representative volume[END_REF][START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][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][START_REF] Garboczi | Computing elastic moduli on 3-D X-ray computed tomography image stacks[END_REF]. Characterizations of the geometry of the microstructures can also be performed (Rintoul andTorquato, 1996, 1998), and complement the former analysis. The concept of the present approach is to explore the ability of the generating presure, the two-point correlation function and various nearestneighbor distribution functions [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF] to efficiently help selecting the more random and isotropic microstructures [START_REF] Chiu | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF]. Applying such a preselection is arguably uninteresting when dealing only with the linear mechanical responses. Nonetheless, it becomes hard to avoid when running finite elements analyses dealing with nonlinear behavior or with cohesive interfaces [START_REF] Park | Cohesive zone models : a critical review of traction-separation relationship across fracture surfaces[END_REF]. Understanding the statistical discrepancy when considering relatively small RVEs, our goal is to reduce the number of numerical microstructures required to still reach accurate estimates of the material mechanical response, and linear elasticity is used only to validate the preselection. Some recent studies [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF] consider spheres volume fractions over 50% and compare finite element results to the differential model [START_REF] Mclaughlin | A study of the differential scheme for composite materials[END_REF], the Torquato model [START_REF] Torquato | Effective stiffness tensor of composite media: II. Applications to isotropic dispersions[END_REF][START_REF] Torquato | Random heterogeneous media: Microstructure and improved bounds on effective properties[END_REF], or the generalized self-consistent model (GSC) [START_REF] Chiu | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Christensen | A critical evaluation for a class of micromechanics models[END_REF]. Since quantitative comparison with experimental data is scarce in the literature, a model composite made of a rubbery matrix highly filled with sifted glass beads was prepared and tested in uniaxial tension in the present study. This was already the case in the work of [START_REF] Gilormini | Stress-strain response and volume change of a highly filled rubbery composite: experimental measurements and numerical simulations[END_REF], where a finite element analysis (FEA) was applied to a single microstructure generated through a much simpler and less validated procedure than here.

The paper is organized as follows. In the next section, an algorithm gen-erating microstructures is described. The impact of its seed parameters in regards of the resulting microtructures and geometries is presented. Based on the requirement of isotropic mechanical properties, a selection procedure is defined in Section 3 from the finite element analysis of the microstructures, including a careful mesh refinement. Finally, the obtained numerical properties are confronted to homogenization models and to original experimental data in Section 4.

Microstructures highly and randomly filled with monosized spheres

Random cell generation

The [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] approach is applied here to generate microstructures. A detailed example of implementation is availaible in [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], and an other version of the algorithm, developed by S. Torquato and his coworkers [START_REF] Skoge | Packing hyperspheres in high-dimensional euclidean spaces[END_REF], is available for download on the website of this author. First, a classical RSA algorithm is applied to disperse a low volume fraction of monosized spheres randomly in a periodic cubic cell. Then, an initial velocity v is given to each particle according to a Maxwell-Boltzmann distribution:

f (v) = m 2πk B T 3/2 exp -m v.v 2k B T (1)
where all particles have the same normalized mass m = 1, k B is the Boltzmann constant, and the only parameter driving the distribution of velocities is the thermodynamic temperature T . To reach the desired high volume fraction, a growing factor g f act is applied to all spheres while their mass is kept constant, resulting in an unaffected distribution of velocities. Due to the growing spheres size, the probability of collision between moving spheres increases, and an event-driven process is defined [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF], where all velocities are decreased uniformly after a given number of events reducing the risks of two simultaneous collisions that would stop the algorithm. Applying the latter constraint also lowers the risk of strong pressure divergence (this phenomenon will be discussed in Section 2.3 below) that may prevent reaching the expected volume fraction. Examples of configurations obtained during a cell generation are presented in Figure 1 with specification of the main generating parameters: temperature T and growing factor g f act .

1% 5% 14% 30% 55%

Figure 1: Examples of microstructures obtained at various intermediate volume fractions during a microstructure generation with a mixed method combining an initial RSA scheme and particles velocities and growth defined by parameters T = 0.2 and g f act = 0.005.

Geometrical analysis

In accordance with the literature [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF][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][START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF], elementary periodic cells containing 64 spheres of equal diameters were generated [START_REF] Chiu | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF]. A microstructure can be analyzed by computing various geometrical characterization functions, but only few of them lead to usable results here for the relatively small number of spheres [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. Nonetheless, the two-point correlation function S 2 and some nearest-neighbor distribution functions have revealed interesting properties to detect poor randomness. Since reference values have been estimated theoretically by [START_REF] Torquato | Nearest-neighbor distribution functions in many-body systems[END_REF] based on [START_REF] Reiss | Statistical mechanics of rigid spheres[END_REF], the neighbor distribution functions H v that characterizes the distance from every point of the matrix to the closest sphere centers, is chosen.

The S 2 function is a common tool for a first characterization of randomness [START_REF] Lee | Computing overall elastic constants of polydisperse particulate composites from microtomographic data[END_REF][START_REF] Chiu | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF]. In the present context, its evaluation will focus on the reinforcement phase and will therefore be defined as the probability for a vector r with its origin, orientation and length chosen randomly to have both its origin and end located in the particle phase, taking the cell periodicity into account. A classical procedure based on fast Fourier transform, recalled in appendix A [START_REF] Ohser | 3D images of materials structures: processing and analysis[END_REF], was implemented as a Matlab (Matlab, 2016) routine.

Figure 2 illustrates the S 2 values obtained when 55% volume fraction of spheres are either organized as a BCC lattice (Figure 2a) or disordered according to a Percus-Yevick (Percus and Yevick, 1958) random distribution (Figure 2b). Note that the later is a rather close representation of a perfect random distribution at the considered volume fractions [START_REF] Katzav | Random Close Packing and the Hard Sphere Percus-Yevick Theory[END_REF]. Both figures, where r 1 /D, r 2 /D and r 3 /D denote the components of the sampling vector normalized by the spheres diameter D, show the S 2 values computed in the three coordinate planes of one eighth (because of periodicity and parity of the S 2 function) of the cell. While usual analyses are restrained to few directions, uncommon 2D plots are proposed here as extensive samples of the full 3D sets of S 2 values given by the above numerical procedure. Circular lines, which correspond to |r|/D = 1, 2, 3, ..., are also plotted to help the eyes recognizing the isotropy symmetry. In Figure 2a, the crystallographic BCC arrangement in the cell is clearly revealed by the periodic variations of the S 2 function. In contrast, the S 2 values in Figure 2b, computed from the results tabulated by [START_REF] Torquato | Microstructure of two-phase random media. V. the n-point matrix probability functions for impenetrable spheres[END_REF] for the Percus-Yevick distribution, obey a spherical symmetry that reveals isotropy. In Figures 2a and2b, the S 2 value at the origin of the axes is the volume fraction of spheres (η), as it stands for the probability for a point chosen randomly (i.e., a zero-length sampling vector) to lie in a reinforcement. In addition, as soon as sampling vectors growth in Figure 2b, S 2 reaches values close to this fraction squared (η 2 ) which denotes a large range of test vectors with decorrelation at their origin and end, as should apply in a perfectly random distribution. Figure 3b shows the S 2 values obtained from the tomographic analysis of a composite material that has been prepared by randomly mixing a volume fraction of nearly 55% of sifted glass beads in a rubbery matrix (Figure 3a). This experimental two-point correlation function exhibit the same features as displayed in Figure 2b, which is in favor of a fairly random and isotropic distribution of reinforcements for the actual composite.

Pressure evolution

The selected generation algorithm provides with an estimate of the reduced pressure p during spheres growth. The latter is obtained by calculating the Clausius virial function and normalized it with the perfect gas pressure (P gas = nk B T /V , where n is the number of particles and V the volume)

p = 1 + 1 2E c n i=1 r i • F Int i , (2)
where the expression of the Boltzmann kinetic energy E c = 3nk B T /2 has been used and brackets correspond to time averaging. Since the only potential considered during the microstructure generation process is the interaction between colliding particles [START_REF] Skoge | Packing hyperspheres in high-dimensional euclidean spaces[END_REF], the scalar product of sphere position r i with internal force F Int i can be estimated from the momentum X m exchanged by particles collisions during duration ∆t,

p = 1 + X m 2E c ∆t . ( 3
)
The pressure in systems of randomly dispersed monosized hard spheres has already been addressed theoretically, and several equations of states (EOS) have been proposed for its evolution with the particles volume fraction η. The Carnahan-Starling-Kolafa model [START_REF] Boublík | P-V-T behaviour of hard body fluids. Theory and experiment[END_REF], as expressed by [START_REF] Robles | Note: Equation of state and the freezing point in the hard-sphere model[END_REF], is used here:

p = 1 + η + η 2 -2 3 (1 + η) η 3 (1 -η) 3 . ( 4
)
Figure 4 shows the evolution that can be obtained for three representative microstructures, depending on the parameters used in the generation process (see Section 2.1). An atypical confrontation to their corresponding geometrical properties is performed. First, a pressure drop may be observed for low values of the g f act /T ratio. The S 2 values for the corresponding microstructures have similarities with the ordered case of Figure 2a, suggesting that some crystallization happened during the generation process [START_REF] Hoover | Melting transition and communal entropy for hard spheres[END_REF]). This result is confirmed by the nearest-neighbor distribution function, which diverges from the reference defined for a random microstructure. This phenomenon is explained by a drop of the number of reachable random states above a 'freezing point' at about 49.4% volume fraction, when the system tends to order [START_REF] Rintoul | Computer simulations of dense hardsphere systems[END_REF]. Figure 4 illustrates also that a divergent pressure evolution may appear if the g f act /T ratio is too large. The corresponding microstructures have S 2 values close to the distribution of Figure 2b, but their nearest-neighbor distribution functions differ from the reference for a random distribution, suggesting a particles clustering consistent with the sudden increase of the number of collisions. Finally, when intermediate values of g f act /T are used, the pressure follows the EOS closely, and some of the obtained microstructures are found satisfactory in terms of both S 2 and H v . Therefore, the two main parameters of the generation process must be balanced carefully in order to ensure a pressure evolution close to the EOS, and random microstructures may be expected after the geometrical analysis is applied.

Tests of the microstructures isotropy by finite element analysis

FEA framework

The generated microstructures are meshed in order to calculate their elastic properties by FEA and to recognize the less anisotropic microstructures as reasonnable RVEs. The Netgen free meshing software [START_REF] Schberl | Netgen an advancing front 2D/3D-mesh generator based on abstract rules[END_REF][START_REF] Schneider | Automatic threedimensional geometry and mesh generation of periodic representative volume elements for matrix-inclusion composites[END_REF] conveniently generates periodic meshes and lists the periodic displacement constraints to be used by the finite element code Abaqus (2014). Most FEAs of monosized spheres systems previously performed in the literature [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Barello | Comparison between the relaxation spectra obtained from homogenization models and finite elements simulation for the same composite[END_REF] were based on RVEs with limited minimum interparticle distances (δ) for improving the mesh quality. However, [START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF] noticed recently that the mechanical properties are underestimated when large interparticle distances are imposed, typically δ/D > 10 -3 . To avoid interparticle distance restrictions without introducing distorted elements, the mesh is refined where spheres are very close (Figure 5). For any given pair of close particles, a set of circles of latitudes is defined around the pole at the point closest to the other sphere. Several points are uniformly distributed on these circles to reduce the element size ensuring a high quality mesh. Although the refinement induces a significant increase of the number of degree of freedom, the computational duration is shortened. For discretizations with C3D10 tetrahedra, one million elements are generated typically, with five millions degree of freedom. Such a method allows automatic and robust meshing of the microstructures without user input, therefore it can be operated for systematic analyses as easily as the technique based on fast Fourier transform [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] which is classically preferred for its straightforwardness [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][START_REF] Dunant | A critical comparison of several numerical methods for computing eective properties of highly heterogeneous materials[END_REF]. FEA is adopted in the present contribution for its proven flexibility to explore further features of the composite such as finite strain, viscoelasticity and interface debounding [START_REF] Park | Cohesive zone models : a critical review of traction-separation relationship across fracture surfaces[END_REF].

Highly contrasted phases with small-strain linear elastic behaviors are considered in the following numerical applications. The composite consists of an incompressible soft matrix with a Young modulus of 6 MPa filled with hard particles with a Young modulus of 69 GPa and a Poisson's ratio of 0.25. Hybrid elements are used in the incompressible matrix phase. Note that the impact of the matrix properties will also be tested by increasing its stiffness and compressibility in Section 4. A perfect adhesion is assumed at the filler-matrix interface and periodic boundary conditions are applied to the cubic cell. When run on eight cluster nodes (2.5 GHz Intel Xeon E5-2640 CPU, 64 GB RAM), the simulation of one uniaxial tensile test requires approximately one hour. In order to test the isotropy of the composite, each microstructure is submitted to six loadings consisting in either a uniaxial tension or a simple shear along a cell axis, so that three Young moduli (E 1 , E 2 , E 3 ) and three shear modului (µ 1 , µ 2 , µ 3 ) are obtained. The lowest dispersions for each set of moduli are expected to select the less anisotropic microstructures. Moreover, an additional and original criterion for isotropy is proposed when an incompressible soft matrix is considered: because of the incompressible overall behavior of the composite due to an incompressible matrix and comparatively rigid inclusions, in addition to E 1 = E 2 = E 3 and µ 1 = µ 2 = µ 3 , the ratio E/3µ, where E = (E 1 +E 2 +E 3 )/3 and µ = (µ 1 +µ 2 + µ 3 )/3, should be exactly equal to 1 for a perfectly isotropic composite. While perfect isotropy is unrealistic for the present microstructures, this test reject ordered microstructure with cubic symmetry where the equalities between moduli apply although E/3µ differs from 1.

Statistical Analysis

Twenty different microstructures containing a 55% volume fraction of spheres have been built with generation parameters providing cell pressure evolutions in good agreement with the EOS (see Section 2.3). Their elastic behaviors have been computed and are presented in Figure 6 in terms of the three Young moduli and shear moduli normalized by the respective matrix values, and the additional isotropy criterion computed as E/3µ -1 . The whole generation process was repeated for each microstructure, starting from a different RSA distribution, so that the number identifying each microstructure is arbitrary and on te figure each color stand for the properties of one realization. One notes in Figure 6 that all microstructures do not perform equally. For instance, microstructure 6 presents low dispersions of Young moduli and shear moduli, with averages close to the values estimated over the twenty microstructures, and a very good value for the isotropy criterion. In contrast, microstructure 7 exhibits a small dispersion of Young modulus but a large dispersion of shear modulus, with a poor value for the isotropy criterion. Therefore, the following question arises: instead of running FEAs on every generated microstructures, can geometrical analyses be used preliminarily to select fewer microstructures to which FEA is worth being applied? Such a procedure might not be crucial for elastic behaviors, but it might become mandatory when finite strain, interfacial damage, or viscoelasticity, for instance, are considered, because of long computation durations.

Consequently, the geometrical analysis described in 2.2 has been applied to each of the above 20 microstructures, in order to show ability to reject the cases where unsatisfactory mechanical properties could have been expected before running FEA. In order to illustrate the results, we will focus on four microstructures, which are representative of the features that can be Isotropy criterion (c) calculated with the average Young modulus and average shear modulus for each microstructure. Insets correspond to the same data but considering only the 7 selected microstructures obtained. Microstructure number 6 presents a reasonable randomness while numbers 2, 7 and 11 are examples of microstructures that should be discarded before being meshed (See Figure 6). Figure 7 presents the S 2 values for these four microstructures, with a reduced color scale centered on the asymptotic value (η 2 ). Microstructures 2 and 7 can be discarded readily as they do not exhibit any asymptotic trend and as their dispersed peaks are inconsistent with a symmetric evolution of S 2 . In contrast, microstructures 6 and 11 present S 2 values with features similar to the reference displayed in Figure 2b, meaning spherical symmetry of S 2 and convergence to η 2 for long test vectors, which suggests promising randomness and isotropy. As a further test, Figure 8 compares the nearest-neighbor distribution function H v computed on microstructures 6 and 11 to the perfect random reference described in 2.2 [START_REF] Reiss | Statistical mechanics of rigid spheres[END_REF]. This comparison allows rejecting microstructure 11 while microstructure 6 shows the best randomness from the nearest neighbor point of view. Considering all 20 microstructures, we have noticed that every microstructure presenting non-negligible anisotropic shows either a poor S 2 periodicity or unsatisfactory nearest neighbor distribution. Therefore, S 2 and H v characterizations provide efficient tools to discard undesired anisotropic microstructures.

The mechanical properties of the seven microstructures selected (namely, microstructures number 3, 4, 6, 8, 14, 17, and 18) are reported in the insets of Figure 6, highlighting the reduction of the moduli dispersions by a factor of nearly 2 and the exclusion of microstructures with unsatisfactory values for the isotropy criterion. It may be noted that averaging the Young modulus and the shear modulus over the seven final microstructures leads to values (80.9 MPa and 27.9 MPa, respectively) similar to the averages computed over the twenty generated microstructures. Two thirds of the computational effort involved in meshing and finite element simulations would have been avoided by applying the geometrical analysis right after microstructures generation. The remaining fluctuations in the insets of Figure 6 (standard deviations of 10.2% for the Young modulus and 7.3% for the shear modulus) highlight the unfeasibility of generating perfectly isotropic microstructures with such a few number of particles and periodicity justifying the requirement for a statistical approach. Convergence is finally facilitated when, despite studying all the generated realizations, the analysis is retrained to the most representatives microstructures, that can be considered as RVEs.

Representativeness of the selected RVEs

The elastic moduli obtained by FEA on the set of selected RVEs are now compared to the predictions of mean field homogenization models, to FEA results taken from the literature, and to original experimental data measured on a rubbery acrylate network filled with sifted glass beads.

In the case of an incompressible matrix filled with a 55% volume fraction of glass beads, the shear moduli given by the FEA simulations compare very favorably to the predictions of the GSC model [START_REF] Chiu | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Christensen | A critical evaluation for a class of micromechanics models[END_REF] for a wide range of matrix Young modulus E m values, as shown in Figure 9a. This result is consistent with recent literature [START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF][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], even when taking into account the reduced uncertainty obtained with the previous selection process, as illustrated for E m = 6 MPa. It should be noted that the predictions of the Torquato model [START_REF] Torquato | Random heterogeneous media: Microstructure and improved bounds on effective properties[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media: II. Applications to isotropic dispersions[END_REF]) and the differential model [START_REF] Mclaughlin | A study of the differential scheme for composite materials[END_REF] would lay roughly 50% below the FEA results at such a large volume fraction. In order to complete our analysis, much larger matrix stiffnesses were considered while the matrix Poisson's ratio was lowered to 0.4, in order to model a polymer in the glassy state. As shown in Figure 9b, the coherence between numerical results and the GSC model is more questionable in these conditions.

Comparing a model to other models provides with interesting fundamental outcomes but does not prove its relevance for actual materials. While most of the reference studies restrain their confrontation to homogenization models [START_REF] Gusev | Controlled accuracy finite element estimates for the effective stiffness of composites with spherical inclusions[END_REF][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][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF], a composite made of an acrylate polymer network higly filled with sifted glass beads has been prepared and tested. The polymer matrix was obtained by mixing 98% molar mass of methacrylate (MA) with 2% molar mass of polyethylene glycol dimethacrylate (PEGDMA) of molar weight 750 g/mol and adding 2,2-Dimethoxy-2-phenylacetophenone (DMPA) as photoinitiator. This composition was chosen to synthetize a rubbery network at room temperature in order to obtain a composite with a soft and incompressible matrix. Glass beads with diameters ranging from 200 to 250 µm were slowly poured into the acrylate mix, avoiding beads agglomeration and ensuring complete wetting. The compound was poured between two glass plates spaced by a 3 mm teflon frame and maintained vertical during polymer crosslinking in an ultraviolet chamber for 40 minutes. A fairly isotropic distribution of spheres is obtained, as already shown with the tomographic analysis shown in Figure 3. Specimens of 40×12×3 mm 3 were punched at the bottom of the plates to ensure similar glass beads volume fractions in every sample. Moreover, 1 mm-thick plates of pure polymer were also prepared to measure the matrix elastic properties and implement them in the numerical simulations. The composite high volume fraction results from the gravity applying on the beads in a low viscosity uncrosslinked matrix, which prevents any adjustment during the synthesis. Densimetry measures revealed an average density of 52% with a rather small scatter of 1%. Tensile tests at low constant crosshead speed (0.1 mm/min) were carried out, during which the strain was measured by video extensometry (Instron AVE2) and the force by a class 0.5 10 kN cell. Six specimens were tested for the pure acrylate and for the composite, showing good reproducibility. The matrix Young modulus was measured at 5.7 ± 0.5 MPa and the composite Young modulus at 55 ± 1 MPa. A comparison between the actual composite modulus, including uncertainties on modulus and on volume fraction, and the average FEA results obtained from 7 selected RVEs at volumes fractions of 55% and 50% (for which the generation and selection procedures presented above at 55% have been repeated) is shown in Figure 10. Results given by GSC model are also plotted when accounting for the matrix modulus uncertainty. Both this model and FEA applied to the set of selected RVEs predict very satisfactorily the stiffness of an actual incompressible rubbery polymer matrix filled with monosized glass beads at the high density experimentally obtained.

The present approach has been developed and tested only for monosized spherical particles with a high contrast between the rigid inclusions and the rubbery matrix. Further extensions may be envisioned for more complex fillers as dedicated approximation of the equation of state are available mostly based on the scaled particles theory [START_REF] Boublík | Statistical thermodynamics of convex molecule fluids[END_REF][START_REF] Reiss | Statistical geometry in the study of fluids and porous media[END_REF][START_REF] Allen | Hard convex body fluids[END_REF]. This would allow selecting complex microstructures based on the analysis of the generated pressure coupled with the geometric characterization, which is the major originality of the current approach. However particles irregularities require an orientational order analysis in addition to present positioning order examinations. That would be exacerbated when large aspect ratio particles are considered, which could also induce considerations for a nematic thermodynamic state [START_REF] Frenkel | The hard ellipsoid-of-revolution fluid: I. Monte Carlo simulations[END_REF]. Therefore, the natural extension of the present approach to monosized ellipsoid for instance, will require a deep revision of the selection process even if Lubechevsly and Stillinger algorithm for such particles as been already developed (Donev et al., 2005a,b). The pressure-based selection will have to consider nematic states (orientationnal order) in addition to solids state (positional order) and the geometrical characterization will have to include orientation characterization functions. Such extensions are beyond the framework of the present study despite the interests that may arouse for a selection process efficient with complex RVEs.

Another extension would have been analyses of the impact of phases contrast on the relevance of the selection process. Reduction of the contrast has been briefly considered in section 4 without any experimental validation of the numerical estimation of the composite mechanical properties. Studying the inversion of the phase contrast (voids instead of rigid inclusions) might be of particular interest especially for modeling rocks behavior.

Conclusion

In this study, a method was proposed for the efficient generation and selection of reliable RVEs that allows reducing of two thirds the statistical analysis required for the numerical simulation of the mechanical behavior of isotropic composites randomly and highly filled with monosized spheres. The method proceeds in two steps. First, an algorithm inspired from molecular dynamics is applied, where spheres collide and grow simultaneously in a periodic cubic cell, starting from an initial RSA distribution at a low volume fraction. The microstructures obtained at high volume fraction are rejected if a reference equation of state is not followed during this process, which allows adjusting the two parameters involved. In a second step, a geometrical analysis is performed to select microstructures further. The two-point correlation function is computed and allows identifying the features expected from random and isotropic microstructures. Finally, a nearest-neighbor distribution function is computed and compared to a reference in order to perform a final selection among the microstructures.

The relevance of the proposed generation and selection of RVEs has been confirmed by the simulation of their elastic behavior with the finite element method, where the mesh was refined carefully. First, the microstructures selected by the geometrical analysis have been shown to be reasonably close to isotropy in terms of the dispersions of the Young modulus and shear modulus. This was also consistent with the ratio of 3 between these moduli that is expected for isotropic incompressible composites. Moreover, the average moduli obtained from the selected microstructures were found very similar to the averages obtained without applying the geometrical analysis, but with much narrower dispersions. Thus, the proposed selection procedure allows reducing the computational effort required to reach reliable elastic moduli by operating on a limited number of suitable RVEs. This confirms the efficiency of comparing pressure evolution with the EOS and of the considered geometrical functions for a rapid discrimination of the reasonably random and isotropic microstructures. While further morphological characterizations might reduce the number of considered realizations and the properties discrepancies, they will mostly be significantly time consuming or even inoperative due to the cell size.

The numerical results have been confronted to homogenization models, and a very good agreement has been obtained with the generalized selfconsistent model, especially when a soft and incompressible matrix contains a large volume fraction of monosized stiff spheres. Very interestingly, a good agreement has also been found with original experimental results obtained on a model composite where a rubbery acrylate matrix was reinforced by sifted glass beads. This contributes to validate the ability of the proposed methodology to predict the elastic behavior of composites with good precision. Keeping the same guidelines, it may also be adapted to nonlinear behaviors, finite strains or damage at matrix-particles interfaces, for instance. and the calculation can then take full advantage of the remarkable properties of circulant matrices (see for instance [START_REF] Van Loan | Computational frameworks for the fast Fourier transform[END_REF], for basic relations between circulant matrices and Fourier transform, and [START_REF] Yongshe | Calculation of average covariance using fast Fourier transform (FFT)[END_REF] for application to 3D correlation). Practically, we applied the following simple numerical procedure to compute the values of the S 2 function:

(i) The elementary microstructure cubic cell is sampled into a N × N × N matrix with a value of 1 for each voxel located in a spherical reinforcement and 0 in the matrix. (ii) This N × N × N matrix of integers is changed into a matrix of complex numbers, where the real part is equal to the integer value (1 or 0) and the imaginary part is zero. (iii) The discrete 3D Fourier transform of this complex matrix is computed by applying a fast Fourier transform (FFT) algorithm. (iv) All elements of the resulting N × N × N matrix of complex numbers are modified, with their real parts replaced by their moduli squared and their imaginary parts set to zero. (v) The inverse discrete 3D Fourier transform of this modified matrix is computed by FFT, which leads to a N × N × N matrix of complex numbers where the imaginary parts are zero and the real parts are the S 2 values.

It is important to note that a division by N 3 is included in the definition of the discrete 3D Fourier transform used above (but not in the inverse transform, of course), and therefore the division by N 3 involved implicitly in the definition of S 2 (to obtain the average over all locations of the origin of a sampling vector with fixed orientation and length) is performed when the moduli are squared in step (iv) above.

Figure 2 :

 2 Figure 2: Maps of the two-point correlation function S 2 for various microstructures at 55% volume fraction for (a) a perfectly ordered microstructure and (b) a Percus-Yevick random microstructure.

Figure 3 :

 3 Figure 3: Experimental microstructure obtained by X-Ray tomography analysis (a) and the corresponding two-point correlation function (b). The composite is made of a 55% volume fraction of glass beads dispersed randomly in an elastomeric matrix.

Figure 4 :

 4 Figure 4: Three types of reduced-pressure evolutions for various g f act /T ratios, compared to the Carnahan-Starling-Kolafa equation of state.

Figure 5 :

 5 Figure 5: (a) A generated microstructure. (b) Its discretization using the meshing procedure. (c) Enlarged view illustrating the local mesh refinement.

Figure 6 :

 6 Figure 6: Values of the Young modulus (a) and shear modulus (b) obtained in the directions of the three cell axes and normalized by the matrix values for 20 microstructures. Isotropy criterion (c) calculated with the average Young modulus and average shear modulus for each microstructure. Insets correspond to the same data but considering only the 7 selected microstructures

Figure 7 :

 7 Figure 7: S 2 values for typical microstructures with a reduced color scale centered on η 2 . Solid lines are guides for eyes to estimate spherical symmetry. Spherical trend of the left diagrams clearly differ from dissipated peaks of the right ones.

Figure 8 :

 8 Figure 8: Nearest-neighbor distribution function for microstructures 6 (top) and 11 (bottom). The dotted line is the reference function.

Figure 9 :Figure 10 :

 910 Figure 9: Shear moduli given by the FEA simulations and the GSC model when the polymer matrix (Young modulus E m ) of the composite is either rubbery (a) or glassy (b).The volume fraction of glass beads is 55% and the results without error bars were obtained with microstructure number 6.

Acknowledgements

The acquisition of the microtomography device of laboratoire Navier (ENPC) has been made possible thanks to grants from the Region Ile de France (SESAME 2007 program) and CNRS. This work was supported by the ANR under contract number ANR-10-EQPX-37 and the authors thank the Délégation Générale de l'Armement (DGA) and ArianeGroup for their financial supports.

Appendix A

A dense set of S 2 values can be computed by letting the origin and end of the sampling vector scan systematically all the nodes of a thin cubic grid,


Keywords: Equatorial area, Ka-Band, Measurements, Propagation, Satellite communications

ONERA, the French Aerospace Lab, and CNES, the French Space Agency, are currently running a Ka-band propagation experiment at the Guiana Space Centre (CSG) in Kourou (French Guiana). A rain gauge and a beacon receiver able to record the 20.2 GHz beacon signal of the Amazonas 3 satellite have been deployed. The equipment is operational since January 1, 2017 and the duration of the experiment has been set to 3 years. This letter addresses some results of the first year of measurements (from January 2017 to December 2017). The annual and monthly Complementary Cumulative Distribution Functions of rainfall rate and rain attenuation are presented as well as a comparison with the rain attenuation prediction method recommended in ITU-R P.618-13.

I. INTRODUCTION

ITH the increasing demand on High Throughput Satellite (HTS) services, there has been a steady move towards using higher frequencies namely the Ka band (20 GHz-30 GHz) as well as Q/V bands 40 GHz-50 GHz for their higher available bandwidth. However, these frequencies suffer from increased signal attenuation when propagating through the atmosphere. These effects are all the more significant in tropical and equatorial areas due to more adverse weather conditions. Thus in order to better design future telecommunication payloads using EHF frequencies covering tropical and equatorial areas, a thorough study of the propagation impairments needs to be carried out to evaluate the distributions of the static (power margins) and dynamic (fade/interfade durations, site diversity, ...) characteristics of Earth-Space communications at these high frequencies.

To this end, ONERA and the French Ministry of Defence (MoD) designed, at the beginning of the XXI st century, a X. Boulanger and Laurent Castanet are with ONERA / DEMR, Université de Toulouse, F-31055 Toulouse, France (e-mail: Xavier.Boulanger@onera.fr, Laurent.Castanet@onera.fr).

B. Benammar is with CNES DSO/RF/ITP, F-31400 Toulouse, France (email: Bouchra.Benammar@cnes.fr).

Color versions of one or more of the figures in this communication are available online at http://ieeexplore.ieee.org.

Digital Object Identifier ZZZZZ propagation experiment to be held in French Guiana with the 20 and 40 GHz beacons of the satellite STENTOR [START_REF] Lemorton | A new opportunity for EHF propagation experiments : the EXPRESS campaign with the satellite STENTOR[END_REF].

Unfortunately, due to the loss of the satellite during its launch, other opportunities have been investigated to study Earthspace propagation in equatorial areas. Since 2015, ONERA and CNES have been conducting a Ka-band propagation campaign over India in collaboration with ISRO, the Indian Space Agency, using the 20 GHz and 30 GHz beacons of the GSAT-14 satellite [START_REF] Castanet | Characterisation of the propagation channel at Ka and Q/V bands for satellite communication and Earth Observation systems[END_REF], back-up of the original GSAT-4 satellite [START_REF] Dasgupta | A new propagation campaign in tropical areas: the Ka-band Propagation experiment over India with the GSAT-4 satellite[END_REF]. On the other hand, Joanneum Research has been conducting a Ka band propagation experiment in two sites in Malaysia under ESA/ESTEC contract using the 20 GHz beacon of the French Syracuse-3A satellite [START_REF] Cuervo | The Joanneum Research SatCom Ka and Q band campaigns in Europe and Malaysia[END_REF]. Finally, since the end of 2017, CNES and ONERA have been conducting a Ka-band propagation experiment in French Guiana using the 20.2 GHz beacon signal of the Amazonas-3 satellite. Section II of this letter describes the main characteristics of this Ka-band propagation experiment and setup. Some details of the space and ground segments are provided and a short description of the data processing procedures is given. Section III addresses the statistical analyses of rainfall rate and rain attenuation at 20.2 GHz for which annual and monthly Complementary Cumulative Distribution Functions of rainfall rate and rain attenuation are presented. Finally, a brief comparison of the results with the rain attenuation prediction method recommended in ITU-R P.618-13 is highlighted.

II. DESCRIPTION OF THE PROPAGATION EXPERIMENT

A. Space Segment

Amazonas 3, launched into space on February 7, 2013, is a HISPASAT satellite located on the geostationary orbit at 61°W of longitude [5]. A Ka-band beacon signal is available on-board. The exact carrier frequency is 20.1995 GHz with a Right-Hand Circular Polarization (RHCP). The beacon EIRP does not seem to be publically available. Starting from a link budget analysis and the experimental assessment of the dynamic range, it has been roughly estimated to be around 20 dBW.

B. Ground Segment

In December 2016, a CNES Ka-band beacon receiver and a rain gauge have been deployed at the Technical Centre of the Guiana Space Centre (CSG) in Kourou (French Guiana). The exact latitude and longitude of the Amazonas 3 propagation experiment are 5.1713°N and -52.6862°E, respectively (see Fig. 1). The climate normals in Kourou, computed over the The rain gauge (see Fig. 2) is a corrosion resistant dualchamber tipping bucket rain gauge able to record the exact instant of tips (1 second resolution). The main technical characteristics are the following:

• collecting area of 316 cm 2 • 0.1 mm per pulse resolution • 0.015 mm accuracy The beacon receiver (see Fig. 2) has similar characteristics than the receiver deployed in Salon de Provence during the French Ka-band site diversity propagation experiment [START_REF] Boulanger | Small and large scale site diversity experiment at Ka-band in the south of France[END_REF]. Basically, it is composed of a 0.9 m prime focus parabolic antenna with a commercial thermally stabilized Ka-band to L-band Low Noise Block-Converter (linear polarization) performing the down-conversion from 20.2 GHz to L-band. Some adjustments have been added to the external parts of the receiver in order to be more resilient to humid outdoor conditions (anti-rust specific painting and sealing of the electronic boxes). Then, a specifically designed softwarecontrolled L-band to 5 MHz board is used. In this way, the design is easily reusable for different experiments between 19.7 GHz and 20.2 GHz. Finally, the signal is digitized and processed with a frequency search and tracking algorithm to be able to integrate the signal power over a very small frequency bandwidth (typically 50 Hz). It has to be noticed that the receiver is linearly polarized while the polarization of the beacon is RHCP, so a 3 dB loss must be added in the link budget. The clear-sky Signal-to-Noise Ratio (SNR) has then been roughly estimated to be around 34 dB for a 50 Hz integration bandwidth. By fixing the SNR detection threshold to 10 dB, the dynamic range of measurements has been experimentally assessed to reach around 24 dB. The receiver is recording data at 10 Hz sampling rate.

Due to the antenna design and the high elevation angle, water can accumulate inside the main reflector resulting in wetting effect on the measurements. To prevent it, a water discharge hole has been drilled. It enables most of the water to be immediately drained out.

Both the beacon receiver and the rain gauge are operational since January 1, 2017. The elevation and azimuth angles are 78.49° and 238.36° respectively for pointing the Amazonas 3 satellite.

C. Data Processing

The objective of data processing is to extract validated concurrent rain attenuation and rain fall rate time series.

The rain gauge data (tips) are converted into rainfall rate measurements (in mm/h) according to the methodology described in [START_REF]On the processing of tipping bucket rain gauge data for study group 3 experimental database[END_REF]. Beacon data processing consists in extracting a reference level to be subtracted from the original signal power level time series to obtain the desired excess attenuation time series (rain attenuation including a part of cloud attenuation when it is raining) [START_REF] Boulanger | Four Years of Total Attenuation Statistics of Earth-Space Propagation Experiments at Ka-Band in Toulouse[END_REF]. The signal level is first low-pass filtered (raisedcosine filter with 0.025 Hz cut-off frequency) in order to remove tropospheric scintillation. Then, when a rain event is detected, this reference level is linearly interpolated. One of the main characteristics of this propagation experiment is that the elevation angle is very high (78.49°). Consequently, there is a very strong correlation between signal fading and rainfall rate measurements. So, the event detection method is based on the use of the concurrent rainfall rate time series. Finally, the beginning and the end of each rain attenuation event are manually adjusted. An example of processed rain attenuation and rainfall rate time series over one full day (May 2, 2017) is shown in Fig. 3.

III. STATISTICAL ANALYSES

One full year (January 2017 -December 2017) of data has been processed. Daily rainfall rate and rain attenuation time series have been produced and statistical analyses have been performed. The availabilities of valid data coming from the rain gauge and from the beacon receiver are 99.75% and 99.65%, respectively. The availability of strictly concurrent valid data coming from both the rain gauge and the beacon receiver is 99.58%.

The collected annual rainfall amount is 3085.1 mm. May 2017 was the wettest month with 685.2 mm and August was the driest month with only 9.3 mm. Note that a record of precipitation has been set in December with 614.7 mm of precipitation which is more than twice the climate normals [6].

Fig. 4 shows the Complementary Cumulative Distribution Functions (CCDF) of strictly concurrent rainfall rate and rain attenuation at 20.2 GHz. It is worth noting that the dynamic range of the beacon receiver is reached for probabilities of exceedance below 0.1%. Fig. 5 and Fig. 6 show the monthly CCDF of rainfall rate and rain attenuation at 20.2 GHz, respectively. A large monthto-month variability can be noticed and corroborates the presence of the two rainy seasons (from April to May and in December), and the dry season (from August to September). The rainfall rate values exceeded for 0.01 % (namely R 0.01 ) range from 23.9 mm/h to 95.4 mm/h. The experimental value of R 0.01 calculated over the full year of measurements is equal to 85.1 mm/h. For information, the value of R 0.01 predicted by ITU-R P.837-7 [START_REF]Characteristics of precipitation for propagation modelling[END_REF] is equal to 102.1 mm/h. Finally, a brief comparison with the rain attenuation statistics prediction method of Recommendation ITU-R P.618-13 [START_REF]Propagation data and prediction methods required for the design of Earth-space telecommunication systems[END_REF] has been performed. Both the value of R 0.01 predicted by ITU-R P.837-7 (102.1 mm/h) and the experimental one (85.1 mm/h) have been used as inputs of the rain attenuation prediction method. The results are highlighted in Fig. 7 where a fair agreement can be noticed especially when using the experimental value of R 0.01 .

Noting the very unfavorable propagation conditions in Kourou, it tends to prove the burning issue of using Ka-band for satellite communications systems in equatorial areas without adequate Propagation Impairment Mitigation Techniques (PIMT).

IV. CONCLUSION

ONERA and CNES are currently running a Ka-band propagation experiment at the Guiana Space Centre (CSG) in Kourou (French Guiana) with a rain gauge and a beacon receiver able to record the 20.2 GHz beacon signal of the Amazonas 3 satellite. This letter has presented some results of the first year of measurements (from January 2017 to December 2017) and preliminary comparisons with ITU-R recommendations. Statistical analyses have been performed and annual and monthly CCDF of rainfall rate and rain attenuation have been presented. A large month-to-month variability has been noticed. It has also been shown that there is a very good agreement of the results with the rain attenuation prediction method recommended in ITU-R P.618-13.

For frequencies above 20 GHz, gaseous attenuation becomes non-negligible. This is particularly the case in equatorial areas for which the water vapour quantity is very high. This is why, at the beginning of 2018, a GNSS receiver has been installed in order to estimate the integrated water vapour content from the zenith tropospheric delay (ZTD) [START_REF] Siles | An Application of IGS Zenith Tropospheric Delay Data to Propagation Studies: Validation of Radiometric Atmospheric Attenuation[END_REF]. It will be very useful to assess the water vapour attenuation during the Ka-band propagation experiment [START_REF] Siles | An Application of IGS Zenith Tropospheric Delay Data to Propagation Studies: Validation of Radiometric Atmospheric Attenuation[END_REF], which duration is planned to be at least 3 years. Future work will address further propagation analyses such as inter-annual variability of excess attenuation and comparisons with Rec. ITU-R P.678-3, as well as statistics of fade dynamics (fade slope, fade duration, time diversity) and tropospheric scintillation.

Fig. 1 .

 1 Fig. 1. Overview of the measurement site.

Fig. 2 .

 2 Fig. 2. Overview of the equipment.

Fig. 3 .

 3 Fig. 3. Example of processed rain attenuation and rainfall rate time series (May 2, 2017).

Fig. 4 .

 4 Fig. 4. CCDF of concurrent rainfall rate and rain attenuation at 20.2 GHz.

Fig. 5 .

 5 Fig. 5. Monthly CCDFs of rainfall rate.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Monthly CCDFs of rain attenuation at 20.2 GHz.

  Propagation Experiment at Ka-Band in French Guiana: First Year of Measurements This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/LAWP.2018.2887280, IEEE Antennas and Wireless Propagation Letters

Xavier Boulanger, Bouchra Benammar and Laurent Castanet W 1536-1225 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

2 period 1981-2010, are the following [6]: the mean number of rainy days is 190.96 days/year and the mean rainfall amount is 2838.4 mm/year. The precipitations are not homogeneously distributed all over the months and are marked by two rainy seasons (from April to June and from December to January) and a dry season (from August to October). The driest month is September with 32.8 mm of precipitation and the wettest month is May with 503.2 mm of precipitation.

ACKNOWLEDGMENT

This work has been carried out under ONERA internal funding and under CNES R&T study R-S15/TC-0005-066 "Ka-Band Propagation Experiment in French Guiana".

The authors are very grateful to Alexia Lapierre, Coline Morel, Ingrid Epailly, and Christophe Ostapek for their excellent on-site support.


F Lespinasse
email: fleur.lespinasse@cea.fr

E Fillastre
email: eric.fillastre@cea.fr

M Prigniau
email: michael.prigniau@cea.fr

F Vaufrey

Criticality safety control at CEA Paris-Saclay Abstract

est destinée au dépôt

The CEA's organization regarding criticality-safety is the same for each of its nuclear center. First, one or more qualified criticality engineer (QCE) is located directly on nuclear installations containing fissile material. The QCE manages criticality safety in the installations. He develops criticality safety analyses for the installations and participates in the implementation of the resulting rules. Second, a center's criticality engineer (CCE) is in charge of the verification for every modification in each nuclear facilities. He is also in charge to support qualified criticality engineer, and capitalise the knowledge of all the installations. Third, a Criticality Specialist (CS) controls the good functioning of the criticality organization in the center. Furthermore, a criticality safety experts group could assist every criticality engineer of each center's organisation for the criticality calculations and analysis.

This paper details the role of the Criticality Specialist, which is to define an independent opinion on nuclear facilities modifications that could affect the criticality analysis. Its role is also to be an interface between CEA and the French nuclear safety authority for criticality risks in the center. Its action is thus to ensure the quality and consistency of the files sent to the authorities, but without participating in the technical solutions choices.

Then, this paper presents an example of a control carried out on one of the nuclear facilities in Paris-Saclay center. This control is about an implementation of a future operation that may present a criticality risk. The criticality specialist did several further requirements for this operation. In the end, this paper presents how this second line control has been taken into account by the nuclear installation.

Finally, this paper establishes a feedback on the positioning of the criticality specialist regarding the controls performed on nuclear facilities.

L'archive ouverte pluridisciplinaire HAL,


Fabien Hubert

Fabienne Testard

Olivier Spalla

Cetyltrimethylammonium Bromide Silver Bromide Complex as the Capping Agent of Gold Nanorods

A complex between cetyltrimethylammonium bromide (CTAB) surfactant and silver bromide (CTASB) is recognized by NMR and X-ray photoelectron spectroscopy (XPS) to be the entity at the surface of gold nanorods, resulting from an in situ formation in the classical scheme of synthesis. It can thus be introduced directly along with the initial reactants in place of silver(I) salt to produce a new effective synthesis of these objects. Complementary XPS and quartz crystal microbalance (QCM) measurements on macroscopic gold surfaces confirm a strong adsorption of CTASB that is higher than that of CTAB and any other CTAX surfactants. The role of CTASB as a rod inducing agent by surface complexation is stressed.

Introduction

One dimensional metallic nanostructures have aroused great attention due to their future potential applications, expected mainly from their size-dependent optical properties 1-3 or their use as building bricks in bottom-up approaches. 4 Biological applications such as detection or gene delivery are other targets of interest. 5 Electrochemical deposition in templates, 6 electrochemical synthesis in solution, 7 photochemical synthesis, 8 and the seedmediated growth approach 9 are nonexhaustive examples of the numerous existing methods for the synthesis of gold nanorods. Wet chemical methods are the prefered routes due to the costeffective, large-scale production of nanorods, and among them the seed-mediated growth approach is particularly interesting due to the high nanorod yields that can be reached. Indeed, since the original work on seed-mediated synthesis of gold nanorods, 9 extensive studies have been devoted to the increase of the yield, monodispersity, and stability of the particles. In fact, achieving control over the growth toward gold nanorods with high yields is a challenging task. High yields were reached through the use of silver nitrate salt along with the other additives. 10 Shorter aspect ratio rods are obtained compared to those obtained without silver(I) nitrate. However, short aspect ratio rods are interesting in that they present both transverse and longitudinal intense plasmon bands in the visible region of the spectrum, making them promising candidates for sensing or imaging applications. 11,12 Several sources have been advanced as the driving force for anisotropy. It has been suggested 13 that the presence of surfactant bilayers around the nanorods promotes the development of anisotropy through a "zipping" process. Long chain surfactants, such as cetyltrimethylammonium bromide (CTAB), are indeed needed to ensure a high yield of nanorods with good anisotropic ratio, while short chain surfactants lead to a higher percent of spherical byproducts and/or shorter rods. 14 The role of the surfactant counterion was also considered, 15,16 concluding that only brominated long chain surfactants lead to significant yields of rod shaped particles. The final mandatory component to obtain high yield syntheses of rods is silver(I) nitrate. 10 The most commonly cited role of silver(I) is a reaction with CTAB to form the silver bromide compound that would deposit onto some preferential gold rods' facets and thus hinder their growth in that direction. 15 However, it has also been suggested that metallic silver deposits preferentially onto certain faces of the rod, due to the more positive underpotential of deposition of silver(I) on (110) facets, thus poisoning these facets. 17 In this study, we propose an alternative mechanism based on spectroscopic characterizations of a seedless variation of the classical route of synthesis using silver nitrate. We recognize the primordial role of an in situ formed complex of CTAB and silver bromide hereafter named CTASB which can also be prepared ex-situ. 18 Accordingly, an improved synthesis involving the use of CTASB as a cosurfactant along with CTAB and without silver nitrate is performed, which leads to a high yield nanorod synthesis for which the anisotropy promoter is specifically identified.

Results and Discussion

In a first stage, gold nanorods were prepared via a seedless synthesis inspired from Jana: 19 a 10 mL solution containing 3.5 mM aqueous tetrachloroauric acid, 0.2 M CTAB, and 1 mM silver nitrate was constituted. Next, 0.12 mL of 0.4 M ascorbic acid was added under stirring. Finally, 64 µL of an ice-cold aqueous 1.6 mM sodium borohydride solution was added. The * To whom correspondence should be addressed. olivier.spalla@cea.fr.

( reaction, followed using UV-visible spectroscopy, was complete in less than 20 min. Transmission electron microscopy (TEM) analysis of the final particles provides the extensive dimensions for the rods, that is, ∼9 nm in diameter for 32 nm in length, with 90% rods along with very few spherical byproducts. As a reference, gold nanospheres were also prepared by the same synthesis method as above but simply without addition of silver nitrate.

Detailed information about the interaction between the surfactant and gold nanorods can be obtained from NMR measurements. 20 Table 1 presents the chemical shifts of the 1 H spectra for a D 2 O solution of CTAB, as well as D 2 O solutions of capped gold nanospheres and capped gold nanorods, extracted from the initial surfactant solution and washed in heavy water. For all the protons of the main chain, except the R ones, the chemical shifts are identical for all three solutions. On the contrary, the chemical shifts of R-CH 2 and N-(CH3) 3 for the capped gold nanorods, gold nanospheres, and CTAB are different. Indeed, chemical shifts go slightly highfield from CTAB to capped gold nanospheres, probably due to the vicinity of the gold surface to the head's protons (Knight shift 21 ). As for capped gold nanorods, the chemical shifts are even lower, indicative of an additional or a different contribution to the chemical shifts other than that of the Knight shift.

A CTAB silver complex (CTASB) was recently identified 18 in the mixture of CTAB and silver nitrate. Accordingly, the 1 H NMR chemical shifts of CTASB are presented in Table 1. The chemical shifts for a solution of CTASB are indeed identical to the ones obtained for capped gold nanorods. The highfield chemical shifts for R-CH2 and N-(CH3) 3 protons as compared to solutions of CTAB and CTAB capped gold nanospheres can now be attributed to the presence of a silver atom near to the N-polar head protons. This influence is complementary to the simple Knight shift of gold observed on capped gold nanospheres. An enlargement of the peaks for capped nanorods protons from the headgroup is also observed, further supporting the surfactants' binding to the gold surface by their headgroup.

The preceding results suggest that CTASB, a stable complex formed in situ, could be the real promoter of the anisotropy of the nanoparticles. In order to confirm this hypothesis, a new synthesis scheme was designed. CTASB was first prepared ex situ as in ref 18. Its purity was checked by X-ray diffraction (XRD) and elemental analysis, providing its formula C 19 H 42 -NAgBr 2 . A 10 mL solution containing 3.5 mM aqueous tetrachloroauric acid, 0.2 M CTAB, and various concentrations (from 0.2 to 2 mM) of CTASB was then constituted. Next, to each solution, 0.12 mL of 0.4 M ascorbic acid was added under stirring. Finally, 64 µL of an ice-cold aqueous 1.6 mM sodium borohydride solution was added. This new scheme does not make use of silver(I) nitrate or another Ag + salt.

Using UV-visible measurements, the reaction was again followed for each CTASB concentration. The reaction is still complete in less than 20 min and provides nanorods. The XRD pattern of the gold nanorods shows that the particles are crystalline. Peaks are consistent with the bulk face-centered cubic (fcc) pure gold crystal. The respective intensities for the [111], [200], [220], and [311] peaks differ from the theoretical ones for bulk gold, indicating the preferential orientation along [111] and [100] planes in the gold nanorods and more to the point, among the two of them, the predominance of the [111] faces, that is, the side ones of the rods. The aspect ratio of the final particles, extracted from TEM analysis, reaches a maximum of 3.2 for 0.55 mM CTASB (Figure 1). A decrease in CTASB concentration leads to less anisotropic rods and an increase, to similar aspect ratios (ranging from 3.2 to 2.8) (Figure 2), but with the presence of more spheres (40% for 1.6 mM CTASB), as seen in the left inset in Figure 1. The evolution in the aspect ratios is confirmed by UV-visible measurements, though, as stated in ref 22, the estimated value is slightly different from the one obtained by TEM. According to TEM analysis, the highest rod length is reached for 0.4 mM CTASB, while for higher concentrations it decreases, from 37.5 to 22.5 nm (Figure 2). On the other hand, the diameter of the rods constantly decreases with an increase in CTASB, from 15 to 8 nm, suggesting the idea of a surface poisoning by CTASB. In conclusion, the optimum CTASB concentration is 0.55 mM, with more than 90% rods along with few spherical and cubical byproducts.

It was checked that the NMR spectrum of the extracted particles is the same as the one for the CTAB/silver nitrate synthesis.

In order to further characterize the CTASB species on the surface of the nanorods, X-ray photoelectron spectroscopy (XPS) measurements were conducted, on both CTAB/CTASB and CTAB/AgNO 3 nanorod syntheses as well as on pure CTASB. Several drops of the extracted nanoparticle solutions or a CTASB solution were deposited on a silicon wafer, dried under nitrogen, and then inserted in the XPS chamber. The XPS spectrum for gold nanorods showed the Au 4f 7/2 and 4f 5/2 doublet typical for Au(0) (Table 2). The Ag 3d 3/2 and Ag 3d 5/2 peaks are indicative of a Ag(I) oxidation state and stand the same as those for CTASB. The Br 3d 5/2 binding energy for both gold nanorods samples is consistent with Br 3d 5/2 in the pure CTASB spectrum and not with that of a silver bromide sol. The N1s peak can be attributed to N + . 23 Its binding energy at 402.45 eV for both gold nanorod samples is slightly higher than that for pure CTASB due to the gold surface vicinity. In fact, the vicinity of nitrogen with silver mainly explains the shift in energy for CTASB and capped gold nanorods as compared with CTAB whose N + 1s peak is at ∼402 eV. [START_REF] Horr | Colloids Surf[END_REF] Therefore, CTASB appears to be the only species present on the surface of the gold nanorods. The C 1s carbon peak for all samples was split into two contributions: one from the aliphatic carbons at 285 eV and another from the supposedly C-N ones at 286 eV. A Br/Ag ratio of 2 could be drawn from the XPS quantitative analysis of capped nanorods, with this value being in agreement with the chemical formula for CTASB. The XPS spectrum of gold nanospheres, synthesized without silver, was also taken. The binding energies for C 1s, Au 4f, and Br 3d present the same values as those for gold nanorods. However, the N + 1s peak is now at 402.15 eV. This value is slightly higher than that for pure CTAB, indicative of the interaction between CTAB and the gold surface, but differs from the value for gold nanorods and CTASB. It thus confirms the capping of the gold nanorods by CTASB and not CTAB.

All these observations are indicative not only of an in situ formation of CTASB in the CTAB/silver nitrate synthesis but also of a stronger adsorption energy for CTASB than that for CTAB, since after washing of the rods CTASB, though in lower concentration initially, as compared to CTAB, is the only species to be detected via XPS as well as 1 H NMR spectroscopy. In order to understand this exclusion of CTAB from the gold surface, the separate adsorption of these two surfactants onto a bulk gold surface was measured by quartz crystal microbalance (QCM) experiments. The gold covered quartz crystal, purchased from Maxtek, was exposed to a 1 mM CTASB solution and to a 1 mM CTAB solution. The adsorption of CTAB(SB) onto the gold is measured in terms of the frequency shift of the QCM (Figure S5, Supporting Information). The adsorbed amount was calculated with the Sauerbrey equation. [START_REF] Sauerbrey | [END_REF] A total adsorbed mass of 150 ng of CTAB and 340 ng of CTASB was obtained, that is, 400 and 615 mmol, respectively, thus indicating that CTASB adsorbs around 1.5 more than CTAB onto a gold surface. As compared with previous QCM experiments led on several CTAX (X ) NO 3 -, OH -, Cl -, Br -, I -) surfactants 16 where the maximum sorption on a gold surface was obtained with CTAB, CTASB appears to have an even increased adsorption to gold. It implies that when both CTAB (or another of the CTAX surfactants quoted above) and CTASB are present in the solution during the synthesis, the strong adsorption of CTASB implies a strong exclusion of CTAB or other CTAX surfactant from the gold surface and, all the more, that for each tested concentration the total amount of CTASB in solution is sufficient to cover the nanorod surface, even considering one CTASB molecule per gold surface atom. Further adsorption experiments were led on sputtered gold surfaces exposed to different surfactant solutions for 2 h and then thoroughly washed with water. The sputtered gold surfaces were then analyzed by XPS. Two groups could be observed. Gold surfaces exposed to a CTASB solution, CTAB/CTASB, or CTAB/AgNO 3 solutions with a CTAB concentration of 0.2 M and silver concentrations of 1.92 mM and 0.76 mM presented the C 1s, Ag 3d, N 1s, and Br 3d binding energies already reported for CTASB and CTASB capped gold nanorods. A gold surface exposed to a 0.2 M CTAB silverless solution has its N 1s peak at the same position as the one for CTAB gold nanospheres. For surfaces exposed to solutions of CTAB/CTASB or CTAB/AgNO 3 with a CTAB concentration of 0.2 M and a silver concentration of 0.48 mM, the N 1s peak appears to be at an intermediary value, at 402.3 eV, indicative of a mixed surface covering of both CTAB and CTASB. In that case, the concentration of CTASB is too low to dominate the surface adsorption. In addition, the quantification shows that the adsorption of bromide as well as nitrogen is increased in the presence of CTASB or silver nitrate as compared to the adsorption where the surfactant solution solely contains CTAB. Calculated atomic ratios between adsorption in the presence (0.76 or 1.92 mM) and absence of silver are ∼1.8 and 2.9 for nitrogen and bromide, respectively. The final curvature of the rod interface must be due to a subtle balance between a strong adsorption of CTASB which wants to stabilize this Au-CTASB complex and the CTA chain association which promotes a directional growth effect.

The behavior of CTASB was studied by tensiometry, reporting it to have a surfactant behavior. However, no micelle formation could be detected, due to the CTASB low solubility limit in water (inferior to 1 mM). This limit explains why the synthesis cannot be led in a pure CTASB solution and has to make use of CTAB as a cosurfactant. In addition to the demonstrated CTASB surface complexation, a possible scenario is thus that CTAB acts as a chemical stabilizing agent for gold salt, through the formation of a CTA + AuBr 4 -complex. In that case, it prevents a direct reduction in solution to metallic gold by ascorbic acid prior to NaBH 4 addition. Further studies along this direction are in progress.

Conclusion

We have set up a new way of synthesis of gold nanorods with the use of a complex cosurfactant along with CTAB. Silver is now inserted in the solution through this complex form. Spectroscopic characterizations of the final products have evidenced the presence of the complex on the surface of the nanorods. Adsorption studies have confirmed the strong adsorption properties of CTASB on gold bulk surfaces and subsequently stressed its primordial role as a rod inducing agent.

(

  20) Terrill, R. H.; Postlethwaite, T. A.; Chen, C.; Poon, C.; Terzis, A.; Chen, A.; Hutchison, J. E.; Clark, M. R.; Wignall, G.; Londono, J. D.; Superfine, R.; Falvo, M.; Johnson, C. S., Jr.; Samulski, E. T.; Murray, R. W. J. Am. Chem. Soc. 1995, 117, 12537-12548. (21) Knight, W. D. Phys. ReV. 1949, 76, 1259-1260.

Figure 1 .

 1 Figure 1. TEM picture of rods obtained with [CTASB] ) 0.55 mM. Right inset: aspect ratio distribution from TEM images analysis. Left inset: rod fraction (percentage in number of objects) as a function of the CTASB initial concentration from TEM image analysis.

Figure 2 .

 2 Figure 2. Dimensions of the gold nanorods as a function of the CTASB initial concentration.

Table 1 . 1 H NMR Chemical Shits in ppm a

 1

		capped gold		capped gold
	protons	nanorods	CTASB	nanospheres CTAB
	R-CH2	3.31	3.32	3.38	3.44
	N(CH 3)3	3.12	3.12	3.15	3.19
	β-CH2	1.79	1.79	1.78	1.77
	CH 2 main chain	1.30	1.30	1.29	1.30
	CH 3 main chain	0.87	0.87	0.87	0.87
	CH 2 main chain does not appear	1.37	1.37	1.37
	a 1 H NMR spectra in the Supporting Information. Values reported for
	capped nanorods stand the same for the two ways of synthesis (initial presence
	of silver nitrate or CTASB).

Table 2 . XPS Binding Energies in eV a

 2 Nanorod Br 3d, N 1s, and Ag 3d XPS spectra in the Supporting Information. Values for nanorod synthesis in the presence of silver nitrate and CTASB stand the same. b This work.

		capped gold		capped gold
	element	nanorods	CTASB	nanospheres CTAB AgBr sol
	Au 4f7/2	84.0		84.0
	Au 4f 5/2	87.7		87.7
	Ag 3d5/2	367.7	367.7 18			367.7 27
	Ag 3d3/2	373.7	373.8 18			373.7 27
	Br 3d5/2	67.8	67.9 18	67.7	67.3 25	68.3 27
					67.5 26
	N 1s	402.4(5)	402.3 18	402.1(5)	402 24
			402.4 b		402.1 b
	a

Acknowledgment. The authors thank Jean-Marc Verbavatz for TEM characterization, Pascale Jegou for XPS analyses, and Aure ´lie Ithurbide for QCM measurements.

Supporting Information Available: NIR-visible, XPS, and 1 H NMR spectra of gold nanorods; XRD patterns of gold nanorods; and QCM adsorption measurements. This material is available free of charge via the Internet at http://pubs.acs.org. LA801711Q


Minou Nabavi

Olivier Spalla

Bernard Cabane

Surface Chemistry of Nanometric Ceria Particles in Aqueous Dispersions

published or not. The documents may come