B Pier

Signalling problem in absolutely unstable systems

Keywords: absolute instabilities, causality, signalling problem, WKBJ

The dynamics of unstable systems crucially depends on the nature of the instability, either convective or absolute. The signalling problem, which is the study of the spatial response to a localized time-harmonic forcing, is generally believed to be relevant only for stable or convectively unstable systems and to be ill-posed for absolutely unstable systems, where the self-sustained perturbations grow faster than the forced harmonic response.

The present investigation shows that the signalling problem may still be well-posed for media displaying absolutely unstable regions. Considering weakly spatially inhomogenous systems, conditions are derived for the validity of the signalling problem. The complete spatial response to harmonic forcing is first analytically derived in terms of asymptotic approximations and then confirmed by direct numerical simulations.

Introduction

In spatially homogenous systems, linear stability characteristics are derived from the dispersion relation ω = Ω(k) between the frequencies ω and the wavenumbers k of normal modes of the form e i(kx-ωt) , with x and t denoting streamwise distance and time respectively. These stability properties can be understood by resorting to different methods: temporal, spatial or spatio-temporal [START_REF] Huerre | Hydrodynamic instabilities in open flows[END_REF][START_REF] Huerre | Open shear flow instabilities[END_REF][START_REF] Schmid | Stability and transition in shear flows[END_REF].

In a temporal approach, a spatially harmonic perturbation of real wavenumber k is considered. This wavelike initial perturbation evolves in time with a complex frequency ω. Its spatial structure, determined by the wavenumber k, remains unchanged while its amplitude grows or decays with time. Growth or decay is determined by the sign of the temporal growth rate 1 ω i , while propagation takes place with a phase speed ω r /k. This analysis based on real wavenumbers and complex frequencies is known as the temporal problem.

In a spatial approach, localized harmonic forcing is applied with real frequency ω, say at x = 0. The spatial response to this forcing yields waves with complex wavenumbers k. The wavelength of the spatial response is determined by k r while the spatial growth or decay depends on k i : for x → +∞, the spatial response grows when k i < 0 and decays when k i > 0; the reverse holds for x → -∞. The analysis based on real frequencies and complex wavenumbers is known as the spatial, or signalling, problem.

The full spatio-temporal stability properties may be investigated by applying an impulsive localized perturbation: the analysis of the resulting wavepacket yields the complete dispersion relation between complex wavenumbers k and complex frequencies ω. While the impulsively started wavepacket decays in stable media, a growing response develops from the impulse location in unstable systems. If the growing wavepacket propagates away from its source and eventually leaves the medium unperturbed, the instability is said to be convective. If, by contrast, the instability grows in place and invades the system both upstream and downstream, the instability is said to be absolute. Convectively unstable (cu) systems do not display intrinsic dynamics and essentially behave as amplifiers: external perturbations are amplified while propagating through the system, and without continuous external input the medium returns to its unperturbed state. By contrast, absolutely unstable (au) systems display nontrivial dynamics without external forcing: perturbations expand in both upstream and downstream directions so as to cover the entire domain and continue to grow exponentially at every point.

These stability concepts remain valid locally for spatially inhomogenous systems, provided the characteristic inhomogeneity length scale is large compared to a typical instability length scale. However, the connection between local stability characteristics and the long-term global dynamics of spatially developing systems is far from obvious. In a linear framework, it has been shown [START_REF] Chomaz | Bifurcations to local and global modes in spatially developing flows[END_REF][START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF][START_REF] Dizès | Linear global modes in spatially developing media[END_REF] that the presence of local absolute instability is a necessary but not sufficient condition for global instability: in general an au region of finite extent is required for a spatially developing medium to become globally unstable. Thus there exists a wide range of parameter settings for which a medium does not support any self-sustained fluctuations despite the presence of a region of local absolute instability. In such a situation, the linear signalling problem is legitimate and this is precisely the class of systems addressed in the present paper.

Globally stable but locally absolutely unstable systems are encountered in a variety of configurations of practical interest, among which wakes and boundary layers: the cylinder wake flow for Reynolds numbers in the range 25 < Re < 49 [START_REF] Monkewitz | The absolute and convective nature of instability in two-dimensional wakes at low Reynolds numbers[END_REF], a class of "synthetic" wake flows [START_REF] Pier | Nonlinear self-sustained structures and fronts in spatially developing wake flows[END_REF][START_REF] Pier | Nonlinear synchronization in open flows[END_REF], the wake of a sphere [START_REF] Pier | Local and global instabilities in the wake of a sphere[END_REF], the three-dimensional boundary layer produced by a rotating disk [START_REF] Davies | Global behaviour corresponding to the absolute instability of the rotating-disk boundary layer[END_REF].

The paper is organized as follows: After formulating the problem in terms of the widely used partial differential complex Ginzburg-Landau equation ( §2), its local ( §3) and global ( §4) stability properties are recalled. In §5, the correspondence between the complex space and frequency planes and the structure of the wavenumber branches are analysed. The complete analytic solution to the signalling problem is obtained in §6 in terms of asymptotic approximations and discussed in §7. These results are confirmed by direct numercial simulations in §8.

Problem formulation

Partial differential model equations account for the dynamics of a variety of physical systems [START_REF] Cross | Pattern formation outside of equilibrium[END_REF] and are often tractable by analytical methods. The linearized complex Ginzburg-Landau model (1) has on many occasions proven to be a convenient testground to recognize and study generic features that have later been identified in a variety of situations. The same strategy is adopted here.

The system under consideration is assumed to be described by a complex scalar field ψ(x, t) in an infinite one-dimensional spatially inhomogenous domain and it is governed by

∂ψ ∂t = -i ω 0 (X) + 1 2 ω kk (X)k 0 (X) 2 ψ + ω kk (X)k 0 (X) ∂ψ ∂x + i 2 ω kk (X) ∂ 2 ψ ∂x 2 + S(x, t), (1)
where the complex functions ω 0 (X), k 0 (X) and ω kk (X) account for the local properties of the medium and only depend on a slow space variable X = ǫx. The coefficients of (1) have been cast in this form for reasons that will become clear in the next section. The weak inhomogeneity parameter ǫ ≪ 1 is defined as the ratio of the typical instability length scale to the inhomogeneity length scale of the medium. The source function S(x, t) represents an externally applied forcing to be specified below. While equation (1) applies to the real x-axis, the functions ω 0 (X), k 0 (X) and ω kk (X) are assumed to be analytic and their continuation in the complex X-plane will be used in the following sections.

Local stability properties

In the subsequent discussion, constant use is made of the local properties of system [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF]. Local characteristics are derived from (1) by freezing X to some arbitrary (possibly complex) value and studying the corresponding spatially homogenous system

∂ψ ∂t = -i ω 0 + 1 2 ω kk k 2 0 ψ + ω kk k 0 ∂ψ ∂x + i 2 ω kk ∂ 2 ψ ∂x 2 + S(x, t), (2)
where the dependence of the coefficients on the parameter X has been omitted. Normal modes of the form e i(kx-ωt) , with complex frequency ω and complex wavenumber k, are then governed by the local dispersion relation

ω = Ω(k) ≡ ω 0 + 1 2 ω kk k -k 0 2 . ( 3
)
To satisfy causality, temporal growth rates must be bounded which requires that Im ω kk < 0. The particular form in which the coefficients of ( 1) and ( 2) have been cast brings to the fore the local complex absolute frequency ω 0 associated with the local complex absolute wavenumber k 0 determined by the zero group velocity condition [START_REF] Briggs | Electron-stream interaction with plasmas[END_REF][START_REF] Bers | Space-time evolution of plasma instabilities -absolute and convective[END_REF]

