Guillaume Noyel
email: guillaume.noyel@cmm.ensmp.fr

Jesus Angulo
email: jesus.angulo@cmm.ensmp.fr

Dominique Jeulin
email: dominique.jeulin@cmm.ensmp.fr

On distances, paths and connections for hyperspectral image segmentation

Keywords: hyperspectral image, connection, quasi-flat zones, η-bounded regions, µ-geodesic balls, top-down aggregation

The present paper introduces the η and µ connections in order to add regional information on λ-flat zones, which only take into account a local information. A top-down approach is considered. First λ-flat zones are built in a way leading to a sub-segmentation. Then a finer segmentation is obtained by computing η-bounded regions and µ-geodesic balls inside the λ-flat zones. The proposed algorithms for the construction of new partitions are based on queues with an ordered selection of seeds using the cumulative distance. η-bounded regions offers a control on the variations of amplitude in the class from a point, called center, and µ-geodesic balls controls the "size" of the class. These results are applied to hyperspectral images.

Introduction

The aim of this paper is to extend and to improve quasi-flat zones-based segmentation. We focus on hyperspectral images to illustrate our developments.

Flat zones, and its generalization, quasi-flat zones or λ-flat zones, initially introduced for scalar functions (i.e., gray-level images) [START_REF] Meyer | The levelings[END_REF] were generalized to color (and multivariate) images [START_REF] Zanoguera | On the implementation of non-separable vector levelings[END_REF] (see def. 2.4 in sect. 2). Classically, λ-flat zones are used as a method to obtain a first image partition (i.e., fine partition). The inconvenience is that the partition typically presents small regions in zones of high gradient (e.g., close to the contours, textured regions, etc.). Several previous works proposed methods to solve it. In [START_REF] Brunner | Iterative area seeded region growing for multichannel image simplification[END_REF], Brunner and Soille proposed a method to solve the over-segmentation of quasi-flat zones. They use an iterative algorithm, on hyperspectral images, based on seeds of area larger than a threshold. Their approach consists in computing an over-segmentation and then merging small regions. Besides, in [START_REF] Salembier | Auto-dual connected operators based on iterative merging algorithms[END_REF], Salembier et al. propose a method to suppress flat regions with area below a given size. They work on the adjacency graph of the flat zones and they define a merging order and criterion between the regions. Moreover, a flat zone, below a given size, cannot be shared by various flat zones. Crespo et al. [START_REF] Crespo | The flat zone approach: a general low-level region merging segmentation method[END_REF] proposed a similar approach for region merging.

However, the main drawback of λ-flat zones for space partition is that they are very sensitive to small variations of the parameter λ. For instance in the example of figure 1, we have only 21 or 1 λ-flat zones depending on a slight variation of λ. In fact, the problem of λ-flat zones lies in their local definition: no regional information is taken into account. For the example of figure 1, with λ = 10 is obtained only one connected region, presenting a local smooth variation but involving a considerable regional rough variation. This problem is more difficult to tackle than the suppression of small regions. The purpose of our study is just to address this issue. We start with an initial partition by λ-flat zones, with a non critical high value of λ, that leads to a sub-segmentation (i.e., large classes in the partition). Then, for each class, we would like to define a second segmentation according to a regional criterion. In fact, two new connections are introduced: (1) ηbounded regions, and (2) µ-geodesic balls; the corresponding algorithms are founded on seed-based region growing inside the λ-flat zones. We show that the obtained reliable segmentations are less critical with respect to the choice of parameters and that these new segmentation approaches are appropriate for hyperspectral images. From a more theoretical viewpoint, the Serra's theory of segmentation [START_REF] Serra | A lattice approach to image segmentation[END_REF] allows us to explain many notions which are considered in this paper.

General notions

In this section some notions necessary for the sequel are reminded.

Hyperspectral images are multivariate discrete functions with typically several tens or hundreds of spectral bands. In a formal way, each pixel of an hyperspectral image is a vector with values in wavelength, in time, or associated with any index j. To each wavelength, time or index corresponds an image in two dimensions, called channel. The number of channels depends on the nature of the specific problem under study (satellite imaging, spectroscopic images, temporal series, etc.).

Definition 2.1 (Hyperspectral image). Let f λ : E → T L (x → f λ (x) = (f λ1 (x), f λ2 (x), . . . , f λ L (x))
), be an hyperspectral image, where: E ⊂ R 2 , T ⊂ R and T L = T × T × . . . × T ; x = x i \ i ∈ {1, 2, . . . , P } is the spatial coordinates of a vector pixel f λ (x i ) (P is the pixels number of E); f λj \ j ∈ {1, 2, . . . , L} is a channel (L is the channels number); f λj (x i ) is the value of vector pixel f λ (x i ) on channel f λj . Definition 2.2 (Spectral distance). A spectral distance is a function d :

T L ×T L → R + which verifies the properties: 1) d(s i , s j ) ≥ 0, 2)d(s i , s j ) = 0 ⇔ s i = s j , 3)d(s i , s j ) = d(s j , s i ), 4)d(s i , s j ) ≤ d(s i , s k )+ d(s k , s j ), for s i , s j , s k ∈ T L
Various metrics distance are useful for hyperspectral points. In this paper, the following two are used:

• Euclidean distance: d E (f λ (x), f λ (y)) = L j=1 (f λj (x) -f λj (y)) 2
• Chi-squared distance:

d χ 2 (f λ (x i ), f λ (x i )) = L j=1 N f .λ j f λ j (xi) fx i . - f λ j (x i ) fx i . 2 with f .λj = P i=1 f λj (x i ), f xi. = J j=1 f λj (x i ) and N = L j=1 P i=1 f λj (x i ).
Definition 2.3 (Path). A path between two points x and y is a chain of points (p 0 , p 1 , . . . , p i , . . . , p n ) ∈ E such as p 0 = x and p n = y, and for all i, (p i , p i+1 ) are neighbours.

Definition 2.4 (Quasi-flat zones or λ-flat zones). Given a distance d :

T L × T L → R + ,
(p i ), f λ (p i+1 )) ≤ λ, with λ ∈ R + .
A path (p 0 , p 1 , . . . , p n ) in an hyperspectral image can be seen as a graph in which the nodes correspond to the points connected by edges along the path. For all i, the edge between the nodes p i and p i+1 is weighted by This definition means that the sum of distances

d (f λ (p i ), f λ (p i+1 )).
n i=0 d (f λ (p i ), f λ (p i+1 )),
along this path, (p 0 , p 1 , . . . , p n ) such as p 0 = x and p n = y, is minimum. It is called the geodesic distance between x and y and noted d geo (x, y).

In order to compute the geodesic path, the Dijkstra's algorithm can be used [START_REF] Dijkstra | A Note on Two Problems in Connection with Graphs[END_REF]. Meanwhile we use an algorithm based on hierarchical queues [START_REF] Vincent | Minimal Path Algorithms for the Robust Detection of Linear Features in Gray Images[END_REF].

For the purposes of segmentation, we need to fix some theoretical notions.

Definition 2.6 (Partition). Let E be an arbitrary

set. A partition D of E is a mapping x → D(x) from E into P(E) such that: (i) for all x ∈ E: x ∈ D(x), (ii) for all x, y ∈ E: D(x) = D(y) or D(x) ∩ D(y) = ∅. D(x) is called the class of the partition of origin x.
The set of partitions of an arbitrary set E is ordered as follows.

Definition 2.7 (Order of partitions). A partition

A is said to be finer (resp. coarser) than a partition B, A ≤ B (resp. A ≥ B), when each class of A is included in a class of B.
This leads to the notion of ordered hierarchy of partitions Π N i=1 D i , such that D i ≤ D i+1 , and even to a complete lattice [START_REF] Serra | A lattice approach to image segmentation[END_REF].

Definition 2.8 (Connection)

. Let E be an arbitrary non empty set. We call connected class or connection C any family in P(E) such that: (0

) ∅ ∈ C, (i) for all x ∈ E, {x} ∈ C, (ii) for each family C i , i ∈ I in C, ∩ i C i = ∅ implies ∪ i C i ∈ C. Any set C of a connected class C is said to be connected.
It is clear that a connection involves a partition, and consequently a segmentation of E. According to [START_REF] Serra | A lattice approach to image segmentation[END_REF], more precise notions than connective criteria (which produce segmentations) can be considered in order to formalize the theory, but this is out of the scope of this paper.

In particular, the λ flat zones can be considered as a connection, λ-flat connection, i.e., λF Z is the partition of the image f λ according to the λ-flat connection.

For multivariate images, where the extrema (i.e., minimum or maximum) are not defined, the vectorial median is a very interesting notion to rank and select the points [START_REF] Astola | Vector Median Filters[END_REF]. Definition 2.9 (Vectorial median). A vectorial median of a set R ⊂ E is any value f λ (k) in the set at point k ∈ R such as:

k = argmin p∈R i/xi∈R d (f λ (p), f λ (x i )) = argmin p∈R δ R (f λ (p)) (1.1)
In order to compute the vectorial median f λ (k), k ∈ R, the cumulative distance δ R (f λ (p)) has to be evaluated for all p ∈ R. Therefore all points p are sorted, in ascending order, on the cumulative distance. The resulting list of ordered points is called ascending ordered list based on the cumulative distance. The first element of the list is the vectorial median (the last element is considered as the anti-median).

η-bounded regions