ω 0 = Ω(k 0 ) with ∂Ω ∂k (k 0 ) = 0.
For impulsive forcing of the form S(x, t) = δ(x)δ(t), where δ denotes the Dirac delta function, the absolute frequency ω 0 characterizes the temporal evolution of the resulting wavepacket observed at a fixed station: when Im ω 0 < 0, the system is stable or convectively unstable and the wavepacket either decays or grows while being swept away by advection; when Im ω 0 > 0, the system is absolutely unstable and the impulse response exponentially grows at each point in space. Now consider a spatially localized but temporally harmonic forcing of (possibly complex) frequency ω f , switched on at t = 0, i.e., S(x, t) = δ(x)H(t)e -iω f t where H denotes the Heaviside unit step function. The spatial response to this harmonic forcing consists in traveling waves of frequency ω f and their wavenumbers are derived from (3) as

k ± (ω f ) = k 0 ± 2 ω f -ω 0 ω kk . (4)
Upon choosing the square-root branch cut in (4) along the positive real axis and with the square-root symbol denoting the root with positive imaginary part, the k + and k -branches pertain to the domains x > 0 and x < 0 on either side of the forcing. However, turning on the forcing at t = 0 also produces a transient wavepacket that grows according to the absolute growth rate Im ω 0 at any fixed spatial location. Hence, two different situations arise depending on the relative values of Im ω f and Im ω 0 :

-When Im (ω f -ω 0 ) > 0 the switch-on wavepacket is overwhelmed for large time by the spatial response tuned to the forcing frequency ω f , and the signalling problem is well-posed. -When Im (ω f -ω 0 ) < 0, the switch-on wavepacket overwhelms the spatial response at ω f and the harmonic forcing does not succeed in tuning the medium to the externally applied frequency. The signalling problem is then ill-posed.

Global stability properties

In spatially inhomogenous systems governed by (1), the above discussion yields stability characteristics prevailing locally at each X. However the associated global behaviour cannot be immediately derived and must be analysed carefully [START_REF] Dizès | Linear global modes in spatially developing media[END_REF][START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF].

In a typical situation of interest, the local absolute growth rate Im ω 0 (X) displays a single maximum over the real X-axis and the medium is stable for X → ±∞. In the same spirit as Huerre & Monkewitz [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF] and in order to keep computational difficulties to a minimum, it is assumed that 1 Correspondence by the analytic mapping ω0(X) between (a) complex X-plane and (b) complex ωplane. Associated curves are rendered by the same symbols. The real X-axis is mapped onto the dotted curve {ω0(Xi = 0)}, and the real frequencies ω ca 0 , ω ac 0 correspond to the boundaries X ca , X ac of the au interval. Horizontal lines in (b) are associated in (a) with the corresponding isolines {X1,2(ωi = c st )} of constant absolute growth rate.

ω 0 (X) = ω s + 1 2 ω XXs (X -X s ) 2 , ( 5
) (a) X i X r {X 1 (ω i = c st )} {X 2 (ω i = c st )} X s X ca X ac (b) ω r ω i {ω 0 (X i = 0)} ω s ω ac 0 ω ca 0 Fig.
where ω s , ω XXs and X s are complex parameters with Im ω XXs < 0. By resorting to asymptotic expansions, it has been shown [START_REF] Dizès | Linear global modes in spatially developing media[END_REF][START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF] that the impulse response of the spatially developing system is then dominated for large time by a global mode of frequency ω s . Hence, whenever Im ω s < 0 the system does not display self-sustained solutions: any perturbation eventually decays even though large amplification may be observed in the transient régime. Systems with Im ω s < 0 are thus globally stable whether or not they exhibit an interval of local absolute instability.

The paper is concerned with globally stable but locally absolutely unstable systems. In this situation all transient behaviour displays eventual decay, and it is legitimate to study the long-term spatial response to an externally applied time-harmonic forcing.

Structure of complex X-plane and spatial branches

When media are considered that are both globally stable and locally absolutely unstable, the local absolute frequency is of the form (5) with Im ω s < 0 and positive absolute growth rate (ω 0,i (X) > 0) occurs over some finite interval X ca < X < X ac . Note that this is possible only if the saddle point X s of ω 0 (X) is not located on the real X-axis. Let ω ca 0 ≡ ω 0 (X ca ) and ω ac 0 ≡ ω 0 (X ac ) denote the real marginal absolute frequencies at the boundaries of the au interval.

The function ω 0 (X) defines an analytic mapping between the complex X-plane and the complex ωplane as shown in figure 1. The real X-axis is mapped onto the locus denoted as {ω 0 (X i = 0)} (dotted parabolic curve in figure 1b) crossing the real frequency axis twice at ω ca 0 and ω ac 0 , corresponding to the marginally absolutely unstable positions X ca and X ac . Note that figure 1 illustrates a configuration where ω ca 0 < ω ac 0 and Im X s > 0; situations with ω ac 0 < ω ca 0 and/or Im X s < 0 yield similar pictures and results.

In the complex X-plane (figure 1a), isolines of constant absolute growth rate Im ω 0 (X) are obtained as preimages by ω -1 0 of horizontal lines in the ω-plane. When such a line {ω i = c st } is located above the {ω 0 (X i = 0)} curve in the ω-plane (thick dashed line in figure 1b), the two corresponding isolines in the X-plane (thick dashed curves in figure 1a) do not cross the real axis and may be labeled {X 1,2 (ω i = c st )}, where the subscripts 1 or 2 correspond to curves confined to the upper or lower half X-planes respectively. When the line {ω i = c st } is lowered onto the real ω-axis, the associated X 1 and X 2 curves move towards each other. Since the real ω-axis is crossed by the {ω 0 (X i = 0)}-curve at ω ca 0 and ω ac 0 , one of the associated isolines (indicated by solid curves in figure 1a) necessarily crosses the real X-axis at X ca and X ac . In figure 1(a), it is the lower {X 2 (ω i = 0)}-isoline that crosses the X-axis while the associated {X 1 (ω i = 0)}-isoline lies entirely within the upper half-plane. Note that the (thin-dashed) curves {X 1,2 (ω i = Im ω s )} pinch at the saddle point X s for the line {ω i = Im ω s }, which is however located in the lower frequency half-plane and will not be considered here.

When solving the dispersion relation (3) with coefficients depending on X for a given frequency ω f , the two spatial branches are obtained as

k ± (X, ω f ) = k 0 (X) ± 2 ω f -ω 0 (X) ω kk (X) . (6)
In the complex X-plane, these expressions display two branch points at X 1,2 (ω f ), where ω 0 (X) = ω f , and branch cuts must be introduced. Upon choosing the branch cut for the square-root function in [START_REF] Cross | Pattern formation outside of equilibrium[END_REF] along positive real values of its argument, the causality condition Im ω kk (X) < 0 guarantees that the branch cuts for k ± lie within the regions where Im [ω f -ω 0 (X)] < 0 and thus do not intersect the region located between the {X 1,2 (ω i = Im ω f )}-curves where Im [ω f -ω 0 (X)] > 0. Hence the branch cut starting at X 1 (ω f ) is located above the {X 1 (ω i = Im ω f )}-curve and the branch cut starting at X 2 (ω f ) is located below the {X 2 (ω i = Im ω f )}-curve (see also figure 2). As indicated in §3, the spatial branches ( 6) may be unambiguously interpreted as downstream k + or upstream k -branches when Im [ω f -ω 0 (X)] > 0, which corresponds to the domain of the X-plane located between the X 1 and X 2 curves. The above choice of branch cuts thus extends the definition of the k ± (X, ω f ) wavenumber branches to the entire complex X-plane for any frequency ω f with Im ω f > Im ω s (region above thindashed line in figure 1b). When considering real forcing frequencies ω f in ( 6), the choice of k + and k -branches corresponds to the downstream and upstream branches derived by traditional causality considerations in the stable or convectively unstable domain of the complex X-plane, characterized by Im ω 0 (X) < 0 and located between the solid {X 1,2 (ω i = 0)}-curves of figure 1.

Signalling problem in spatially developing medium

For the signalling problem, the system is driven at some position on the real axis, say X f , with a real frequency, say ω f , and the spatial response to this localized time-harmonic forcing is considered. Since the medium is assumed to be globally stable, any switch-on transients will eventually decay and it is legitimate to seek a long-term response of the form ψ(x, t) = φ(x)e -iω f t , globally tuned to the forcing frequency.

Under the slowly-varying medium hypothesis characterized by ǫ ≪ 1, time-periodic solutions may be obtained in terms of wkbj approximations [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF] where the spatial structure is described by a rapidly varying complex phase, accounting for the local wavelength and spatial growth/decay rate, and a slowly varying envelope. For a given global frequency ω f , a wkbj approximation is of the form

ψ(x, t) = A(X) exp i ǫ X k(u, ω f )du -iω f t (7)
where the envelope A(X) is expanded in powers of ǫ as A ∼ A 0 + ǫA 1 + . . . and governed by amplitude equations that can be computed recursively up to any order. In [START_REF] Davies | Global behaviour corresponding to the absolute instability of the rotating-disk boundary layer[END_REF], the local wavenumber k(X, ω f ) satisfies the local dispersion relation and follows one of the two spatial branches [START_REF] Cross | Pattern formation outside of equilibrium[END_REF]. Hence, in different domains of the complex X-plane, the solution may be approximated by different wkbj expansions of the form (7) pertaining to different spatial wavenumber branches.

The above wkbj expansions are singular [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF] at each of the two turning points

X t = X 1,2 (ω f ) of the dispersion relation, where k + (X t , ω f ) = k -(X t , ω f ).
These turning points are located on their respective {X 1,2 (ω i = 0)}-contours, see figure 2. From a turning point, three Stokes lines emerge, defined by Im

X X t [k + (u, ω f ) -k -(u, ω f )]
du = 0, and partition the complex plane into three different sectors. Along these Stokes lines both wkbj approximations remain of the same order of magnitude, while inside the sectors one approximation is exponentially large with respect to the other.

Consider first the turning point X t 1 ≡ X 1 (ω f ), located on the upper {X 1 (ω i = 0)}-contour in the complex X-plane (upper dashed curve in figure 2), and the two associated linearly independent wkbj approximations

Ψ ± 1 = A ± (X) exp i ǫ X X t 1 k ± (u, ω f )du -iω f t . (8)
Three Stokes lines and one branch cut emanate from X t 1 . The branch cut (dotted line), along which k + and k -branches get exchanged, is located above the {X 1 (ω i = 0)}-curve and extends from X t 1 towards

Consider now the turning point X t 2 ≡ X 2 (ω f ), located on the lower {X 2 (ω i = 0)}-contour (lower dashed curve in figure 2). The same reasoning as above holds for the expansion of the spatial response in terms of the two associated wkbj approximations

Ψ ± 2 = A ± (X) exp i ǫ X X t 2 k ± (u, ω f )du -iω f t . ( 12
)
The associated branch cut (dotted curve) is now located below the {X 2 (ω i = 0)}-curve. With Stokes lines and sectors around X t 2 labeled as in figure 2 (left zoom), the solution to the signalling problem is then approximated by

ψ ∼    C 2 Ψ + 2 in S 2+ 2 , S 3 2 , S 1 2 and along ℓ 1 2 , ℓ 2 2 , C 2 (Ψ + 2 + Ψ - 2 ) along ℓ 3 2 , C 2 Ψ - 2 in S 2- 2 , ( 13
)
and crossover between the dominant k + and k -wavenumber branches takes place across the Stokes line ℓ 3 2 , located opposite the sector containing the region near X = +∞. In the complex X-plane, the solution is thus dominated by the k + branch everywhere except in the two sectors issuing from the turning points X t 1 and X t 2 respectively and delimited by the Stokes lines ℓ 3 1 or ℓ 3 2 and the corresponding branch cuts; these regions are indicated in grey in figure 2. Note that, although the approximations [START_REF] Monkewitz | The absolute and convective nature of instability in two-dimensional wakes at low Reynolds numbers[END_REF] are formally identical to [START_REF] Huerre | Open shear flow instabilities[END_REF], the spatial branches in [START_REF] Monkewitz | The absolute and convective nature of instability in two-dimensional wakes at low Reynolds numbers[END_REF] are integrated from X t 2 and the dominant/subdominant character depends on the behaviour for X radiating away from this turning point.

The relationship between the constants C 1 and C 2 is readily derived by comparing the wkbj expansions ( 11) and ( 13) in regions where they overlap, for example for X → +∞. From

ψ ∼ C 1 A + (X) exp i ǫ X X t 1 k + (u, ω f )du -iω f t and ψ ∼ C 2 A + (X) exp i ǫ X X t 2 k + (u, ω f )du -iω f t it follows that C 2 = C 1 exp i ǫ X t 2 X t 1 k + (u, ω f )du -iω f t . (14)
The above results [START_REF] Dizès | Linear global modes in spatially developing media[END_REF][START_REF] Pier | Open-loop control of absolutely unstable domains[END_REF][START_REF] Pier | Primary crossflow vortices, secondary absolute instabilities and their control in the rotating-disk boundary layer[END_REF] entirely specify the asymptotic wkbj approximation to the unique solution of frequency ω f satisfying a causal boundary condition, i.e. made up of the k + branch near X = +∞. The spatial response of system (1) to harmonic forcing of frequency ω f thus follows this solution over the real X-axis on the right side of the forcing location, for X f < X < +∞. When the external forcing is applied with O(1) amplitude, the response is also O(1) near X f which yields the magnitude of the remaining unknown constant as

C 1 ∼ exp i ǫ X t 1 X f k + (u, ω f )du.
Note however that in this linear setting the precise value of C 1 , and thus the exact spatial location of the forcing, does not influence the spatial structure of the response to the right side of the forcing.

The spatial response pertaining to the domain -∞ < X < X f on the left side of the forcing is similarly derived from the unique solution of frequency ω f made up of the subdominant k -branch for X → -∞. The corresponding details are straightforward and will not be presented here. In a typical situation of interest, forcing is applied upstream of the au domain, i.e., X f < X ca , so that all relevant features concerning the structure of the spatial response over the au interval are derived from the solution valid over X f < X < +∞.

Discussion

The previous analysis yields the spatial response to harmonic forcing over the entire X-axis, including the au interval X ca < X < X ac .

For harmonic forcing of real frequency ω f applied at some station X f upstream of the au domain, local causality considerations suggest that the spatial response follows the k + (X, ω f ) wavenumber display a jump somewhere in the au interval X ca < X < X ac if the k + branches prevailing for both X f < X < X ca and X ac < X < +∞ are to be reconciled. In this situation indeed, the real X-axis is also crossed by the ℓ 3 2 Stokes line and, as shown in the previous section, an exchange in the dominant wkbj solution approximating the spatial response precisely takes place across this Stokes line. Hence the region where the spatial response is dominated by the k -branch (shaded regions in figures 3b,c) overlaps the real X-axis, and along the au interval both the branch cut and the ℓ 3 2 Stokes line are encountered. As a result, the spatial response to harmonic forcing is then obtained in terms of both wkbj approximations and the jump of the dominant local wavenumber occurs at the intersection of the ℓ 3 2 Stokes line with the real X-axis.

Confirmation by direct numerical simulation

The structure of the spatial response derived by analytic methods in the previous sections has been confirmed by direct numerical simulations of system [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF]. The results presented in figure 4 have been obtained with a spatially varying local absolute frequency (5) defined by ω s = 2 -0.5i, X s = 2 + i and ω XXs = -1 -2i and ǫ = 0.1, while the parameters k 0 = 1 -i and ω kk = 1 -0.5i were kept at constant values. These parameter settings correspond to an au interval characterized by X ca ≃ 1.6, ω ca 0 ≃ 3.2 and X ac = 3.3, ω ac 0 ≃ -1.1. Harmonic forcing is applied upstream at X f = 1 with frequencies ω f = 5 (left sequence of figure 4) and ω f = 2 (right sequence). The thick curves in figure 4 illustrate the different characteristics of the numerically computed spatial response: envelope |ψ| and real part ψ r (first row), envelope on a logarithmic scale (second row), real part k r (third row) and imaginary part k i (fourth row) of local wavenumber. The thin curves in the wavenumber plots represent the wavenumber branches k ± (X, ω f ) analytically computed via [START_REF] Cross | Pattern formation outside of equilibrium[END_REF] with the corresponding ω f .