One of the main idea is to understand that on quasi-flat zones the distance or slope between two neighbouring points must be inferior to the parameter λ. We can establish a comparison with a hiker, from one point who will only deal with the local slope and not on the cumulative difference in altitude on the λ-flat zones. To consider this limitation we propose a new kind of regional zones: η-bounded regions, according to the following connection. Definition 3.1 (η-bounded connection). Given an hyperspectral image f λ (x) and its initial partition based on λ-flat zones, λF Z, where λF Z i is the connected class i and R i ⊆ E (with cardinal K) is the set of points p k , k = 0, 1, 2, ..., K -1, that belongs to the class i. Let p 0 be a point of R i , named the center of class i, and let η ∈ R + be a positive value. A point p k belongs to the η-connected component centered at p 0 , denoted ηBR p0 i , if and only if d(f λ (p 0 ), f λ (p k )) ≤ η and p 0 and p k are connected.

For each class λF Z i the method is iterated using different centers p j (j = 0, 1, • • • J) until all the space of the λ-connected component is segmented:

∪ J j=0 ηBR pj i = λF Z i , ∩ J j=0 ηBR pj i = ∅,
where the η-bounded regions are also connected. Each center p j belongs to λF Z i \ ∪ j-1 l=0 ηBR p l i . The new image partition associated to η-bounded connection is denoted ηBR. It is evident that this second-class connection is contained in the λF Z initial connection, i.e., ηBR(x) ≤ λF Z(x). As shown by Serra [START_REF] Serra | A lattice approach to image segmentation[END_REF], a center or seed p j is needed to guarantee the connectiveness, thus being precise, the region is bounded with respect to the center. The η-connection is a generalization of the jump connection [START_REF] Serra | Set connections and discrete filtering[END_REF] with the difference that working on hyperspectral images, the seeds p j cannot be the minima or maxima. We propose to compute the median value as initial seed.

Note that the method can be also applied on the space E of the initial image f λ , without considering an initial λF Z (which is equivalent to take a value of λ equal to the maximal image distance range). The advantage of our approach is that we have now a control of the local variation, limited by λ, and the regional variation, bounded by η. Moreover, the computation of seeds is more coherent when working on relatively homogenous regions.

In practice, for all the K points on each λF Z i , the ascending ordered list based on cumulative distance δ i is computed. Then, the first element of the list, i.e. the vectorial median, is used as first seed k of the η-bounded region ηBR 1 i . The distance, from the seed to each point p ∈ λF Z i , such as p ∈ N eigh(q), q ∈ ηBR 1 i , is measured. If this distance is less than η then p is added to ηBR 1 i and removed from δ i . For all the others points q ∈ λF Z i \ ηBR 1 i , the first point of the list δ i is the seed of the second η-bounded region ηBR 2 i . Then we iterate the process until all the points of the λ-flat zone λF Z i are in an η-bounded region (algorithm 1).

Algorithm 1 η-bounded regions

Given a distance d, the λ-flat zones, λF Z, of an image f λ , δ a list of cumulative distance, Q a queue, imOut an output scalar image Initialize the value of η currentlabel ← 0 for all λF Z ∈ f λ do for all point p ∈ λF Z do distance ← P q∈λF Z d (p, q) δ ← add the pair (p, distance) end for Ascending sort on parameter distance of δ while δ is not empty do

k ← first point of δ Q ← push(k) while Q is not empty do p ← pop(Q) imOut(p) ← currentlabel
Remove q and its distance in δ for all q ∈ N eigh(p) and q ∈ δ do if d (k, q) ≤ η then Q ← push(q) end if end for end while currentlabel ← currentlabel + 1 end while end for Using again the mountain comparison, this can be compared to a hiker starting from a point with a walk restricted to a ball of diameter 2 × η centered on the starting point. He cannot go upper or lower than this boundary. The points to be reached by the hiker, given the previous conditions, are connected and constitute an η-bounded region.

To understand the effect of these regions, they are applied to the image "tooth saw" with an Euclidean distance d E (fig. 2). For the sake of simplicity, in this image only the grey levels of the first channel have a shape of "tooth saw" (fig. 1 (b)), the others being constant.

We notice from the figure that the smallest the parameter η, the smallest the area of η-bounded regions. Besides, η-bounded regions are very sensitive to the peaks when working on scalar functions. In fact, if η is less than the difference of altitude between the seed and the peak, the points in between are in the same η-bounded region. However, if η is larger than the difference of altitude between the seed and the peak, the points behind the peak, for which the difference of altitude from the seed is less than η, are in the same η-bounded region (fig. 3).

(a) (b) (c) (d) (e) (f) (g) (h) (i)

µ-geodesic balls

As for η-bounded regions, we have created µ-geodesic balls, µGB, to improve the λ-flat zones, but now introducing a control of the dimension of the zone. First of all, λ-flat zones are built. Then, the cumulative difference in altitude is measured from a starting point, the seed, in each λ-flat zone. From this point a geodesic ball of radius µ is computed. This zone is a µ-geodesic ball. It corresponds to the maximum of steps (points) that can be reached by a hiker starting from the seed inside a λ-flat zone, for a given cumulative altitude. µ-geodesic balls are defined using the following connection:

Definition 4.1 (µ-geodesic connection). Given an hyperspectral image f λ (x) and its initial partition based on λ-flat zones, λF Z, where λF Z i is the connected class i and R i ⊆ E (with cardinal K) is the set of points p k , k = 0, 1, 2, ..., K -1, that belongs to the class i. Let p 0 be a point of R i , named the center of class i, and let µ ∈ R + be a positive value. A point p k belongs to the µ-connected component centered at p 0 , denoted µGB p0 i if and only if

d geo (f λ (p 0 ), f λ (p k )) ≤ µ.
It is important to notice, that the geodesic paths imposed the connectivity to the µ-geodesic ball. Formally, for each class λF Z i the method is iterated using different centers p j (j = 0, 1, • • • J) until the full segmentation of the λ-connected component.

Algorithm 2 µ-geodesic balls

Given a distance d, the λ-flat zones, λF Z, of an image f λ , δ a list of cumulative distance, Q a queue, imOut an output scalar image Initialize the value of µ currentlabel ← 0 for all λF Z ∈ f λ do for all point p ∈ λF Z do distance ← P q∈λF Z d (p, q) δ ← add the pair (p, distance) end for Ascending sort on parameter distance of δ while δ is not empty do k ← first point of δ for all point p inside the geodesic ball of center k and radius µ inside the λF Z do imOut(p) ← currentlabel Remove p and its distance in δ end for currentlabel ← currentlabel + 1 end while end for

The new image partition associated to µ-geodesic connection is denoted µGB. As for η-bounded connection, this second-class connection is contained in the λF Z initial connection, i.e., µGB(x) ≤ λF Z(x). The advantage of this approach is that we have now a regional control of the "geodesic size" of the classes by measuring the geodesic distance, limited by µ inside the local variation limited by λ. In practice, µ-geodesic balls are built as η-bounded regions, except that from each seed the geodesic ball is computed inside the λF Z.

µ-geodesic balls are computed with an Euclidean distance d E for the image "tooth saw" (fig. 2). We notice that the smaller is µ, the smaller the area of µ-geodesic balls is. Besides, these balls are not very sensitive to the peaks. In fact, if µ is less than the difference of altitude between the seed and the peak, the points in between are in the same µ-geodesic ball. However, if µ is larger than the difference of altitude between the seed and the peak, only the points behind the peak, for which the cumulative altitude from the seed is less than µ, are in the same µ-geodesic ball (fig. 3).

Results and discussions

In order to illustrate our results on real images, we extracted η-bounded regions and µ-geodesic balls in the image "woman face" of size 45 × 76 × 61 pixels (fig. 4). The channels are acquired between 400 nm and 700 nm with a step of 5 nm. Moreover, to reduce the number of flat zones in this image, a morphological leveling is applied on each channel, with markers obtained by an ASF (Alternate Sequential Filter) of size 1. Then, λ-flat zones are computed. Besides, the computation time in our current implementation with Python is moderate, i.e. a few minutes. However, note that queues algorithms in C++ are very fast. It is important to choose an appropriate distance with respect to the space of the image. In the spectral initial image of "woman face" we choose the Chi-squared distance. For η-bounded regions and µ-geodesic balls, the number of zones minus the number of λ-flat zones is measured. In the figure 6, we notice that the number of zones decrease with the parameters η or µ. From figures, we notice that η-bounded regions are less sensitive than µ-geodesic balls to small variations on distances between points. However, the area of µ-geodesic balls is more controlled than the area of η-bounded regions. Therefore, η-bounded regions are better to find the details and µ-geodesic balls are better to build smoother zones. Besides, in order to evaluate the influence of choosing the vectorial median as a reference seed, we have tested the use of the reverse order for the ascending ordered list based on cumulative distance. In fact, this order corresponds to the vectorial anti-median (fig. 7). Comparing these figure to the zones obtained with the median seed (fig. 5), we notice that almost the same zones are obtained. Consequently, η-bounded regions and µ-geodesic balls have a small dependence to the chosen seeds. Moreover, the same segmentations can be obtained in factor space using an Euclidian distance because it is equivalent to Chi-squared distance in image space. We have also computed it, keeping three factorial axes with relative inertia: 87.2 %, 10.2 % and 1.5 %. By reducing the volume of data, the computation is more efficient on 3 channels than on 61.

(a) (b) (c) (d) (e) (f) (g) (h) (i)
• • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • • 0

Conclusion and Perspectives