The forcing frequency ω f = 5 lies outside the range ω ca 0 -ω ac 0 , so that the associated turning points X t 1 ≃ 2.7 + 2.5i and X t 2 ≃ 1.3 -0.5i are located on either side of the real X-axis, which is thus free from any branch cuts. Figure 3(a) has been computed with exactly these parameter values. The critical Stokes lines ℓ 3 1 and ℓ 3 2 do not cross the real axis and the spatial response for X > X f is thus predicted to be entirely made up of the k + wavenumber branch. Inspection of the numerical results given in figures 4(a1-a4) reveals that this is indeed the case. From figures 4(a1,a2) it is seen that the spatial response grows from X f through the entire au interval to reach maximum amplitude beyond X ac , before eventual decay further downstream. The local wavenumber in the spatial response, numerically computed as -i∂ x ψ/ψ and corresponding to the thick curves in figures 4(a3,a4), very closely follows the analytically computed wavenumber branches (thin curves). The salient feature is that the spatial reponse follows the k + branch in the entire domain to the right of the forcing location X f , inluding the au interval X ca < X < X ac (delimited by dashed vertical lines).

The forcing frequency ω f = 2 lies within the range ω ca 0 -ω ac 0 , so that the associated turning points X t 1 ≃ 1.8 + 1.7i and X t 2 ≃ 2.2 + 0.3i are both located in the upper complex half plane, as illustrated in figure 3(b) obtained with these parameter values. Now the Stokes line ℓ 3 2 crosses the real axis for X 3 2 ≃ 1.9 where a jump in the dominant wavenumber of the spatial reponse is predicted. This phenomenon is indeed observed in figures 4(b3,b4). The spatial response is seen to be dominated by the k + branch from the forcing location X f to the Stokes line at X 3 2 , beyond which the k -branch dominates. This exchange in dominant wkbj approximations at X 3 2 materializes by a cross-over of the numerically computed local wavenumber in figures 4(b3,b4) and by a kink of the envelope in figure 4(b2). The spatial response then follows the k -branch from X 3 2 to the branch cut at X ⋆ ≃ 3.0 (vertical dotted line), beyond which the k -branch is seemlessly relabelled as k + . As anticipated from the previous section, the spatial response displays the k + wavenumber branch in the two at most cu domains X f < X < X ca and X ac < X < +∞. However, with ω f in the range ω ca 0 -ω ac 0 these branches are not continuously connected through the central au interval. Hence the solution is there made up of both branches and cross-over of the dominant wavenumber occurs at the intersection with the Stokes line.

Conclusion

The signalling problem in spatially developing systems is fairly well understood for media that are at most convectively unstable, and it is generally believed to be ill-posed for absolutely unstable

X f X ca X ac (a1) x |ψ| ψ r (a2) x log |ψ| (a3) x k r k + r k - r k + r (a4) x k i k + i k - i k + i k - i X f X ca X 3 2 X ⋆ X ac (b1) x |ψ| ψ r (b2) x log |ψ| (b3) x k r k + r k - r k - r k + r k + r (b4) x k i k + i k - i k + i k - i k + i k - i
Fig. 4 Envelope |ψ|, real part ψr, local wavenumbers kr and ki of spatial response to harmonic forcing obtained by direct numerical simulation. Forcing is applied at X f , upstream of the absolutely unstable interval X ca < X < X ac . (a) With ω f outside the range ω ca 0 -ω ac 0 , the spatial response follows the k + wavenumber branch for X f < X < +∞, including the au interval. (b) For ω f within the range ω ca 0 -ω ac 0 , the real X-axis intersects the Stokes line ℓ 3 2 at X 3 2 , where a jump in the dominant wavenumber is observed. Beyond the branch cut at X⋆, the k -wavenumber branch prevailing for X 3 2 < X < X⋆ becomes k + and extends to X = +∞.

systems. The present investigation has shown that this problem remains well-posed for a certain class of absolutely unstable systems: whenever the system is globally stable (i.e. all global modes decay with time) all transient behaviour is eventually damped and a spatial response tuned at the forcing frequency establishes for large time. By resorting to asymptotic wkbj expansions, the spatial structure of the time-harmonic response has been obtained in the entire domain, including the au interval. In the au interval, the dominant local wavenumber of the spatial response has been shown to depend on the positions in the complex X-plane of the turning points associated with the forcing frequency. The positions of these turning points depend on the relative values of the forcing frequency ω f and the (real) marginal absolute frequencies ω ca 0 and ω ac 0 prevailing at the end points of the au interval. The detailed analysis in the complex X-plane of the Stokes lines and branch cuts has revealed that two distinct situations arise:

-For forcing frequencies ω f outside the range ω ca 0 -ω ac 0 , the local wavenumber observed in the spatial response continuously depends on streamwise distance and follows the k + -branch everywhere downstream of the forcing location.

-For forcing frequencies within the range ω ca 0 -ω ac 0 , the spatial response displays a jump in the dominant local wavenumber branch at the location within the au domain where the real X-axis is crossed by a Stokes line.

The above results have been analytically derived in a general setting and confirmed by numerical integration of the Ginzburg-Landau model equation. Work currently in progress concerns application of the present findings to the spatial response in the rotating-disk three-dimensional boundary-layer and its implications for an open-loop control strategy [START_REF] Pier | Open-loop control of absolutely unstable domains[END_REF][START_REF] Pier | Primary crossflow vortices, secondary absolute instabilities and their control in the rotating-disk boundary layer[END_REF], based on localized harmonic forcing.

Fig. 2 Turning points and Stokes lines in the complex X-plane. The two turning points X t 1,2 are located on their respective {X1,2(ωi = 0)}-isocontours (dashed lines). Each turning point is associated with a branch cut (dotted curves) and gives rise to a network of three Stokes lines (solid curves).

infinity in the upper complex X-plane. The three Stokes lines ℓ 1 1 , ℓ 2 1 and ℓ 3 1 divide the complex plane into three sectors S 1 1 , S 2 1 and S 3 1 , and sector S 2 1 is further divided into S 2- 1 and S 2+ 1 by the branch cut, see zoom to right of figure 2. Along these Stokes lines both approximations Ψ ± 1 are of the same magnitude, and the dominant (resp. subdominant) solution becomes subdominant (resp. dominant) whenever a Stokes line is crossed.

Since the region near X = +∞ is at most convectively unstable, the spatial response to the external forcing there necessarily features the k + wavenumber branch and is thus made up of the subdominant Ψ + 1 approximation. Hence the dominant approximation Ψ - 1 is not present in the sector S 3 1 , defined as the sector containing the region near X = +∞, and the solution of the signalling problem in S 3 1 and along its bordering Stokes lines ℓ 1 1 and ℓ

Similarly, when continuing the solution from S 3 1 into S 2+ 1 across ℓ 1 1 , approximation Ψ + 1 becomes dominant and a solution of the form

When crossing the branch cut, the approximations Ψ + 1 and Ψ - 1 get exchanged so that the same wkbj superposition reads

1 , with Ψ - 1 now the dominant term. On the Stokes line ℓ 3 1 , both approximations are again of the same order, and the solution is then approximated by

Finally, comparing both expansions ( 9) and [START_REF] Huerre | Hydrodynamic instabilities in open flows[END_REF] shows that all coefficients C = 1 , C - 1 and C + 1 are identical to, say, C 1 . In the different regions around turning point X t 1 , the following approximations to the signalling problem are then obtained

Note that nothing special happens across the branch cut, except relabeling of the superscripts of k ± and Ψ ± 1 . Hence a result similar to the above holds even in situations where the orientation of the Stokes lines is such that the branch cut is not confined between ℓ 1 1 and ℓ 3 1 (see example in figure 3a). Thus it is only across the Stokes line ℓ 3 1 , defined as the Stokes line opposite the sector containing the region near X = +∞ (see also figure 3 for further examples), that a change in the dominant wkbj approximation occurs: a crossover between k + and k -as dominant wavenumber branches takes place across ℓ 3 1 while everywhere else the dominant local wavenumber continuously depends on X.

(a)

Fig. 3 Spatial structure of the signalling problem in the complex X-plane for different forcing frequencies ω f . The two turning points X t 1,2 (thick dots) are located on their respective {X1,2(ωi = 0)}-isolines (dashed curves) and give rise to a network of Stokes lines (solid curves) and branch cuts (dotted lines). The spatial response is dominated by the k -wavenumber branch in the shaded regions and by the k + branch everywhere else. The au interval of the real X-axis is hatched. (a) With ω f outside the range ω ca 0 -ω ac 0 , the turning points X t 1,2 are located in opposite half X-planes and the k + wavenumber dominates over the entire real X-axis. (b,c) With ω f within the range ω ca 0 -ω ac 0 , the turning points X t 1,2 are located in same half X-planes, and the spatial response displays a jump in the dominant wavenumber on the real axis at the intersection with the ℓ 3 2 Stokes line.

branch over the stable or cu interval X f < X < X ca , between the forcing location and onset of absolute instability. Similarly, in the at most convectively unstable downstream domain extending from X = X ac to X → +∞, the solution is obtained as the downstream response to some upstream located cause and is there also expected to follow the k + (X, ω f ) wavenumber branch. However, local considerations fail to predict the structure of the spatial response in the au interval X ca < X < X ac since the k ± wavenumber branches ( 6) are there only formal solutions of the local dispersion relation and cannot be interpreted in terms of upstream or downstream spatial branches.

Based on the developments of the previous section, it will now be shown that in the au region X ca < X < X ac the structure of the spatial response and the dominant wavenumber essentially depend on the positions of the turning points X t 1,2 ≡ X 1,2 (ω f ) relative to the real X-axis and thus on the relative values of ω f , ω ca 0 and ω ac 0 . Indeed, a forcing frequency ω f outside the range ω ca 0 -ω ac 0 lies above the {ω 0 (X i = 0)} curve in the complex ω-plane (dotted parabolic curve of figure 1b). The associated turning points are then located on either side of the real X-axis: X t 1 in the upper half-plane and X t 2 on parts of the {X 2 (ω i = 0)}isoline in the lower half-plane, as illustrated in figure 3(a). In this situation, the branch cuts emanating from X t 1,2 and respectively extending into the far upper and lower half X-planes do not cross the real X-axis. As a result, the spatial branches k ± (X, ω f ) continuously depend on X over the entire real axis, and the k + (X, ω f ) branches prevailing on either side of the au interval (hatched in figure 3) are continuously connected over the real X-axis. Note that the branch cut starting at X t 2 , while remaining below the {X 2 (ω i = 0)}-isoline, could cross the real axis twice before heading toward the far lower X-plane; however, the conclusions are not affected as long as an even number of crossings occurs since the branch cut could then be modified so as to avoid the X-axis. At the same time, the Stokes lines ℓ 3 1 and ℓ 3 2 over which an exchange of the dominant wkbj approximations occurs do not cross the real X-axis either. The regions in the complex X-plane where the spatial response is dominated by the k - branch (shaded regions in figure 3a) then do not overlap the real X-axis. Hence the spatial response to harmonic forcing is approximated by the wkbj approximation featuring the k + wavenumber branch over the entire domain X f < X < +∞ including the au interval X ca < X < X ac .

By contrast, for a forcing frequency ω f within the range ω ca 0 -ω ac 0 , i.e. below the {ω 0 (X i = 0)}-curve in the complex ω-plane, both turning points are located on the same side of the real X-axis: turning point X t 2 has moved to the part of the {X 2 (ω i = 0)}-isoline located in the upper complex half-plane, as illustrated in figures 3(b,c). Thus the branch cut emanating from X t 2 crosses the real X-axis and the spatial k ± (X, ω f ) branches display a discontinuity across the branch cut. The k + (X, ω f ) branches prevailing on either side of the au interval are no longer continuously connected over the real X-axis. It follows that the local wavenumber of the spatial response to harmonic forcing applied at X f must


François Willot
email: francois.willot@cmm.ensmp.fr

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

Elastic behavior of composites containing Boolean random sets of inhomogeneities

Keywords: Homogenization, Bulk modulus, Representative volume element, Integral range, Boolean model of spheres, Linear elasticity, Random media PACS: 46.04.+b, 46.15.-x, 46.65.+g

The overall mechanical response as well as strain and stress field statistics of an heterogeneous material made of two randomly distributed, linear-elastic phases, are investigated numerically. The Boolean model of spheres is used to generate microstructures consisting of either porous or rigid inclusions, at any volume fraction of the phases. Stress and strain field integral ranges, or equivalently the representative volume element, are computed and linked to features of the field statistics, and to the microstructure geometry.

Introduction

This work is concerned with the overall (i.e. mean) elastic response and field statistics of an heterogeneous medium subjected to an applied macroscopic deformation. The usual approach for the prediction of the elastic properties of random media is based on the calculation of bounds from a limited amount of information on the microstructure. Hashin and Shtrikman [START_REF] Hashin | A variational approach to the theory of the elastic behavior of multiphase materials[END_REF] rely on a specific variationnal principle and on the assumption of statistical morphological isotropy of the medium. Beran and Molyneux (1) worked out third order bounds of the bulk modulus based on the three-points statistics of random media. These bounds were specialised to the case of two-phase media by Milton [START_REF] Milton | Bounds on the elastic and transport properties of two component composites[END_REF], and to some specific models of random media (26; 10). However, bounds are of limited usefulness when the contrast between components of an heterogeneous medium is too large (namely for porous media or for rigid components). In the absence of any general analytical formula to predict the overal elastic properties of random media, an alternative way is based on the use of numerical techniques to solve the elasticity problem on images of the microstructure. Recent progress in computer science and in numerical implementation of the solution of partial differential equations make it possible to operate on 3D images. However, the size of the images that can be handled still remains (and will remain) limited, so that must be addressed the question of the representativity of the effective properties estimated by this approach (14; 11; 13).

In this paper, we consider the numerical homogenization of the bulk modulus of the Boolean model of spheres, in the case of porous or rigid inclusions. We first give a short reminder of the Boolean model and of its periodic version. In a second step, the numerical technique used to estimate the stress field, based on iterations of Fast Fourier transforms, is detailed. It is validated in the case of a periodic array of spheres, for which an accurate estimate is known (3; 4). Results obtained for Boolean models of spheres, covering the effective properties as well as some statistical properties of the field are given in the case of porous and rigid inclusions.

The Boolean model of spheres in linear elasticity

General framework

This study is focused on a medium made of two linear-elastic phases, which consists of spherical overlapping inclusions (phase 2) embedded in a "matrix" (phase 1). The limit of infinitely-contrasted properties between the phases is considered, i.e. the material is either reinforced with rigid particles or weakened by pores. More precisely, the microstructure is generated from a Boolean model of spheres [START_REF] Matheron | Random sets and integral geometry[END_REF], with the number of sphere centers in any given region following a Poisson point process. The intensity of the Poisson process is chosen so that the volume fraction of the inclusions is equal, on average, to the prescribed volume fraction. Consequently, the volume fraction is different for each configuration, except in the limit of an infinite sample. This work is restricted to the case of an isotropic macroscopic deformation, which yields results pertaining to the effective bulk modulus, whereas the effective shear modulus is, on average, null. It is assumed that the applied loading is small enough that it induces very small matter displacement (compared to the smallest typical length of the microstructure), such that the phase geometry may be fixed for each given sample, and the homogenization problem be treated within the theory of small perturbations.