We have presented two new connected zones: η-bounded regions and µgeodesic balls. They improve the λ-flat zones, which deals only with local information, by introducing regional information. Moreover, these new connections are of second order because they are built, and included, in the λ-flat zones which are already connected. The approach consists in selecting a sufficiently high parameter λ to obtain first a sub-segmentation. Then, ηbounded regions or µ-geodesic balls are built, leading to a segmentation by a top down aggregation. The η-bounded regions introduce a parameter controlling the variations of distance amplitude in the λ-flat zones, meanwhile µ-geodesic balls introduce a parameter to control the size by controlling the cumulative amplitude inside the λ-flat zones.

Besides, these two second order connections produce pyramids of partitions with a decreasing number of regions when the value of η or µ is increased, until the partition associated to the last level is equal to the partition defined by the λ-flat zones. However, it is important to notice that this pyramid is not an ordered hierarchy in the meaning that the classes are not ordered by increasing level.

Furthermore, we have proposed algorithms for the construction of partitions associated to both new types connections, η and µ, which are based on queues with an ordered selection of seeds using the cumulative distance.

About the perspectives, we notice that the proposed method does not solve the problem of small classes of the initial partition of the λ-flat zones. However, we can combine our method with the approaches aggregating smaller regions to these of larger area [START_REF] Brunner | Iterative area seeded region growing for multichannel image simplification[END_REF][START_REF] Crespo | The flat zone approach: a general low-level region merging segmentation method[END_REF][START_REF] Salembier | Auto-dual connected operators based on iterative merging algorithms[END_REF]. For the future, we are thinking on more advanced methods in order to select the seeds, and to determine locally, for each λ class, the adapted value of η or µ.

Figure 1 .

 1 Figure 1. λ-flat zones of image "tooth saw" (21×21 pixels). (a) Image. (b) Row profile. (c) λ = 9.9 (21 zones). (d) λ = 10 (1 zone).

Definition 2 . 5 (

 25 Geodesic path). The geodesic path between two points x and y in E is the path of minimum weight.

Figure 2 .Figure 3 .

 23 Figure 2. η-bounded regions and µ-geodesic balls of image "tooth saw" (21 × 21 × 4 pixels) for λ = 10. Seeds are marked with a white (or black) point, when they are not trivial. (a) λF Z. λ = 10. (b) η = 0. (c) η = 10. (d) η = 20. (e) η = 30. (f) µ = 0. (g) µ = 20. (h) µ = 40. (i) µ = 100.