Given a finite sample of the medium, materialized by a cube of size L, an overall strain loading with periodic boundary conditions (defined below) is applied. Previous results [START_REF] Kanit | Apparent and effective physical properties of heterogeneous materials: Rep-resentativity of samples of two materials from food industry[END_REF] obtained for random linear composites have shown that, as far as the effective elastic moduli are concerned, faster size convergence is achieved with periodic boundary conditions than with classical "kinematic" or "static" types boundary conditions, which tend to perturbate fields along boundaries. To take advantage of this property, a Boolean model suited to periodic media is used i.e. each randomly-chosen sphere center, or seed, is periodicized before being replaced with spheres.

Constitutive law

The constitutive law relates at each point x in the medium the stress and strain tensors σ and ε. In the matrix, it is characterized by its bulk and shear elastic moduli, noted repectively κ and µ, and is defined by:

σ(x) = L : ε(x), L ijkl = κδ ij δ kl + µ(δ ik δ jl + δ il δ jk ) ( 1
)
where L is a fourth-order elastic tensor, the double point : denotes a double contraction between its right and left-hand terms, and δ is the Kronecker symbol. In the following, the shear and bulk moduli in the matrix have been chosen as: κ = 1/3, µ = 1/2. It is assumed that the results of this work apply qualitatively to more general materials as long as both elastic moduli in the matrix are finite. Expanded in terms of the mean and deviatoric strain and stress components:

σ m = 3κε m , σ ′ = 2µε ′ , ε ′ ij = ε ij -ε m δ ij , σ ′ ij = σ ij -σ m δ ij , (2)
where the mean strain and stress are the scalars noted ε m = ε ii /3 and σ m = σ ii /3, respectively, whereas the deviatoric parts ε ′ and σ ′ are tensors. In phase 2, the constitutive law in porous (resp. rigid) inclusions is σ ≡ 0 (resp. ε ≡ 0) with strain (resp. stress) undefined.

Stress equilibrium and strain compatibility

The strain must derive from a displacement field u (referred to as the "strain compatibility"condition). Assuming small deformation, this reads:

ε ij = 1 2 (δ i u j + δ j u i ), (3)
Symetrically, the stress field satisfies an "equilibrium equation", which reads, in the absence of external forces:

∂ i σ ij = 0. ( 4
)
Note that interfaces between inclusions and matrix are perfect, i.e. the stress boundary conditions along matrix-pores interfaces read

σ • n = 0, ( 5
)
where n is a vector normal to the boundary. In the case of a rigid inclusion, the boundary condition along interfaces read, in terms of the displacement field:

u • n = 0, (6)
Assuming sufficiently smooth interfaces (i.e. without angles or corners), Eq. ( 5) may be derived using ( 4), (3), writing each derivative as a small difference between two points accross the interface. However, in the numerical computations undertaken in this work, it is easier to appeal to the constitutive law and to the equilibrium and compatibility equations, defined at each point in the heterogeneous material, rather than interface boundary conditions.

Full-fields computations by means of Fast Fourier Transforms

The "Fast Fourier method" introduced by Moulinec, Suquet ( 22) is a breakthrough in numerical methods for computing the mechanical response of composite materials. It has successfully been applied to many other homogenization problems, e.g. conductivity or permittivity [START_REF] Jeulin | Moreaud Statistical representative volume element for predicting the dielectric permittivity of random media[END_REF]). It makes use of an iterative fixed-point algorithm derived from the periodic Lippman-Shwinger equation (see e.g. Kroner ( 16)). Given an homogeneous elastic reference material of elastic tensor L 0 , the strain field ε is a solution of the implicit equation:

ε(x) = ε -d 3 x ′ G 0 (x -x ′ ) : τ (x ′ ), τ (x) = σ(x) -L 0 : ε(x), (7)
where the Green function G 0 depends only on the reference tensor L 0 , and its average is chosen equal to zero. The polarization field with respect to the reference L 0 is noted τ and ε is the prescribed macroscopic strain field. The original scheme of the FFT method is a straightforward application of the Lippman-Schwinger equation. Given an image of the strain field (e.g. ε ≡ 0), the polarization field τ on the right side equation in ( 7) is computed at each point in the real space, making use of the constitutive law, whereas the convolution on the left side equation is computed in the Fourier space, as a product, yielding a new strain field. The original FFT algorithm consists in iterating the two equations, computing forward and backward fast Fourier transforms at each step, on a cubic grid of L × L × L voxels. Iterations are stopped when the given convergence criterion is lower than a fixed value. Authors commonly check for stress equilibrium in the Fourier space e.g.:

q |q • σ(q)| 2 1/2 |σ(0)| < η, (8)
where η is typically of order 10 -5 . The use of an optimal reference elastic tensor L 0 strongly improves convergence, however the number of iterations to achieve convergence is roughly proportional to the contrast between the phases, so that this method can not be applied to infinitely contrasted (e.g. porous) media.

Subsequently, Eyre and Milton (8) and Michel et al. [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] have developped refined algorithms which proved more efficient in the case of a high contrast between the phases (see Moulinec,Suquet (24) for comparisons between the three methods). In this work, we use Michel's "agmented Lagrangian" algorithm, which, among its advantage, is the only one to converge in the limit of an infinitely contrasted medium. However, no formula for the optimum elastic reference tensor L 0 is known. Worst, numerical computations show that the optimum choice is very sensitive to the contrast, geometry and the required precision, although less to the system size. Consequently, a sub-optimal reference field is numerically estimated at small grid sizes, and applied to larger systems.

The Augmented Lagrangian scheme is derived from the minimization of an elastic potential over an auxiliary variable e with the additional condition that e = ε (see Eq. ( 6) in Moulinec, Suquet [START_REF] Moulinec | Comparison of FFT-based methods for computing the response of composites with highly contrasted mechanical properties[END_REF]). Contrary to the strain field ε in the the original scheme, the compatibuility for e is not enforced at each step. Since the Fourier transform of the stress is not directly available in the augmented Lagrangian algorithm, convergence is checked as a relative difference in the real space between the tensors ε and e:

x ||ε(x) -e(x)|| 2 1/2 /||σ(0)|| < η, (9)
with η = 10 -6 or η = 10 -7 . Attention must be paid to the highest frequencies (e.g. q = ±(L/2 + 1), where L is the number of voxels along each direction, see Moulinec, Suquet ( 23)), at which point the Green function lacks its usual symmetry (i.e. G * (q) = G(-q)), inducing an edge effect. To check this effect, we have compared the strain fields as computed with two different methods "FFT-I"and "FFT-II"that ensure, respectively, σ(q) = 0 and G(q) = 0 at the highest modes. In Fig. 1, both numerical methods are compared to 2D exact solutions recently derived (Willot, Pellegrini, Ponte ( 27)) in the case of a strongly anisotropic medium weakened by an array of porous cylinders.

As expected, the difference between results obtained with the two methods "FFT-I" and "FFT-II" is negligible compared to the size effect. The exact solution is recovered when size increases, with high oscillations observed at low resolution.

In the case of a rigid phase, high but finite values of the local elastic moduli must be introduced. In the numerical computations, we have chosen κ ′ = µ ′ = 10 3 ≫ κ, µ as the local bulk and shear moduli in "rigid" inclusions. In pores, however, κ ′ = µ ′ = 0.

Homogenization and overall loading

Given a sufficiently large sample of a microstructure, the material behavior is governed by effective properties which relate the macroscopic strain and stress tensors. When all phases are locally linear elastic, the effective constitutive law is also linear elastic:

σ = L 0 : ε , (10)
where • denotes an average over the sample. Conversely, a given sample of a material is said to be statistically representative of a microstructure if it is large enough to satisfy Eq. [START_REF]Mechanics of Random and Multiscale Microstructures[END_REF]. Consistently with Eq. ( 10), we apply a macroscopic strain field ε 0 = ε in a given direction, and derive the effective elastic tensor from the resulting macroscopic stress field σ 0 = σ . Averages are computed from the complete images of the strain and stress fields that satisfy periodic boundary conditions, equilibrium and compatibility as well as the microscopic strain-stress law.

Validation on a periodic array of spheres

In this section, we investigate the accuracy of the FFT method against an analytical approximation (3) of the effective elastic moduli of a cubic array of spheres which draws on previously developped methods by Bergman(2). This estimate gives, in terms of the sphere volume fraction f , an analytical formula for the effective bulk modulus which is valid up to order 3 in f . First, the effective elastic moduli is numerically computed at a high level of accuracy, up to 8 digits, for finite contrast, with the FFT tool (we have chosen elastic moduli in the spheres that are two times higher than in the matrix). Comparing these results with analytical estimates, we have checked that the difference between the Bergman and FFT computations is of order ∼ f 4 (see Fig. 2). Secondly, at infinite contrast (here, with porous inclusions), the Bergman estimate compares well to FFT results for the effective bulk and elastic shear modulus (Fig. 3).

6 Results

Effective moduli

In this section and hereafter, the Boolean model of spheres is considered for the microstructure. The effective bulk modulus as a function of the fraction f of porous and rigid inclusions are given in Figs. ( 4) and [START_REF] Delarue | Prévision du comportement électromagnétique de matériaux composites à partir de leur mode d'élaboration et de leur morphologie[END_REF]. For comparison the self consistent [SC] estimates, as well as standard Hashin and Shtrikman [HS] bounds, and tighter third-order Beran bounds [TOB] are provided. It appears that in the present case, bounds are not useful, being too much separated. In addition, the self-consistent estimate, does not provide a pertinent estimate over the entire range of volume fractions.

The effects of percolation appear clearly on these figures: the bulk modulus vanishes for a pore volume fraction at f = f ′ ≈ 0.95 in Fig. [START_REF] Cohen | Simple algebraic approximations for the effective elastic moduli of cubic arrays of spheres[END_REF]. According to various studies the geometrical percolation threshold of the matrix of a Boolean model ranges from 0.946 [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF] to 0.97 [START_REF] Elam | Critical properties of the void percolation problem for spheres[END_REF]. In the porous Boolean medium, we observe the scaling law: κ 0 /κ ∼ 6.5(f ′ -f ) 3.0 (see Fig. 4), for

|f ′ -f | < 0.1.
Inversely, in the rigid Boolean model, the bulk modulus grows according to a power law beyond the percolation threshold (close to f = f ′′ ≈ 0.29, see e.g. (25; 12)). It turns out that our results for f > f ′′ may be fited with a good approximation to a power law dependence on f , for almost the full range of volume fractions beyond percolation. Normalizing the effective bulk modulus with the bulk modulus in the spheres, we find κ 0 /κ 2 ≈ 2.39(f -f ′′ ) 2.35 for f > f ′′ ≈ 0.29. Interestingly, power law dependencies that extend well beyond the percolation region have been found [START_REF] Deptuck | Measurement of elasticity and conductivity of a three-dimensional percolation system[END_REF] for the Young modulus of a different medium.

Fluctuations of effective properties

In this section the representativity of the effective property estimated on a bounded domain of a microstructure is addressed. To compute the so-called "Representative Volume Element" (14; 11) (RVE) we make use of a geostatistical approach based on the experimental determination of the integral range (17) from numerical simulations.

The integral range

We consider fluctuations of average values over different realizations of a random medium inside the domain B with the volume V . In geostatistics, it is well known that for an ergodic stationary random function Z(x), with mathematical expectation E(Z), one can compute the variance D 2 Z (V ) of its average value Z(V ) over the volume V as a function of the central covariance function C(h) of Z(x) (17) by :

D 2 Z (V ) = 1 V 2 V V C(x -y) dxdy, (11)
where

C(h) = E{(Z(x) -E(Z)) (Z(x + h) -E(Z))}
For a large specimen (with V ≫ A 3 ), equation ( 11) can be expressed to the first order in 1/V as a function of the integral range in the space R 3 , A 3 , by

D 2 Z (V ) = D 2 Z A 3 V ( 12
)
with

A 3 = 1 D 2 Z R 3 C(h) dh (13)
where D 2 Z is the point variance of Z(x) (computed by taking averages from numerical estimates) and A 3 is the integral range of the random function Z(x), defined when the integral in equations ( 11) and ( 13) is finite. The asymptotic scaling law [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF] is valid for an additive variable Z over the region of interest V .

In the case of a two-phase material modeled by a random set A with the characteristic function k(x) (with the volume fraction p), the point variance D 2 Z of the random function k(x) is given by p(1-p), and the central covariance function is obtained from the set covariances C(h

) = P {x ∈ A, x + h ∈ A} or Q(h) = P {x ∈ A c , x + h ∈ A c } by : C(h) = C(h) -p 2 = Q(h) -q 2 , with q = 1 -p.
To estimate the effective elasticity tensor L 0 defined by Eq. 10, we have to compute the averages σ or ε . For the applied boundary conditions (a hydrostatic loading in the present case), the bulk modulus is obtained from the estimations of a scalar, namely the mean stress field. Therefore the variance of the effective property κ 0 follows the equation ( 12) when the integral range A 3 of the relevant field is known.

Since the theoretical covariance of the fields (σ or ε) is not available, the integral range can be estimated according to the procedure proposed by Georges Matheron for any random function [START_REF] Matheron | Estimating and Choosing[END_REF]: working with realizations of Z(x) on domains B with an increasing volume V (or in the present case considering subdomains of large simulations, with a wide range of sizes), we can estimate the parameter A 3 by fitting the obtained variance according to the expression [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF]. The point variance D 2 Z of the corresponding field is directly estimated from the experimental variance of the field.

Practical determination of the size of the RVE

When considering a material as a realization of a random set or of a random function RF, the idea that there exists one single possible minimal RVE size must be left out. Instead, the size of a RVE can be defined for a physical property Z, a contrast and, above all, a given precision in the estimation of the effective properties depending on the number of realizations that one is ready to generate. By means of a standard statistical approach, the absolute error ǫ abs and the relative error ǫ rela on the mean value obtained with n independent realizations of volume V are deduced from the interval of confidence by:

ǫ abs = 2D Z (V ) √ n ; ǫ rela = ǫ abs Z = 2D Z (V ) Z √ n (14)
The size of the RVE can now be defined as the volume for which for instance n = 1 realization is necessary to estimate the mean property Z with a relative error ǫ rela = 1%, provided that the function D Z (V ) is known. Alternatively, we can decide to operate on smaller volumes (provided no bias is introduced), and consider n realizations to obtain the same relative error.

In the case of effective elastic moduli, the exact mean value and variance for a given domain size are a priori unknown. The absolute error on the mean value is evaluated using Eq. ( 14). This methodology was applied to the case of the dielectric permittivity of various random media (5; 13), and to the elastic properties and thermal conductivity of a Voronoï mosaic (14).

Covariances and integral range for the Boolean model

In this section, we consider the covariance field of the characteristic function of phase ( 2) in the Boolean model C micro (h) and of the mean stress field C σm (h), estimated from FFT computations. In the case of the characteristic phase, the theoretical covariance of the Boolean model

Q(h) = P {x ∈ A c , x + h ∈ A c } is given by : Q(h) = q 2-r(h) , (15)
where r(h) = K(h) K(0) , K(h) being the geometrical covariogram of the primary grain A ′ of the model (here a sphere with radius a). In the present case we have:

K(h) = V (A ′ ∩A ′ -h ) = π 6 (2a) 3   1 - 3 2 h 2a + 1 2 h 2a 3   for h < 2a, else K(h) = 0 (16)
Starting from Eqs. ( 13) and ( 16), the integral range of the volume fraction of the Boolean model is easily computed as a function of p or q. Alternatively, it can be estimated from the simulations of the microstructure, which gives a validation of the representativity of the simulations.

When a hydrostatic loading is applied, the function C X (h) does not depend on the direction of the vector h and is approximated by computing the mean of the three values C X (|h|e 1 ), C X (|h|e 2 ) and C X (|h|e 3 ) in the three privileged directions of the FFT grid.