Figure 4 .

 4 Figure 4. Channels of image "woman face" f λ (45 × 76 × 61 pixels) and λ-flat zones. (a) f λ 30 . (b) f λ 45 . (c) f λ 61 . (d) λ = 0.003. (e) λ = 0.006. (Source: Spectral Database, University of Joensuu Color Group, http://spectral.joensuu.fi/)

Figure 5 .

 5 Figure 5. η-bounded regions and µ-geodesic balls of image "woman face" for λ = 0.005 (Chi-squared distance d χ 2 ). (a) λF Z. λ = 0.005. (b) η = 0.007. (c) η = 0.009. (d) η = 0.011. (e) η = 0.02. (f) µ = 0.01. (g) µ = 0.02. (h) µ = 0.03. (i) µ = 0.05.

Figure 6 .

 6 Figure 6. Variations of the number of ηBR or µGB minus the number of λF Z versus the parameter η or µ in image "woman face" for λ = 0.005 (Chi-squared distance d χ 2 ). (a) ηBR. (b) µGB.

Figure 7 .

 7 Figure 7. η-bounded regions and µ-geodesic balls with an anti-median seed in image "woman face" for λ = 0.005 (Chi-squared distance d χ 2 ). (a) λF Z, λ = 0.005. (b) η = 0.009. (c) η = 0.011. (d) µ = 0.02. (e) µ = 0.03.

  two points x, y ∈ E belongs to the same quasi-flat zone of an hyperspectral image f λ if and only if there is a path (p 0 , p 1 , . . . , p n ) ∈ E n such as p 0 = x and p n = y and, if, for all i, (p i , p i+1 ) ∈ E 2 are neighbours and d (f λ


B Figliuzzi

D Jeulin

M Faessel

F Willot

M Koishi

N Kowatari

Modelling the microstructure and the viscoelastic behaviour of carbon black filled rubber materials from 3D Simulations

Technische Mechanik 32(1-2), 22-46 (2016). "Postprint" version with larger figures, text and results unchanged.

Volume fraction and spatial repartition of fillers impact the physical properties of rubber. Extended percolating networks of nano-sized fillers significantly modify the macroscopic mechanical properties of rubbers. Random models that describe the multiscale microstructure of rubber and efficient Fourier-based numerical algorithms are combined to predict the material's mechanical properties. From TEM image analysis, various types of multiscale models were proposed and validated, accounting for the non-homogeneous distribution of fillers: in the present work, aggregates are located outside of an exclusion polymer simulated by two families of random models. The first model generates the exclusion polymer by a Boolean model of spheres. In the second model, the exclusion polymer is a mosaic model built from a Johnson-Mehl tessellation. Here the exclusion polymer and the polymer containing the filler show a similar morphology, contrary to the Boolean model. Aggregates are then described as the intersection of a Boolean model of spheres and of the complementary of the exclusion polymer. Carbon black particles are simulated by a Cox model of spheres in the aggregates. The models rely on a limited number of parameters fitted from experimental covariance and cumulative granulometry. The influence of the model parameters on percolation properties of the models is studied numerically from 3D simulations. Finally, a novel Fourier-based algorithm is proposed to estimate the viscoelastic properties of linear heterogeneous media, in the harmonic regime. The method is compared to analytical results and to a different, time-discretized FFT scheme. As shown in this work, the proposed numerical method is efficient for computing the viscoelastic response of microstructures containing rubbers and fillers.

Introduction

Mechanical properties of rubbers are often improved by including nanoscopic fillers in their composition. Usually, silica and carbon black materials are used for this purpose. The geometrical properties of the fillers are an important feature of the global material. In particular, their volume fraction and spatial distribution largely determine the physical properties of the rubber. Extended percolating networks of fillers can for instance significantly enhance the mechanical properties of rubbers on macroscopic scales.

This work investigates the influence of the fillers geometry on the viscoelastic properties of rubber. An important step toward the completion of this research is the development of random models describing their complex microstructure. Following a now well-established method in material science as in [START_REF] Jeulin | Modèles morphologiques de structures aléatoires et de changement d'échelle[END_REF], random realizations of these models can in turn be used to perform an extensive study of the mechanical properties of the rubber.

As reviewed by [START_REF] Jean | A multiscale microstructure model of carbon black distribution in rubber[END_REF], several models have been developed in the past 20 years to describe the microstructure of rubbers. [START_REF] Bergström | Constitutive modeling of the large strain time-dependent behavior of elastomers[END_REF] developed a model describing the carbon aggregates as squares or dodecahedra in an elastomeric volume. Even if their model only accounts for the scale of aggregates, they were able to conduct extensive finite elements calculations aimed at describing the large deformation behavior of the rubber. Relying on a similar approach, [START_REF] Jha | Micro-structural finite element modelling of the stifness of filled elastomers: the effect of filler number,shape and position in the rubber matrix[END_REF] developed a cuberbille model to describe carbon aggregates in a rubber matrix. In their model, they notably added a third phase accounting for a layer of bound rubber surrounding the aggregates. [START_REF] Naito | 3D modeling and simulation of micro to macroscopic deformation behavior of filled rubber[END_REF] improved the microstructure description by considering that aggregates can be represented as the union of carbon black spherical-shaped particles, eventually opening the way to the development of more accurate models. Carbon black composites have finally been widely investigated at the Center for Mathematical Morphology by [START_REF] Jeulin | Morphological modeling of random composites[END_REF], [START_REF] Savary | Morphological analysis of carbon-polymer composite materials from thick sections[END_REF], [START_REF] Delarue | 3D morphological analysis of composite materials with aggregates of spherical inclusions[END_REF] and [START_REF] Moreaud | Multi-scale simulation of random spheres aggregates: Application to nanocomposites[END_REF].

In this paper, we first introduce the studied material nanostructure from TEM micrographs. After a morphological segmentation of images, measurements are performed to quantitatively characterize the microstructure. Then, two random set models to reproduce the multiscale microstructure of rubbers are developed. They are identified from the available experimental data obtained on TEM images and on simulations. In Section 4, we focus on the issue of percolation, which is of particular interest regarding the microstructure optimization. In the last part, viscoelastic properties of the material are predicted by means of Fast Fourier Transforms performed on simulated 3D images generated from the random model.

TEM Images analysis 2.1 Material

This work concerns a rubber material reinforced by carbon black particles as fillers, that are embedded in a matrix of rubber and can be geometrically well-approximated by spheres of diameter D 30 -40 nm, as illustrated in Figure 1a. The volume fraction of the fillers is around 14 percents in the rubber. A key feature of the material is that the rubber matrix is constituted of two distinct polymers: an exclusion polymer, which cannot contain any filler, and its complementary polymer, which can contain fillers. Fillers tend to agglomerate together within the rubber matrix. The filler particles indeed exhibit a turbostratic structure which enable them to merge and create aggregates [START_REF] Jean | A multiscale microstructure model of carbon black distribution in rubber[END_REF].

TEM images segmentation

To analyze the rubbers microstructure, we dispose of a set of 50 TEM micrographs of size 1024 by 1024 pixels with resolution 2.13 nm per pixel. These images were obtained in the Yokohama Rubber Research Center. The slices of material probed by the microscope have a thickness around 40 nm. The first task is to segment these images, in order to identify the spatial distribution of the fillers. The segmentation procedure and its results are presented in the first two paragraphs of this section. The next step is to extract morphological characteristics from the segmented TEM images, as described in the third paragraph.

Algorithm The segmentation process aims at extracting the carbon black phase from the 8 bits TEM images. White regions can be observed in the images, which correspond to aggregates prints that have been taken off by the blade of the microtome during the samples preparation. These prints have to be included in the set of regions to be extracted. The main difficulty here is the presence of an illumination gradient and of contrast difference. Therefore, a simple thresholding is not sufficient to discriminate the black visible aggregates from the rubber matrix.

The first step of the segmentation consists in denoising the TEM images. The denoising is performed using an alternate sequential filter, a classical filter in mathematical morphology [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. The second step of the segmentation process is the extraction of the white aggregates prints. This task is easily performed by thresholding the filtered image with a grey level higher than 170. Finally, the prints with area smaller than 3 pixels are removed from the thresholded image.

The third step of the segmentation process is the extraction of the visible black carbon aggregates. First, the white prints of aggregates are extracted by an opening with a disc of radius 50 pixels to the filtered image, followed by a geodesic reconstruction. Once the white prints have been removed, a black top hat is performed, to account for contrast difference in the image and to detect the darkest zones. A closing of size 20 pixels is used in this purpose. The top hat introduced by [START_REF] Meyer | Iterative image transformations for an automatic screening of cervical smears[END_REF] is a morphological operator which returns the difference between an original image and that same image after opening by some structuring element. This operator is classically used to remove the large scale contrast in an image. A black top hat is simply a top hat operator applied on the inverse image. At this point, a threshold on the resulting image extracts the visible black aggregates, with a grey level around 15. This threshold level is specified to recover an average apparent area fraction of fillers close to 28 percents for the available TEM images. This apparent area fraction results from the projection of the particles on the image during the acquisition by TEM on slices with a thickness approximately twice the radius of a black carbon particle. Finally, the aggregates with area lower than 3 pixels are removed by an area opening applied to the thresholded image. The results of the segmentation procedure are shown step by step for an initial TEM image in Figure 1.

Morphological analysis An efficient way to keep track of the information embedded in the segmented images is to rely on a morphological characterization of the material. In this work, use is made of the covariance, to measure scales in the images, and of size distributions of the filler by opening transformations.

The covariance of a random set A is the probability that x and x + h both belong to A: where h is some vector of R n . For a stationary random set, the covariance is a function of the distance h only:

C(x, x + h) = P {x ∈ A, x + h ∈ A}, (1)
C(x, x + h) = C(h). (2)
The covariance C(h) provides useful information about the spatial arrangement of the random set A and can account for the presence of several scales in the studied set or for periodicity. By definition, C(0) simply corresponds to the volume fraction p of the set A. In addition, the covariance C(h) reaches a sill at the distance or range d. At this distance, we have

C A (d) = p 2 . ( 3
)
The average covariance measured on segmented binary images is given in Figure 2a.

Another important statistical characteristic is the granulometry of the material. The axiomatic of granulometries was first formulated in [START_REF] Matheron | Random sets and integral geometry[END_REF]. In this work, we will consider a granulometry relying on a family of openings by disks. Let K be a convex set. We consider the family {K λ , λ > 0}, where K λ = λK. The opening operator (the symbols and ⊕ representing the morphological erosion and dilation)

Φ λ (A) = (A Ǩλ ) ⊕ K λ , (4)
defined for all closed set A, is a granulometry. Therefore, the granulometry by openings provides a convenient way to characterize the size distribution of the aggregates in the material. The cumulative granulometry curve of the carbon black filler, deduced from the changes of volume fraction measurement with the size λ, is shown in Figure 2b.

Using the covariance and the granulometry curves, we can select a proper set of parameters for the identification of random sets models and control their validity, as explained in the next section.

Random models of multiscale microstructures

As suggested by the intrinsic nature of the material and as observed after the segmentation, the fillers tend to regroup in the material to create aggregates. As a consequence, we have to rely on a two-scale model to properly describe the microstructure [START_REF] Jeulin | Multi scale random models of complex microstructures[END_REF][START_REF] Jeulin | Morphology and effective properties of multi-scale random sets: A review[END_REF]. The first scale corresponds to the aggregates, while the second one describes more specifically the single particles inside the aggregates.

In addition, several alternatives can be considered to locate the aggregates outside of the exclusion polymer. scale. In this work, two distinct models are used in this purpose. The first one relies on a Boolean model of spheres of constant radius to model the exclusion polymer distribution, the second one on Johnson-Mehl tesselations. In both cases, the location of the centers of carbon black spheres are located according to a Poisson point process restricted to the permitted zones, namely outside of the exclusion polymer and inside the aggregates, generating a Cox Boolean model of spheres [START_REF] Jeulin | Multi scale random models of complex microstructures[END_REF][START_REF] Jeulin | Morphology and effective properties of multi-scale random sets: A review[END_REF]. Both models are presented and discussed in this section. All simulations presented here are realized using the software VtkSim, developed in [START_REF] Faessel | VtkSim software[END_REF][START_REF] Faessel | 3D multiscale vectorial simulations of random models[END_REF] at the Center for Mathematical Morphology. The method used for the implementation is briefly introduced later in this section.

Exclusion polymer based on a Boolean model of spheres

Theoretical model

The first model relies on the Boolean model [START_REF] Matheron | Random sets and integral geometry[END_REF]. As eluded to earlier, the studied material contains two kinds of polymers:

• an exclusion polymer P e with volume fraction v e ,

• a second polymer P i with volume fraction v i that can contain fillers.

The exclusion polymer is modelled by a Boolean model of spheres of constant radius R e (Figure 3a) The model depends on two parameters: the intensity θ e of the underlying Poisson point process and the exclusion radius R e . These two parameters can be related to the volume fraction of exclusion polymer to yield:

1 -v e = exp -θ e 4π 3 R 3 e . (5)
Therefore, if the volume fraction of the exclusion polymer is known, the exclusion polymer modelling is controlled by one single parameter. The second polymer is obtained as the complementary of the exclusion polymer. At this point, the spatial distribution of aggregates still has to be simulated. Aggregates are modelled as the intersection of the polymer P i and of a Boolean model of spheres with radius R a , independent from the model used to simulate the exclusion polymer (Figure 3b). Again, the model is described by two parameters: the intensity θ i of the underlying Poisson point process, and the radius R i of the implanted spheres. The volume fraction v a of the aggregates is given by

v a = (1 -v e )   1 -e -θi 4π 3
R 3 i   . (6)
The second scale of the material is made of the carbon black particles themselves. These particles are simulated using a Boolean model of spheres of constant radius R cb implanted on the realization of a Cox point process driven by the aggregates. Since the carbon black particles exhibit a turbostratic structure, these particles can interpenetrate. Therefore, in the model, the particles can partially overlap by introduction of a hard-core distance equal to the radius of the particle (Figure 3c).

The intensity of this Cox process is the random measure defined by

θ cb (x) = θ if x ∈ A aggr , 0 otherwise. ( 7
)
The overall construction of the model is summarized in Figure 4. A vectorial simulation of a domain with a 1µm 3 volume is generated in 2 minutes. Concerning the parameters of the model, the radius of the spheres can be measured experimentally. Therefore, the modelling of the second scale is entirely controlled by the single parameter θ cb . Overall, if the volume fraction v e of exclusion polymer and if the radius R cb of the carbon black particles are known, the global model is controlled by the following set of parameters: (v) the hard-core distance for the Cox point process, which is equal to the radius of the particle.

Finally, note that a bound rubber of constant thickness t surrounds the carbon black particles. In the model, the bound rubber is obtained by considering the Boolean difference of a boolean model of sphere of constant radius R cb + t implanted on the centers of the black carbon particles and of the Boolean spheres describing the black carbon particles. This version of the model is close to the version given in [START_REF] Jean | A multiscale microstructure model of carbon black distribution in rubber[END_REF]. However in the previous work, aggregates were obtained from a Boolean model of spheroids, instead of spheres, and no hard-core process was introduced for the location of carbon black particles. Furthermore, the identification process of the parameters of the model was completely different.

Model fitting from experimental data

The model can be identified and validated with the available experimental data. The four free parameters R e , R a , θ a and R cb are estimated by fitting the covariance and the granulometry curves of the simulated images to the corresponding experimental curves.

v e R e (nm) θ a (nm -3 ) R a (nm) θ cb (nm -3 ) R cb (nm) 0.2 30
4.0 × 10 -7 60 3.0 × 10 -5 18

Table 1: Optimized parameters for the Boolean model of exclusion polymer

Optimization algorithm The procedure to fit the model parameters is as follows. First, we start from the experimental covariance C(h) for all distances 0 ≤ h ≤ N and the experimental granulometry curve G(s) for all sizes 1 ≤ s ≤ K. Then, for each realization of the model the least square function is computed

F(R e , θ a , R a ) = N n=0 (C[R e , θ a , R a , n] -C[n]) 2 + K k=1 (G[R e , θ a , R a , k] -G[k]) 2 , (8)
where C (resp. G) is the covariance function (resp. granulometry curve) measured at all distances {0 ≤ n ≤ N } (resp. for all sizes {0 ≤ k ≤ K}) on the realization. The optimization problem is thus the determination of the parameters that minimize the functional F.

Considering the lack of analytical formulae for the covariance in the model, we have to rely on numerical procedures to perform the optimization. Most numerical optimization algorithms rely on steepest-descent or on Gauss procedures. A popular method in numerical optimization is thus the Levenberg-Marcquardt algorithm [START_REF] Levenberg | A method for the solution of certain problems in least squares[END_REF], which combines these two approaches. This algorithm has been successfully used in [START_REF] Jean | A multiscale microstructure model of carbon black distribution in rubber[END_REF] to optimize a similar model for the microstructure of carbon black particles embedded in rubbers. The main issue with the Levenberg-Marcquardt algorithm is that it relies on a gradient calculation. Estimating the gradient of the covariance with respect to the simulation parameters is a difficult task here, since a small variability cannot be avoided due to the intrinsic probabilistic nature of the model. As a consequence, in the present work another classical optimization algorithm is implemented to perform the optimization, namely the Nelder-Mead procedure [START_REF] Nelder | A simplex method for function minimization[END_REF]. The advantage of this algorithm is that it does not rely on the gradients of the functional with respect to the parameters.

Results and discussion

The values of the parameters obtained after optimization are given in the Table 1. Note that a particular emphasis was given to the covariance fitting, as evidenced in Figure 5. A simulated image is shown in Figure 6 and compared to a segmented TEM image.

It is clear from Figure 6 that the model accurately reproduces the reference covariance curve for the studied material. However, a small but still significant discrepancy can be observed between the reference granulometry and the simulations. This is a direct consequence of the model chosen to describe the microstructure. By using a Boolean model to describe the exclusion polymer, strict limitations are imposed on the aggregates geometry. In addition, the model does not preserve the geometrical symmetry between the exclusion polymer and its complementary. As a consequence, the model fails to properly simulate the smallest aggregates, which explains the discrepancy between the granulometry curves. This discrepancy can also result from the aggregate prints removed during the sample preparation: they were considered to be originally mostly constituted of black carbon particles, but this remains an approximation. Hence, the granulometry of the larger grains is largely overestimated, which strongly impacts the granulometry features.

It is finally interesting to take a closer look at the optimal parameters, and in particular at the intensity of the Cox model for the carbon black particles. One can note that the intensity of the Cox model yields a significantly high value. This is explained by the high density of carbon black particles in the aggregates.

Exclusion polymer based on a Johnson-Mehl mosaic model

As indicated previously, a major issue regarding the Boolean model is that it does not preserve the geometrical symmetry between the exclusion polymer and its complementary. It is therefore of interest to consider other models to describe the exclusion polymer. For this purpose, a model for the exclusion polymer based upon Johnson-Mehl tesselations is introduced in this section.

Theory

Let Ω denote a given volume in R 3 . A tessellation in Ω is a subdivision of Ω into three dimensional subsets referred to as cells. A general theory of random tessellation has been established over the years, notably by [START_REF] Matheron | Introduction aux ensembles aléatoires[END_REF]), [START_REF] Stoyan | Stochastic geometry and its applications[END_REF], [START_REF] Zähle | Random cell complexes and generalised sets[END_REF] and [START_REF] Møller | Random tessellations in R d[END_REF]. Another useful reference on this topic is the book of [START_REF] Schneider | Stochastic and integral geometry[END_REF]. Finally, note that the paper of [START_REF] Møller | Random Johnson-Mehl tessellations[END_REF] provides a unified exposition of random Johnson-Mehl tessellations

The Johnson-Mehl tessellation model was introduced by Johnson and Mehl, mainly to describe crystallization processes in metallography. Their model can be seen as a variation of the Voronoï model. A Voronoï tessellation is a tessellation built from a Poisson point process P in the space R 3 . Every point x of R 3 is associated to the class C i containing all points of R 3 closer from the point x i of P than from any other point of P. Johnson-Mehl tesselations can be seen as a sequential version of the Voronoï model, where the Poisson points are implanted sequentially with time. All classes grow then isotropically with the same rate, and the growth of crystal boundaries is stopped when they meet. All Poisson points falling in an existing crystal are removed.

From a mathematical viewpoint, a Johnson-Mehl tessellation is constructed from a sequential Poisson point process where the points x i , i = 1, .., N are implanted sequentially at a time t i , i = 1, .., N . The classes C i , i = 1, .., N corresponding to the points x i , i = 1, .., N are defined by

C i = y ∈ R 3 , ∀j = i, t i + x i -y v ≤ t j + x j -y v . (9)

Exclusion model

In our multiscale model, we rely now on a Johnson-Mehl tesselation to describe the exclusion polymer and its complementary: each class C i is attributed to the exclusion polymer with the probability v e (= 0.2 in the present case) and to the other polymer with the probability 1 -v e , generating a binary Johnson-Mehl mosaic (Figure 7a), that allows us to preserve the symmetry between the geometry of the exclusion polymer and its complementary. (iv) the intensity θ cb of the Cox point process on which relies the modelling of the spatial distribution of the black carbon particles.

Finally, the bound rubber is obtained as before by considering the Boolean difference of a Boolean model of sphere of constant radius R cb + t implanted on the centers of the black carbon particles and of the Boolean spheres describing the black carbon particles.

Model fitting from experimental data

As previously, the parameters of the Johnson-Mehl exclusion model are estimated to match the morphological descriptors obtained on the experimental TEM images. The values of the parameters obtained after optimization are given in the Table 2. The experimental covariance and granulometry are compared with the corresponding curves obtained by the simulations with fitted parameters on Figure 9. It turns out that the agreement between experimental data and measurements made on simulations is better than in the previous version of the model. A simulated image is shown in Figure 10 and compared to a TEM image.

Implementation of the models

Simulations of random structures are generally performed on a grid of points (i.e. 2D or 3D images), using primary grains based on combination of pixels. The implementation in the software VtkSim relies on a completely different approach based upon level sets and implicit functions. In this approach, primary grains are described by implicit functions [START_REF] Faessel | 3D multiscale vectorial simulations of random models[END_REF][START_REF] Faessel | VtkSim software[END_REF], which are real valued functions defined in the 3D space. The level sets of an implicit function Φ are described by an equation of the form Φ(x, y, z) = c, for some constant c. A surface is described as a level set of the function Φ, most commonly the set of points for which Φ(x, y, z) = 0. In this case, the points for which Φ(x, y, z) < 0 correspond to the interior of the primary grain associated to the implicit function, the points for which Φ(x, y, z) > 0 to its complementary and the level set Φ(x, y, z) = 0 to the boundary of the primary grain. Any primary grain can be implemented, whatever its shape, as long as it can be represented using an implicit function.

In the implicit function approach, complete simulations are generated, using Boolean combinations of primary implicit functions: the union and the intersection of two objects A 1 and A 2 are defined to yield the minimum and the maximum, respectively, of their corresponding implicit functions. Thus, we have

Φ(A 1 ∪ A 2 ) = min{Φ(A 1 ), Φ(A 2 )} (10) and Φ(A 1 ∩ A 2 ) = max{Φ(A 1 ), Φ(A 2 )}. (11)
Similarly, the complementary A c of set A is defined to be the opposite function

Φ(A c ) = -Φ(A). (12)
Several classical random models can be generated using this method, including among others soft-core or hard-core Boolean models of spheres. Interestingly, since a combination of implicit functions is itself an implicit function, the realization of a first random simulation enables us to generate multiscale simulations. An interesting feature of the software is that it allows to directly export the simulated microstructures as a mesh for further use, like finite elements based computations. Overall, using implicit functions to perform the simulation allows us to build complex combinations of simulations that we could not process with a pixel based method. Furthermore, vectorial simulations do not require a large amount of computer resources, which allows us to obtain fast generations of microstructures with very large sizes and numbers of primary grains.

Percolation of the carbon black network

It is known that there is a large impact of the connectivity of elements in a random system on its macroscopic behavior. Controlling the percolation of a network is of major interest in the optimization of the microstructure. In this section, we focus on the percolation of the clusters constituted by the filler particles. Percolating networks of filler can indeed significantly change the mechanical properties of the material. It is therefore of interest to be able to select simulation parameters that will (or not) create percolating paths in the material.

General approach

A simulation is said to percolate when a continuous path can be drawn in the filler phase that connects two opposite faces of the simulation window (a cube in the 3D space). Our aim is to determine the probability p that a simulation in a cubic window of length L percolates. To estimate p, N simulations of the random model are performed, and the number of simulations that percolate is recorded. A set of parameters leads to a percolating medium when more than 50 percents of the simulations percolate. In this approach, p is estimated by

pN = 1 N N n=1 P k , (13)
where P k is the random variable equal to 1 when the simulation percolates and 0 otherwise. The estimator pN is without bias. We denote by σ 2 = p(1 -p) the variance of the random variables {P i , i = 0, .., N }. The variance of pN is then given by

Var(p N ) = σ 2 N . (14)
Hence, Var(p N ) → 0 when N → +∞ and the estimator p converges. For N large enough, according to the Central Limit Theorem, we have pN

∼ N (p, σ 2 /N ), (15)
where N (p, σ 2 /N ) denotes the normally distributed law with mean p and standard deviation σ/ √ N . Therefore, it can be show that, for N = 100, with confidence 95%,

p -1.96 p(1 -p) N ≤ p ≤ p + 1.96 p(1 -p) N . (16)
Since the quantity p(1 -p) is maximal for p = 1/2, yielding p(1 -p) = 1/4, we have, in the worst case,

p -0.098 ≤ p ≤ p + 0.098. ( 17
)
In mathematical terms, we have P p -0.098 ≤ p ≤ p + 0.098 ≥ 0.95. [START_REF] Jeulin | Morphological modeling of random composites[END_REF] This result provides an estimation of the precision of the results obtained in the percolation study.

Percolation algorithm

The aim of the percolation algorithm is to label the connected clusters of filler particles. For each connected cluster, it is then determined if it is percolating. The percolation algorithm proceeds as follows:

(i) A realization of the model is simulated in a volume of fixed size. At the end of the simulation, the algorithm returns a list of the coordinates of the filler particles centers.

(ii) Following the approach of [START_REF] Moreaud | Multi-scale simulation of random spheres aggregates: Application to nanocomposites[END_REF], the volume is subdivided into cubes of size R, R being the radius of the filler particle. The coordinates of the filler particles centers are then sorted according to the subvolumes.

(iii) The classical two-pass algorithm is used to label all clusters of filler particles. For a given particle, the number of test of intersection is restricted to the adjacent cubes, thus yielding a fast method of construction.

(iv) For each labelled cluster, the algorithm checks if the cluster is percolating.

Results

The effect of various morphological parameters on the proportion of percolating simulations, estimated on 100 realizations for each set of parameters, was investigated for the two previous models of multiscale random microstructures: volume fraction of filler, size of aggregates, size of the exclusion polymer clusters. Two sizes of simulations were used, namely (400 nm) 3 and (800 nm) 3 . We report here results obtained on the largest simulations.

Volume fraction of fillers

The volume fraction of fillers is controlled in the simulations by means of the intensity θ cb for the centers of the carbon black spheres. The remaining parameters of the simulations are as follows: v e = 0.2; R e = 30 nm (Boolean model for the exclusion polymer), θ e = 4.0 × 10 -7 nm -3 (Johson-Mehl mosaic); v a = 0.47; R a = 60 nm; R cb = 18 nm; simulated volume: (800 nm) 3 . The effect of θ cb is illustrated by Figures 11a and11b. From these results, it turns out that the percolation threshold p c (i.e. the lowest value of p for which more than 50% of realizations percolate) is given by v cb = 0.15 for the Boolean model of exclusion and by v cb = 0.16 for the Johnson-Mehl mosaic. Therefore, the second model shows a slightly higher percolation threshold. Since the nominal carbon black volume fraction of the real material is v cb = 0.14, it is expected that its filler network does not percolate, which prevent the material to show a rigid macroscopic behavior.

Aggregates size The effect of the size R aggr of the filler aggregates of the material was studied, while keeping constant the volume fraction of aggregates. We remind that the volume fraction of aggregates is given by

v a = (1 -v e )   1 -e -θi 4π 3 R 3 i   . ( 19
)
By properly varying the intensity θ aggr , the volume fraction of aggregates remains constant. The simulation parameters are v e = 0.2, R e = 30 nm (Boolean model), , θ e = 4.0 × 10 -7 nm -3 (Johnson-Mehl mosaic); v a = 0.47; θ cb = 3.5 × 10 -5 nm -3 , v cb = 0.14 and R cb = 18 nm.

The results are displayed in Figures 11c and11d. Interestingly, we note that, all things otherwise equal, increasing the size of the aggregates results in a monotonous increase of the percentages of simulations that percolate. For instance the percentage of percolating simulations climbs from 40 percents when R aggr = 55 nm to 60 percents when R aggr = 70 nm for the two types of models, so that percolation in enhanced for larger aggregates.

Exclusion polymer clusters size From simulations generated by variations of the intensity θ e of the Poisson point process, the number of exclusion polymer clusters increases (and thus their size is reduced since the volume fraction of exclusion polymer is kept constant). This results in a very slight decrease of the number of percolating simulations for both models. Overall, the effect of the exclusion polymer clusters size on percolation properties remains rather small.

Prediction of viscoelastic properties of simulated 3D microstructures by Fast Fourier Transform (FFT)

In this section a preliminary study of the viscoelastic behavior of the models of random microstructure is presented. For simplicity, the approach is restricted to the case of two-components media, namely a viscoelastic matrix and an elastic filler. In the present, accordingly, the mechanical behavior of the various polymers (exclusion polymer, polymer containing the filler, thin layer around carbon-black particles) is supposed to be the same. In addition, only small deformations of the polymer are considered. Using a Fourier decomposition in time, strain and stresses are specified as a series of harmonics. This is commonly accounted for by the introduction of complex elastic moduli. A Fourier decomposition in space is introduced, making use of FFT methods to compute the elastic response of highly-contrasted composites. This method is first validated using analytical results for periodic media and compared with explicit timediscretization. It is efficiently applied to arbitrarily complex microstructures provided by a regular grid of voxels, such as images of real or of simulated microstructures, and does not require any further meshing of images. Extended to the viscoelastic case, it provides full-field reconstructions in time-space, from which the effective viscoelastic properties of a composite medium are estimated.

In what follows, a complex FFT scheme for viscoelastic composites is introduced and applied to 3D simulations of the two previously introduced versions of the multiscale model.This section presents a complex Fourier (FFT) scheme for computing the viscoelastic response of composites. Linear strain-stress relations with elastic and viscoelastic phases are considered. The problem to solve is presented in Section (5.1). The FFT method used to solve the problem is described in Section (5.2). It is validated in Section (5.3) and used to predict the properties of rubber containing fillers in Section (5.4).

Problem setup

Consider a domain V representing a heterogeneous material subject to steady-state oscillatory conditions. In the harmonic regime, the time-dependence of the strain ε and stress fields σ at point x is given by the complex fields:

ε(x; t) = ε(x) e iωt , σ(x; t) = σ(x) e iωt , (20)
where the amplitudes ε and σ are second-rank tensor fields, and ω is the frequency applied to the system. In the above, ε(x) and σ(x) are complex and do not depend on the time t. The physical strain and stress fields are the real parts of ε and σ:

ε(x; t) = R [ ε(x; t)] , σ(x; t) = R [ σ(x; t)] . (21)
For simplicity, it is assumed that the material contains elastic inclusions (phase α = 2) embedded in a viscoelastic matrix (phase α = 1). Referring to C 2 as the inclusions's isotropic elastic stiffness tensor, we have, in phase 2:

σ(x; t) = C 2 : ε(x; t), σ(x) = C 2 : ε(x), (22)
or equivalently,

σ kk (x) = 3κ 2 ε kk (x), σ (x) = 2µ 2 ε (x). ( 23
)
where µ 2 is the inclusion's shear modulus, κ 2 its bulk modulus, and where the strain and stress deviatoric parts are defined as:

ε ij = ε ij -(ε kk /3)δ ij , σ ij = σ ij -(σ kk /3)δ ij , (24)
δ being the Kronecker symbol. Inclusions are embedded in a viscoelastic matrix, where the local strain-stress relation is written generally as:

σ(x; t) = t -∞ dτ C 1 (t -τ ) : d ε(x; τ ) dτ . (25)
In the above, C 1 (t) designates a time-dependant fourth-order isotropic tensor with shear and bulk (strain-rate) moduli µ 1 (t) and κ 1 (t). Following [START_REF] Christensen | Theory of viscoelasticity: an introduction[END_REF] the time-dependence is removed using Eq. ( 20):

σ kk (x) = 3κ * 1 (iω)ε kk (x), σ (x) = 2µ * 1 (iω)ε (x), (26)
for x in the matrix. The complex elastic moduli κ * 1 (iω) and µ * 1 (iω) are:

κ * 1 (iω) = κ 1 + iω ∞ 0 dη [κ 1 (η) -κ 1 ] e -iωη , (27a)
µ * 1 (iω) = µ 1 + iω ∞ 0 dη [µ 1 (η) -µ 1 ] e -iωη , (27b)
with:

κ = lim t→∞ κ 1 (t) ≥ 0, µ = lim t→∞ µ 1 (t) ≥ 0. ( 28
)
The local response of the material is given by Eqs. ( 23) and [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. Viscous effects are taken into account by the imaginary parts of the elastic moduli µ 1 (iω) and κ 1 (iω). When the latter are positive, the stress fields's oscillations are advanced by an amount of:

∆t = arg [σ(x)] -arg [ε(x)] ω , (29)
compared to that of the strain field. For instance, if ε(x) is real:

σ(x; t) = R [σ(x)] e iω(t-∆t) . (30)
Such effects are absent in the purely elastic case [START_REF] Matheron | Introduction aux ensembles aléatoires[END_REF].The material equations above are completed by admissibility equations for the strain and stress fields. Assuming first that the displacement field u(x; t) is infinitesimal compared to the microstructure characteristic sizes, we have:

ε ij (x; t) = ∂ j u i (x; t) + ∂ i u j (x; t) 2 . (31)
Second, stress equilibrium is met at every point:

∂ i σ ij (x; t) = 0. (32)
We define the complex displacement amplitude u(x) by:

u(x; t) = u(x) e iωt (33)
so that Eqs. ( 31) and (32) become:

ε ij (x) = ∂ j u i (x) + ∂ i u j (x) 2 , ∂ i σ ij (x) = 0. (34)
Applying periodic boundary conditions with time-harmonic macroscopic strain loading as:

ε(x; t) = ε e iωt , ε(x; t) #, σ(x; t) #, (35)
where ε is a constant, • designates a spatial average over the domain V and # denotes a periodic field in space with elementary cell V . The above is indifferently written as:

ε(x) = ε, ε(x) #, σ(x) #. ( 36
)
Hereafter it is assumed, without loss of generality, that all components of ε are real. The composite's effective response is given by the stiffness tensor C eff :

σ = σ(x) = C eff : ε. ( 37
)
When the composite is isotropic and V is large enough, the tensor C eff is isotropic as well. Its shear and bulk moduli are noted µ eff and κ eff . Eqs. ( 23), ( 26), [START_REF] Savary | Morphological analysis of carbon-polymer composite materials from thick sections[END_REF] and [START_REF] Serra | Image analysis and mathematical morphology[END_REF] define the response of the material to a harmonic solicitation. In theory, these solutions allow one to reconstruct the local and macroscopic responses of the material subjected to non-harmonic (arbitrary) strain loading, in the steady-state regime [START_REF] Christensen | Theory of viscoelasticity: an introduction[END_REF]. The first equation in [START_REF] Schneider | Stochastic and integral geometry[END_REF] can be replaced by:

α(t) = ε(x; t) , (38)
so that a time-varying macroscopic strain loading α(t) is applied. Taking the Fourier transform of the latter:

α(t) = 1 2π ∞ -∞ dω α(ω) e iωt . ( 39
)
The strain field acting in the material is the superposition of the following harmonic responses:

ε(x; t) = 1 2π ∞ -∞ dω ε ω (x) e iωt , (40)
where ε ω (x) is the complex field solution of Eqs. ( 23), ( 26), ( 34) and ( 36) with ε = α(ω).

FFT methods

The purpose of this section is to solve numerically Eqs. ( 23), ( 26), ( 34) and ( 36) on a grid of voxels, using "Fourier methods" as developed by [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF][START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. The "direct scheme" [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF] and "augmented Lagrangian" [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] Fourier methods apply to linear elasticity. The "accelerated scheme" [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], originally given in the context of conductivity, has been extended to linear elasticity by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] . We remark that the problem ( 23), ( 26), [START_REF] Savary | Morphological analysis of carbon-polymer composite materials from thick sections[END_REF] and ( 36) is mathematically identical to that of linear elasticity, except that complex instead of real fields are concerned. Thus, FFT methods apply to the viscoelastic problem considered here, as will be seen later. This is true of other more recent FFT methods like the "rotated" [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] and "variational FFT" [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] schemes. Accordingly, any of these schemes can be implemented for this purpose.

Hereafter, the "discrete scheme" proposed by [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise-disordered porous media[END_REF], coupled with the algorithm of [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF], is used. A similar scheme has already been investigated in the context of conductivity by [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF] leading to much improved convergence properties, especially for highly-contrasted composites. The scheme gives for the k th -iteration:

ε k+1 (q) = ε -G 0 (q) : σ(q) -C 0 : ε k (q) , (41)
where q is the Fourier wave vector, C 0 is a reference elastic tensor and G 0 (q) its associated Green operator, discretized using a forward and backward finite difference scheme on the rotated grid. The stress field is computed at each step in the direct space using the material's constitutive law, and the strain field is initialized by ε k=0 ≡ ε. Convergence strongly depends on C 0 , an arbitrary (complex) reference tensor, which has to be optimized. An absolute convergence criterion is defined on stress equilibrium as:

max x |divσ(x)| ≤ η ( 42
)
where η is the required precision. Care is needed when choosing C 0 . To ensure that the physical strain and stress fields are symmetric, i.e.:

ε kl (x; t) = ε lk (x; t), σ kl (x; t) = ε lk (x; t), (43)
it is necessary to enforce:

ε kl (x) = ε lk (x), σ kl (x) = σ lk (x). (44)
Hence, the tensors ε and σ are not symmetric in the complex sense. To enforce this, the Green tensor must satisfy:

G 0 ij,kl (q) = G 0 ji,kl (q) = G 0 ij,lk (q) = G 0 kl,ij (q) * . (45)
In turn, this implies:

C 0 ij,kl = C 0 ji,kl = C 0 ij,lk = C 0 kl,ij * = C 0 kl,ij , (46)
so that C 0 must be real. This is in contrast with optics where we are free to choose complex references (see [START_REF] Azzimonti | Optical properties of deposit models for paints:full-fields FFT computations and representative volume element[END_REF]).

In the rest of this work, we choose κ 0 = 0.51 [Re(κ Figure 13: Macroscopic strain ε 12 = ε 12 (x; t) and resulting stress component σ 12 = σ 12 (x; t) vs. time t, for the time-discretized scheme. The macroscopic strain is 0 for t < π/(2ω). In violet: comparison with the harmonic regime as predicted by the complex Fourier scheme.

20%. For the ratio Im(µ p,s )/Re(µ p,s ), good agreement is obtained for quite smaller volume fractions (not larger than 10%). Finally, the Fourier method developed in Sec. (5.2) is compared to an explicit time-discretized scheme. Assume that: ε(x; t) = 0, σ(x; t) = 0 (52) for t < t 0 . For t ≥ t 0 , the material is subject to harmonic strain loading:

ε(x; t) = cos(ωt)ε (53)
Hereafter we take t 0 = π/(2ω) and again ω = 1 [rad/s] for simplicity. The stress field σ at time t + dt is then computed by discretizing Eq. (47). For instance, for the deviatoric parts:

∆ σ (x, t) = 2µ 0 ∆ ε (x, t) - σ (x, t)∆t t 1 , (54)
∆ σ (x, t) = σ (x, t + ∆t) -σ (x, t), ∆ ε (x, t) = ε (x, t + ∆t) -ε (x, t).
The above is mathematicaly equivalent to a thermoelastic stress-strain relation with unknown ∆ ε and ∆ σ and with applied strain loading given by (53):

∆ ε (x, t) = -ω sin(ωt)∆tε (55)
In the inclusions, the fields ∆ σ(x; t) and ∆ ε(x; t) are related by Hooke's law. The FFT method for thermoelasticity proposed by [START_REF] Vinogradov | An accelerated fft algorithm for thermoelastic and non-linear composites[END_REF], [START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF] is then applied at each time step. We set ∆t = 0.01 and apply shear loading with ε 12 = ε 21 = 1 and ε ij = 0 for (i, j) = (1, 2), [START_REF] Azzimonti | Optical properties of deposit models for paints:full-fields FFT computations and representative volume element[END_REF][START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF]. A comparison between the FFT scheme with time-discretized strain-stress relation and that with complex elastic moduli, is represented in Figure 13 for the macroscopic response. The macroscopic stress σ 12 = σ 12 is given as a function of time t. It is emphasized that, for the time-discretized sheme, σ 12 (t) (top dashed line) is computed at many loading step whereas, for the FFT method involving complex moduli one computation is carried out which gives the harmonic response at t → ∞ (bottom dashed line). As expected, beyond an intermediary regime, the time-explicit scheme converges to a harmonic response identical to that predicted by the complex FFT scheme.

Computations on random two-scales media

To illustrate the efficiency of the method, viscoelastic computations are performed on the Boolean and Johnson-Mehl two-scales models. The local properties are set to 78 400 and 30 000 MPa for the bulk and shear moduli, resp., in the filler. In the polymer, a viscoelastic law with Prony series given by [START_REF] Laiarinandrasana | Modelling the effects of various contents of fillers on the relaxation rate of elastomers[END_REF] is used. Complex elastic moduli are used for 23 frequencies. Each 3D computation required about 3 hours on a 12-cores computer and 800 iterations with η = 3.10 -5 . No time discretization is used. Domains with size (0.8 µm) 3 discretized over (400) 3 voxels are simulated. FFT maps of the mean stress field (opposite of pressure) are shown in Figures 14a and14b, showing a 2D section for the exclusion zones modelled by a Boolean model. The material is submitted to a hydrostatic strain loading, the frequency being ω = 117 Hz. Strong internal stresses in the polymer appear in regions surrounded by fillers. Maps of the shear stress component with shear strain loading are shown in Figures 14c and14d.

The effective moduli are obtained by averaging the fields calculated by FFT, using a standard homogenization approach. The relative precision of the estimated moduli is calculated from the variance of the fields at different scales, as developed in the statistical RVE approach of [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. At a frequency ω = 117 Hz, the relative precision for the real part and the imaginary part of the bulk modulus is equal to 0.2% and to 0.3%. For the the real part and the imaginary parts of the shear modulus, it is equal to 17% and 3.8% respectively.

The frequency dependence of the effective complex bulk modulus and of the effective shear modulus are illustrated in Figure 15. FFT data for the Boolean and Johnson-Mehl models are compared to the viscoelastic response of the Hashin sphere assemblage [START_REF] Hashin | Complex moduli of viscoelastic composites-I. general theory and application to particulate composites[END_REF]. The inner spheres and outer layers in the Hashin coating have the same properties as the fillers and polymer respectively, and the volume fraction are the same as in the model (about 10.6%). It has been shown that the bulk and shear moduli of isotropic two-phases viscoelastic composites are constrained to a region in the complex plane [START_REF] Gibiansky | On the effective viscoelastic moduli of two-phase media. i. rigorous bounds on the complex bulk modulus[END_REF][START_REF] Milton | On the effective viscoelastic moduli of two-phase media. ii. rigorous bounds on the complex shear modulus in three dimensions[END_REF] that generalizes the Hashin and Shtrikman bounds. The Hashin coating has extremal properties for the bulk modulus in the sense that its bulk modulus lies on the frontier of the region. For the material considered here, the contrast of properties in terms of bulk moduli is moderate (about 25), and so the Hashin coating is very close to that of the Boolean and Johnson-Mehl models, as computed using FFT (Figures 15a and15b). The situation is sensibly different for the real part of the shear modulus, which is 3.10 4 higher in the filler than in the polymer (Figure 15c). The real part of the shear modulus of the Johnson-Mehl model is higher than that of the Boolean model.

We note that the viscoelastic response of the models is much softer than that of the Hashin assemblage with stiff outer coating, which is not represented in figure 15. Our results are also consistent with that obtained by [START_REF] Jean | Large scale computations of effective elastic properties of rubber with carbon black fillers[END_REF] in the context of elasticity, where elastic moduli close to the Hashin and Shtrikman lower bounds were predicted.

Conclusion

Multiscale morphological models of random 3D microstructures were developed to simulate the distribution of fillers in rubber materials, accounting for the presence of aggregates and of an exclusion polymer. These models enable us to study the percolation of carbon black filler as a function of its multiscale distribution. It was proved that the two models involving a different morphology for the exclusion polymer do not percolate for the volume fraction of carbon black contained in the material.

The effective viscoelastic moduli are predicted by Fast Fourier Transforms on 3D simulated images of the microstructure. The simple method proposed in this work predicts the harmonic response of composites with arbitrary microstructures, frequency per frequency and is restricted to linear materials. The effective bulk modulus and shear modulus are predicted in a large frequency range (about from 10 -1 to 10 4 rad/s).

Further studies will involve a systematic change of the microstructures by means of the parameters of the model, and its impact on the effective properties, to find optimal microstructures with respect to the macroscopic behavior.

Figure 1 :Figure 2 :

 12 Figure 1: (a) TEM Image of rubber matrices filled with carbon black, obtained in the Yokohama Rubber Research Center. Different steps of the segmentation (see text): (b) Filtered TEM Image; (c) Extracted white prints; (d) Dual Top-hat; (e) extracted aggregates; (f) segmented image.

Figure 3 :

 3 Figure 3: The three steps of the construction of the model: (a) exclusion polymer (Boolean model of spheres); (b) aggregates; (c) carbon black particles network (v e = 0.2; R e = 30 nm; R a = 60 nm; R cb = 20 nm; simulated volume: (500 nm) 3 ).

Figure 4 :

 4 Figure 4: Scheme of the multiscale random model (exclusion zones made of a Boolean model of spheres).

Figure 5 :

 5 Figure 5: Covariance (a) and granulometry (b) of the simulated material. The simulated curves are represented by the discontinuous lines, the experimental curves by the continuous lines.

Figure 6 :Figure 7 :

 67 Figure 6: Comparison between an experimental (a) and a simulated TEM image (b) for the Boolean model of exclusion polymer.

  Here, the two previous parameters of exclusion zone model (R e and θ e ) are replaced by a single one, namely the intensity θ e of the Poisson point process associated to the Johnson-Mehl tesselation. The remaining parts of the model are the same as previously. Simulation parameters are the intensity of the Poisson point processes attached to each polymer. A realization of the Johnson-Mehl based model is displayed in Figure 7, as an illustrative example. A vectorial simulation of a domain with a 1µm 3 volume is generated in 10 minutes. The overall construction of the model is summarized in Figure 8. If the volume fraction v e of exclusion polymer and the radius R cb of the carbon black particles are known, the global model is controlled by the following set of parameters: (i) the intensity θ e of the Poisson point process associated to the Johnson-Mehl tesselation (ii) the intensity θ a of the Poisson point process on which relies the modelling of the aggregates, (iii) the radius R a of the boolean spheres used to model the aggregates,

Figure 8 :

 8 Figure 8: Scheme of the multiscale random model (exclusion polymer made of a Johnson-Mehl mosaic).

18 Table 2 :

 182 v e θ e (nm =3 ) θ a (nm -3 ) R a (nm) θ cb (nm -3 ) R cb (nm) 0.2 4.0 × 10 -7 4.0 × 10 -7 60 3.5 × 10 -5 Optimized parameters for the Johnson-Mehl model of exclusion polymer.

Figure 9 :

 9 Figure 9: Covariance (a) and granulometry (b) of the simulated material (Johnson-Mehl mosaic for the exclusion polymer). The simulated curves are represented by the discontinuous lines, the experimental curves by the continuous lines.

Figure 10 :

 10 Figure 10: Comparison between an experimental (a) and a simulated TEM image (b) for the Johnson-Mehl mosaic model of exclusion polymer.

Figure 11 :

 11 Figure 11: (a,b): Influence of the density θ cb of filler (with V e and V a fixed) on the proportion of percolating simulations for the two models of exclusion zones: (a) Boolean model of spheres; (b) Johnson-Mehl mosaic. Simulated volume: (800 nm) 3 . (c,d): Influence of the aggregates size (with V e , V a and v cb fixed) on the proportion of percolating simulations for the two models of exclusion zones: (c) Boolean model of spheres; (d) Johnson-Mehl mosaic. Simulated volume: (800 nm) 3 .

Figure 14 :

 14 Figure 14: (a-b): Real and imaginary parts of the mean stress field σ m for a hydrostatic strain loading (ε m = 1%) at a frequency ω = 117 Hz. The minimal and maximal values are -47 and 140 Gpa (real part (a)) and -1.0 and 6.2 GPa (imaginary part (b)). The stress values are thresholded between 8.68 and 8.73 GPa (real part) and 1.5 and 1.6 (imaginary part) to highlight the field patterns in the polymer. (c-d): Real and imaginary parts of the shear stress field σ 12 for a shear strain loading (ε xy = 1%) at a frequency ω = 117 Hz. The minimal and maximal values are -17 and 44 GPa (real part (c)) and -1.2 and 0.7 GPa (imaginary part (d)). The stress values are thresholded between 0 and 6 MPa (real part) and 0 and 1 (imaginary part).

Figure 15 :

 15 Figure 15: Comparison between the Johnson-Mehl and Boolean sphere models for the exclusion polymer and Hashin sphere assemblage: frequency dependence of the effective real part (a) and imaginary part (b) of the complex bulk modulus. (c): Real part of the complex shear modulus. (d): Imaginary part.

  Comparison between analytical estimates given by[START_REF] Cohen | Simple algebraic approximations for the effective elastic moduli of cubic arrays of spheres[END_REF] (solid and dashed lines) and FFT results (symbols) for the effective elastic response of a periodic array of spheres with varying volume fraction: imaginary part of the effective shear moduli µ p,s (a) and ratio Im(µ p,s )/Re(µ p,s ).

	2 3 Im(µ s,p )	Im(µ s ) (Cohen) Im(µ p ) (Cohen) Im(µ s ) (FFT) Im(µ p ) (FFT) σ 12 , ε 12 2 0 σ 12 (complex 1 scheme)	0.01 1-Im(µ s,p )/Re(µ s,p ) 0.015 Im(µ s ) (Cohen) Im(µ p ) (Cohen) Im(µ p ) (FFT) Im(µ s ) (FFT)
			0.005
	1	-1
	0 Figure 12: 0 0.1 0.2 0.3 0.4 0.5 spheres volume fraction (a) 5 -2 σ 12 (time-explicit)	0 0.1 0.2 0.3 0.4 0.5 (b) 15 t (s) spheres volume fraction 0 ε 12 (time-10 explicit)
	* 1 ) + Re(κ * 2 )], µ 0 = 0.51 [Re(µ * 1 ) + Re(µ * 2 )], for the bulk and shear moduli of the reference tensor, respectively, where µ * 1,2 , κ * 1,2 are the complex elastic moduli in phases 1
	and 2.

Acknowledgements F. Willot is indebted to S. Forest for helpful suggestions regarding the FFT-based computations.

Validation and comparison with an explicit time-discretized FFT scheme

In this section, a periodic array of spheres with cubic symmetry is considered. The spheres are purely elastic with shear and bulk moduli µ 2 = κ 2 = 1000 GPa, and, for simplicity, the embeding medium is supposed to be Maxwellian. Outside the spheres, the stress σ(t) and strain fields ε(t) satisfy:

with µ 0 = κ 0 = 1 GPa, t 1 = 1 s. This is equivalent to [START_REF] Christensen | Theory of viscoelasticity: an introduction[END_REF]:

with:

In the above, H(t) is the Heaviside step function. Using [START_REF] Meyer | Iterative image transformations for an automatic screening of cervical smears[END_REF] the complex elastic moduli in the embedding medium read, as functions of the angular frequency ω:

Heareafter the frequency is set to ω = 1[rad/s]. A sphere is discretized in a cubic domain of size 64 × 64 × 64 voxels. The FFT method described in Sec. (5.2) is then used to solve the viscoelastic problem in the complex domain. The spheres composite has a macroscopic response with cubic symmetry of bulk moduli κ s and shear moduli µ p , µ s :

in Voigt notation where ε = (ε xx , ε yy , ε zz , 2ε yz , 2ε xz , 2ε xy ) and σ = (σ xx , σ yy , σ zz , σ yz , σ xz , σ xy ). FFT Results are compared to the analytical estimates in [START_REF] Cohen | Simple algebraic approximations for the effective elastic moduli of cubic arrays of spheres[END_REF] (Figure 12). Excellent agreement between the FFT results and analytical estimates is found for the complex shear moduli µ p,s when the volume fraction of spheres is less than