In this section we also make use of stress field maps as shown in Figs. [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF] and [START_REF] Jeulin | Moreaud Statistical representative volume element for predicting the dielectric permittivity of random media[END_REF]. They represent random 2D sections of the mean component of the stress field σ m , for the porous and rigid Boolean models, respectively, at various increasing volume fractions f of the spheres. Colors have been adjusted to the maximum and minimum of the field in each section to highlight the field patterns. A detailed study of these fields is out of the scope of the present paper.

Plotted as a function of the distance h, the covariance function exhibits a characteristic length of the microstructure. For a continuous field, it makes conspicious the typical length related to fluctations of the apparent average of a field over subdomains. Applied to the mean stress field, this is the minimal length of the representative volume element. Normalizing the covariance as

C(h) = [C(h) -σ 2 0 ]/[C(0) -σ 2 0 ]
, covariance functions C(h) of the mean stress are plotted as a function of the length h, for various concentrations f of the rigid phase in Fig. 8. From this figure it appears that the correlation length (or range) of the mean component of the stress field (i.e. the length for which the central correlation function reaches its asymptotic value 0) compares well to the diameter of the spheres (18 voxels in the present case) for the lower pore volume fraction (f = 0.25). For higher porosities, the range of the stress field is larger than the range of the microstructure.

Normalized covariance functions C(h) of the mean stress in a Boolean model with spherical rigid inclusions are given in Fig. [START_REF] Elam | Critical properties of the void percolation problem for spheres[END_REF]. Contrary to the porous case, it is observed that the correlation length of the mean component of the stress field is significantly higher to the diameter of the spheres except for volume fractions f lower than the percolation threshold.

The integral ranges obtained from the numerical simulations (using the procedure explained in subsection 6.2.1) are given as a function of the volume fraction of inclusions f in Fig. (9. They are normalised for a sphere radius a equal to 1. Error bars are obtained from the quality of the fit obtained by means of the theroretical change of the variance D 2 Z (V ) with the volume V (Eq. 12).

It appears that the integral range (and therefore the RVE) of the microstructure is smaller than that of the mean stress component in the composites, except in regimes where the softer phase has percolated or is close to the percolation threshold (when f > 0.8 and f < 0.35, for the porous and rigidlyreinforced media, respectively).

For porosities lower than the percolation threshold, the integral range of the mean stress σ m is close, but slightly larger than that of the microstructure. As expected, high values of the mean stress component are concentrated around voids (see e.g. maps (b) and (c) in Fig. 12), inducing anti-correlations and lowering the integral ranges. Such phenomenon is not observed in the rigidly-reinforced media, where the stress field in the rigid inclusions is high. In this material, the integral range of the mean stress is, consequently, much larger than that of the microstructure beyond the percolation threshold. Such large integral ranges have already been observed (13; 14) for another model of microstructure, in cases of finite contrast between the phases.

Field histogram

This section is concerned with the field histograms of the hydrostatic stress component. The field histogram (i.e. probability distribution) for the mean stress in the matrix (phase 1) is defined as:

P σm (t) = 1 1 -f (1) d 2 xδ(σ m (x) -t), (17)
Field histograms are computed numerically by means of the FFT solutions, and are plotted in Figs. [START_REF]Mechanics of Random and Multiscale Microstructures[END_REF] and [START_REF] Jeulin | Random Structures in Physics, in: "Space, Structure and Randomness[END_REF], for the porous and rigid media respectively. Each figure shows a set of curves corresponding to particular values of the concentration f of the spheres. As in the previous section, we appeal to maps of the hydrostatic stress component as given in Figs. [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF] and ( 13), for porous and rigid Boolean models, respectively, to discuss the effects of the morphology of the microstructure on the local distribution of strains and of stresses in the material.

In general, the field distributions P σm of the mean stress component are not Gaussian, and non-symmetric. In the porous material, the field distribution P σm (Fig. 10) takes high values around the average of the mean stress in the matrix σ m = σ m (1) for f < 0.3. For higher porosities, large values of the distribution are found at σ m = 0. At a particular volume fraction f ≈ 0.7 (see e.g. map (d) in Fig. 12), one of the local maxima of P σm disappears, at which point the distribution field is highly non-symmetric. This change of the shape of the distribution corresponds to the appearance of "dead zones" of the matrix where the strain and stress fields are null.

In the rigid material, however, the distribution P σm of the mean stress in the matrix, as shown in Fig. [START_REF] Jeulin | Random Structures in Physics, in: "Space, Structure and Randomness[END_REF], always has a unique maximum. A change of the shape of the distribution is observed near f = 0.3, close to the percolation threshold. At this volume fraction, it is conspicuous that the distribution P σm is most non-symmetric, whereas for f > 0.3, when the "matrix" is made of isolated clusters and is much softer than the spheres, the distribution P σm is close to a Gaussian distribution.

Conclusion

The coupling between 3D simulations of random sets and the numerical resolution of the equations of elasticity by means of FFT iterations is a powerful tool to estimate the effective properties of random media, even in the case of a high contrast between its components.

Using this methodology, we could, for the first time, estimate the variation of the bulk modulus of Boolean models of spheres in the case of porous or rigid inclusions. In addition to the effective properties, fluctuations of the field averages over finite domains and histograms of the mean component of the stress field were investigated. The same approach may be applied to multiscale heterogeneous materials [START_REF] Jeulin | Percolation d'agrégats multi-échelles de sphères et de fibres : Application aux nanocomposites[END_REF], and to multi-component media, allowing a comparison with the prediction of standard models, such as the self-consistent estimates or bounds of effective properties.

The representative volume element (RVE) depends on the properties under consideration. For porous media, for most volume fractions, the size of the RVE for the mean stress component is close to that of the RVE of the microstructure. It is either smaller, or much larger, in the case of rigidly-reinforced media.

For porous materials, bimodal distributions of the mean stress were observed. Investigations regarding the shear modulus of Boolean models are currently undertaken using similar numerical techniques.

-0.

Fig. 1 . 4 Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 14234567 Fig.1. Difference os the strain components ε xxε yy along a half-line y = 1/2 located on the border of the unit cell, for a periodic array of cylinders: comparison between the exact solution in the continuum (solid line) and two FFT methods, "FFT-I" and "FFT-II" with slightly different treatment of the highest frequencies. Two grid sizes are shown for the first method. The material is strongly anisotropic (see Willot, Pellegrini, Ponte (27)) and "pure shear" loading applied (i.e. only the macroscopic components ε xx =ε yy are non zero).

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig.8. Normalized covariance C σm (h) of the mean component of the stress field σ m as a function of the distance h (number of voxels), for a Boolean medium at various different concentration of rigid inclusions f = 0.2, 0.4, 0.6 and 0.8. Along the x-axis, the length h is expressed in voxel units (the diameter of a sphere is 18 voxels). A hydrostatic load is applied as a macroscopic deformation ε m = 1, on a 256 3 FFT grid.

Fig. 11 .

 11 Fig.11. Histograms P σm (t) of the mean component of the stress field σ m in the matrix, in the rigidly-reinforced random medium at increasing volume fractions of the rigid particles f = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.8, 0.9 (right to left, as indicated on the figure). The values of the stress on the x-axis correspond to a hydrostatic macroscopic deformation ε m = 1, on a 256 3 FFT grid.

Fig. 12 .

 12 Fig. 12. 2D sections of the stress component σ m in the random porous medium, at various porosities f = 0.05 (a), 0.25 (b), 0.45 (c), 0.6 (d). A hydrostatic deformation ε m is applied. The color scale goes from black (smallest values) to white (highest values) with red and yellow in-between (color online).

Fig. 13 .

 13 Fig. 13. 2D sections of the stress component σ m in the random rigidly-reinforced medium, at various rigid inclusions concentrations f = 0.1 (a), 0.4 (b), 0.7 (c), 0.9 (d). A hydrostatic deformation ε m is applied. The color scale goes from black (smallest values) to white (highest values) with red and yellow in-between (color online).

The authors are indebted to Institut Carnot M.I.N.E.S. for a financial support during this study.