T Cambonie
email: cambonie@fast.u-psud.fr

J Bares

M L Hattali

D Bonamy

‡ V Lazarus

H Auradou

Effect of the porosity on the fracture surface roughness of sintered materials: From anisotropic to isotropic self-affine scaling

Keywords: number(s): 46, 50, +a

To unravel how the microstructure affects the fracture surface roughness in heterogeneous brittle solids like rocks or ceramics, we characterized the roughness statistics of postmortem fracture surfaces in homemade materials of adjustable microstructure length scale and porosity, obtained by sintering monodisperse polystyrene beads. Beyond the characteristic size of disorder, the roughness profiles are found to exhibit self-affine scaling features evolving with porosity. Starting from a null value and increasing the porosity, we quantitatively modify the self-affine scaling properties from anisotropic (at low porosity) to isotropic (for porosity >10%).

I. INTRODUCTION

Fractography, i.e., the morphological characterization of postmortem fracture surfaces, is a classical tool used to identify the mechanisms and the damaging processes (fatigue, stress corrosion, cleavage, plastic cavitation, crazing, etc.) responsible for failure (see [START_REF] Hull | Fractography: Observing, Measuring and Interpreting Fracture Surface Topography[END_REF] for a recent textbook). Since the pioneer work of Mandelbrot et al. [START_REF] Mandelbrot | Fractal character of fracture surfaces of metals[END_REF], numerous studies have evidenced the existence of self-affine scaling invariance over a wide range of length scales (see [START_REF] Bouchaud | Scaling properties of cracks[END_REF][START_REF] Bonamy | Failure of heterogeneous materials: A dynamic phase transition?[END_REF] for reviews). Early measurements in various materials including metallic alloys, ceramics, and rocks [START_REF] Bouchaud | Fractal dimension of fractured surfaces: A universal value?[END_REF][START_REF] Måløy | Experimental measurements of the roughness of brittle cracks[END_REF][START_REF] Schmittbuhl | Field measurements of the roughness of fault surfaces[END_REF] reported values for the roughness exponent ζ close to 0.8, suggesting the existence of a universal value [START_REF] Bouchaud | Fractal dimension of fractured surfaces: A universal value?[END_REF][START_REF] Måløy | Experimental measurements of the roughness of brittle cracks[END_REF], independent of the propagation mode and material. These experimental observations yielded intense theoretical [START_REF] Hansen | Roughness of crack interfaces[END_REF][START_REF] Bouchaud | Models of fractal cracks[END_REF][START_REF] Schmittbuhl | Interfacial crack pinning: Effect of nonlocal interactions[END_REF][START_REF] Larralde | The shape of slowly growing cracks[END_REF][START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF][START_REF] Hansen | Origin of the universal roughness exponent of brittle fracture surfaces: Stress weighted percolation in the damage zone[END_REF][START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF] and numerical [START_REF] Termonia | Formation of fractal cracks in a kinetic fracture model[END_REF][START_REF] Raisanen | Fracture of three-dimensional fuse networks with quenched disorder[END_REF][START_REF] Kumar | Crack surface roughness in three-dimensional random fuse networks[END_REF][START_REF] Nukala | Fracture roughness in three-dimensional beam lattice systems[END_REF] research.

Recently, the picture has been made more complex. First, the roughness exponent has been found to depend (slightly) on the anisotropy of the material microstructure [START_REF] Morel | Roughness analysis of the cracked surfaces of a face centered cubic alloy[END_REF] and on the fracture speed [START_REF] Mallick | Discrepancy between subcritical and fast rupture roughness: A cumulant analysis[END_REF]. Second, fracture surfaces were shown to exhibit anisotropic scaling features, characterized by two roughness exponents whether observed along the direction of the crack front or along the direction of crack growth [START_REF] Ponson | Two-dimensional scaling properties of experimental fracture surfaces[END_REF][START_REF] Ponson | Anisotropic self-affine properties of experimental fracture surfaces[END_REF]. Third, fracture surfaces exhibit anomalous scaling [START_REF] Lopez | Anomalous scaling of fracture surfaces[END_REF][START_REF] Morel | Anomalous roughening of wood fractured surfaces[END_REF]; the introduction of an additional global roughness exponent ζ glob is then necessary to describe the scaling between the global crack width and the specimen size. Fourth, multiaffinity, disappearing at large scales, was invoked [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF][START_REF] Santucci | Statistics of fracture surfaces[END_REF]. Finally, the scale invariance properties were found to depend on the propagation mode: The seminal self-affine feature with ζ 0.8 is to be linked with quasibrittle fracture; lower values for ζ , around 0.4-0.5 [START_REF] Plouraboué | Experimental study of fracture surface roughness on rocks with crack velocity[END_REF][START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF][START_REF] Bouchbinder | Fracture surfaces as multiscaling graphs[END_REF][START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone[END_REF] or even logarithmic roughness [START_REF] Dalmas | Crack propagation through phase separated glasses: Effect of the characteristic size of disorder[END_REF], are to be associated with brittle fracture; and a multiaffine regime with a roughness exponent close to 0.5 is characteristic of ductile fracture [START_REF] Bouchaud | Fracture through cavitation in a metallic glass[END_REF][START_REF] Ponson | Statistics of ductile fracture surfaces: The effect of material parameters[END_REF][START_REF] Srivastava | Effect of inclusion density on ductile fracture toughness and roughness[END_REF]. Note that different regimes can be observed on the same fracture surface [START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF][START_REF] Morel | Transient damage spreading and anomalous scaling in mortar crack surfaces[END_REF][START_REF] Gjerden | Universality classes in constrained crack growth[END_REF], depending on the scale of observation and its position with respect to the relevant sizes associated with the various dissipation mechanisms (plastic zone size, damage zone size, fracture process zone size). This enables us to infer the fracture process zone from a statistical characterization of the crack roughness [START_REF] Ponson | Cleaved surface of i-alpdmn quasicrystals: Influence of the local temperature elevation at the crack tip on the fracture surface roughness[END_REF].

Most past studies dedicated to the roughness of cracks aimed at characterizing the scale-invariant properties of fracture surfaces. This quest for universality classes, i.e., features independent of the fracturing conditions (loading and material parameters), strayed from the initial metrologic purpose of the fractography science. Here, we go back to this primary purpose and seek to characterize how the microstructure affects the fracture surface roughness. Rather than using existing materials like ceramics or rocks, we use homemade porous solids obtained by sintering spherical monodisperse grains. The advantage is that both the grain size and the porosity can be easily adjusted by modifying, respectively, the bead diameter and the sintering parameters, while presenting a structure of cemented grains very comparable to the structures found in natural heterogeneous brittle solids like rocks (sandstones, for instance) and concrete and other artificial heterogeneous brittle solids like pharmaceutical pills and sintered ceramics. The procedure for preparing the samples is described in detail in Sec. II A. The samples are finally broken using a technique detailled in Sec. II B that permits us to grow stable mode I cracks. The fractography of the surfaces (see Sec. II C) is used to determine the propagation mode of the fracture, i.e., interor intragranular propagation. The spatial correlation of the roughness of the crack surfaces is, finally, characterized by the structure function of profiles taken along and perpendicular to the direction of crack propagation. The evolution of this function as a function of the bead size, the direction of measurements, and the porosity is reported in Sec.III. Finally, in Sec. IV we discuss the change in the self-affine exponents observed when the porosity is modified.

II. SAMPLE PREPARATION AND EXPERIMENTAL SETUP

The experiments reported here consist (i) in driving a stable tensile crack with a controlled velocity in a material obtained by sintering polystyrene beads and (ii) subsequently, in analyzing the spatial distribution of the crack roughness. The material fabrication, fracture setup, and topographical recording of the postmortem roughness profiles are presented thereafter.

A. Sintering protocol

Samples are obtained by sintering monodisperse polystyrene beads (Dynoseeds) of various diameters. Four steps are achieved: (i) a homemade mold, whose geometry is given in Fig. 1(a), is filled with beads and heated up to 90% of the temperature at glass transition to soften the beads; (ii) a slowly linearly increasing compressive force is then applied up to the prescribed sintering load F s for nearly 1 h to achieve sintering; (iii) an annealing stage is performed, in which the mold is unloaded, unscrewed, and loosened, keeping the temperature high enough to avoid thermal shocks; and (iv) the sample is finally slowly cooled down to ambient temperature. The whole sintering protocol is sketched with more details in Fig. 1(b), and the relevant temperatures are listed in Fig. 1(c) (see also [START_REF] Barés | Failure in brittle heterogeneous materials: Intermittency, crackling and seismicity[END_REF] for an extensive presentation). The bead diameter d 0 was varied between 21 and 583 μm. By changing F s between 0.1 and 8 kN (corresponding to a loading pressure ranging from 6 to 490 kPa), we are able to fabricate specimens of porosities between ∼0% and 19% (see Table I). Smaller values of F s lead to unusable friable samples. The porosity is measured as the ratio = 1ρ/ρ 0 of the density of the sample ρ and a reference density ρ 0 corresponding to the density of the void-free samples obtained with the highest sintering force. It is found that ρ 0 is slightly bead size dependent, which may be explained by the slightly different initial bead material.

B. Fracture wedge-splitting tests

Specimens [Fig. 1(a)] are parallelepipeds of dimension 140 × 125 × W mm 3 in the x × y × z direction, where the specimen width W depends on the material porosity and ranges from 15 mm (for the most compacted specimens) to 20 mm (for the most porous specimens). A rectangular notch is milled by cutting a 42 × 30 × W mm 3 rectangle from the sample at one of the two lateral faces. A 10-mm-long, 2-mm-thick groove is then introduced in the middle of the rectangular notch using a diamond saw, and a seed precrack is added at the groove tip by means of a razor blade. The latter operation prevents the propagation of a dynamic fracture and allows us to grow slow stable cracks. Thereafter, we define the axes x parallel to the direction of crack propagation, y perpendicular to the mean crack plane, and z parallel to the mean crack front. The sample is then placed between the two jaws of a loading machine (a homemade or a commercial electromechanical machine). A homemade metallic wedge (semiangle of 15 • ) is placed between the jaw and the specimen notch [Fig. 2(a)]. It has been designed to minimize friction dissipation by means of two rollers and to convert the compression along x into a tension along y [START_REF] Scheibert | Brittle-quasibrittle transition in dynamic fracture: An energetic signature[END_REF][START_REF] Guerra | Understanding fast macroscale fracture from microcrack post mortem patterns[END_REF]. The sample is then compressed by lowering the upper jaw at a constant velocity. In the experiments reported here, the wedge speed V wdg varied from 1.6 to 1600 nm s -1 . During each test, the force F (t) applied by the wedge to the specimen was monitored in real time by a static load cell. As soon as the wedge starts to push on the specimen, F increases. When F (t) reaches the critical loading [Fig. 2(b)], the crack starts to propagate. Prior to fracture the load-displacement curve is linear, confirming the elastic behavior of the material. The crack velocity is observed first to increase (over ∼20 mm), then to stabilize at a roughly constant value, ∼100 times the wedge speed V wdg (over 40 to 60 mm), and, finally, to decrease over the last 20 mm. Our study focuses on the interval in which the velocity is constant and all measurements are done in that region.

C. Postmortem analysis of crack surfaces

Once the samples were broken into two halves, the morphology of the postmortem crack surfaces was analyzed. First, the surfaces were visualized by means of a Keyence numerical microscope. This microscope automatically performs a vertical scan at different heights and reconstructs a well-focused image pixel by pixel over the whole topography range, irrespective of the level of height variations.

Topography profiles are also recorded by a Bruker mechanical profilometer. Three profiles are scanned along the x direction and three others along the z direction. Their lengths are respectively 55 mm along x and 14 mm along z. Each profile has been positioned in the central part of the fracture surface away from the sample borders. They are located between x = 40 mm and x = 95 mm and at z = {-3,0,3} mm for profiles in the x direction and between z = -7 and z = 7 mm and at x = {-40,60,80} mm for profiles in the z direction. The out-of-plane and in-plane resolutions are δy = 0.1 nm and δx = δz = 0.8 μm, which corresponds, per profile, to 68 750 points in the x direction and 16 250 in the z direction.

For d 0 = 583 μm, the fracture surfaces have been measured using another profilometer due to the vertical displacement limitation of the Bruker profilometer [START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone[END_REF]. Its out-of-plane and in-plane resolutions are δy = 100 nm and δx = δz = 1 μm.

To have a better statistic, ten 75-mm-long profiles have been measured in the x direction. Yet, in this case, due to the small number of beads in the thickness, no measurements were made in the z direction.

III. RESULTS

A general overview of the fracture surfaces obtained by changing the bead size, the porosity, and the crack velocity is first given. In most cases, the propagation is intergranular. We focus on this case, investigating the effects of the microstructure size and porosity on the fracture surfaces. Figures 3(a)-3(c) correspond to a null porosity and to samples broken at the same wedge velocity but with different bead diameters. The particle surfaces are observed to be completely deformed, resulting in sharp edges and corners, in agreement with the zero value of the porosity. Also, the surfaces are found to be made of smooth facets. Such a faceted morphology is characteristic of a brittle intergranular fracture, with a crack growing along the grain-grain interfaces. Note that the facets are slightly more elongated along x than along z [visible in Fig. 3(b)]. This anisotropy is due to the sintering process. In fact, the compression was applied along z, or in a granular packing, the two components σ xx and σ zz of the stress tensor are classically related [START_REF] Janssen | Versuche ber getreidedruck in silozellen[END_REF] via σ xx = Kσ zz , where the Janssen constant K is less than 1, typically around 0.6-0.8. During sintering, the beads, hence, contract more along z than along x. Figure 3(d), corresponding to a lower wedge speed, V wdg = 1.6 nm/s, presents a qualitatively different morphology: The facets are observed to be blurred, with a multitude of small-scale fragments. They betray an intragranular propagation mode, with a crack propagating throughout the grains, without necessarily following the graingrain interfaces. Figures 3(e) and 3(f) show microscopic visualizations of fracture surfaces of sintered samples realized with a lower sintering force F s . At low F s [Fig. 3(e)], the beads globally maintain a spheric shape. The contact areas between adjacent beads are small and disk shaped. As a result, an important porosity is found in such samples. When F s gets larger, the contact surfaces grow and intersect, resulting in the appearance of sharp edges of increasing length [Fig. 3(f)], while the undeformed spherical parts of the particles shrink, so that the volume of the pores, delimited by these parts, decreases. When F s is large enough to close all the pores [Fig. 3(c)], the undeformed parts of the particles disappear and the corners are sharp.

A. Fractography

Figure 4 presents typical topographical profiles H (z) as measured along the z direction on the fracture surfaces shown in Fig. 3. In Figs. 4(a)-4(c), associated with intergranular fracture, the corrugations in the y and z directions are of the same order, that of the bead size. Conversely, the roughness observed for the intragranular mode [Fig. 4(d)] appears much flatter and exhibits plateaus much longer than the bead diameter. These plateaus correspond to zones where the crack has cut throughout the grains, without being perturbed by the interfaces. Finally, Figs. 4(e) and 4(f) show typical profiles at finite porosity. Note the increase in roughness with amplitudes higher than the bead diameter d 0 . The shift from inter-to intragranular fractures occurs for a fracture velocity of 1.6 nm/s. In the present study, we focus on the intergranular case, hence on experiments corresponding to V wdg > 1.6 nm/s. Works on the transition between intra-and intergranular propagation and on intragranular fracturation are in progress.

B. Effect of the microstructure length scale

We first restrict the analysis to the effect of the microstructure length scale and, hence, only consider specimens (nos. 1-17 in Table I) with a very low porosity (<0.5%). To quantify the spatial distribution of the crack roughness, we computed the structure functions [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF] S z ( z) along z and S x ( x) along x, defined by

S z ( z) = (H (z + z) -H (z)) 2 ,
(1) where the operator refers to an average over all possible positions z and x, respectively.

S x ( x) = (H (x + x) -H (x)) 2 , ( 2)
Figure 5 displays in a log-log plot the structure functions S x measured along x for different bead diameters. Two scaling regimes can be distinguished. At large scales ( x 10 μm), S x slowly increases with x, with a prefactor increasing with d 0 . This betrays the fact that the roughness scales with d 0 . At small scales ( x 10 μm), all curves collapse, and the roughness in this regime is weakly dependent on d 0 .

The diameter d 0 appears to be the only length scale in the problem. Therefore, we made dimensionless the structure functions by making x → x/d 0 and S x → S x /d 2 0 (Fig. 5, inset). A good collapse is obtained, confirming that d 0 is the length scale governing the spatial distribution of the roughness.

Nevertheless, the collapse shown in the inset in Fig. 5 is not perfect and a clear scattering remains visible. We applied the global minimization technique described in Ref. [START_REF] Cambonie | Experimental study of counter-rotating vortex pair trajectories induced by a round jet in cross-flow at low velocity ratios[END_REF] to reduce the dispersion of the data. For each curve, an optimal length is obtained. The values are listed in Table II for each sample. Figures 6(a) and 6(b) show the structure functions in the x and z directions normalized by x and z . As shown, the sole parameters x and z are sufficient to achieve a very good collapse. We checked that these lengths correspond to the microstructure sizes in both directions (as, e.g., shown in Fig. 3). In other words, the slight differences observed among x , z , and d 0 result from the sintering process. Moreover, two power-law regimes, characterized by two different exponents, are observed. To determine these exponents, we performed linear fits of the data displayed in Fig. 6. The roughness exponents ζ +,- x and ζ +,- z , for length scales above and below d 0 , are defined by the relations S x ∝ x 2ζ x and S z ∝ z 2ζ z . Table III lists the roughness exponent values along the propagation and front crack directions. In the small-scale regime, the exponents ζ -are close to 0.8 in both directions. At larger scales, the exponents ζ + are functions of the direction of measurement, the roughness along the crack being more important than that in the propagation direction. Altough the exponents are relatively small with respect to usual TABLE II. Optimal lengths x and z along the x and z directions used to collapse (in Fig. 6) the structure functions. III.

values of power-law behaviors, the insets in Fig. 6 clearly prove that they are finite and exclude a logarithmic behavior.

C. Effect of the porosity

We now investigate the influence of the porosity . Figure 7 shows the dimensionless structure functions. Contrary to the previous case, the data do not collapse. First, the amplitude of the structure functions is found to increase with at all scales. This is compatible with the observations of rougher profiles as the porosity increases (Fig. 4). At small length scales ( x/d 0 < 1), the exponents are similar for all porosities. A linear regression of the data gives a roughness exponent ζ - x = 0.88 ± 0.04, very close to that otained for = 0. Yet, for x/d 0 > 1, the exponents increase with the porosity.

In Fig. 8, we analyze more quantitatively the effect of on the roughness. The roughness amplitudes, defined as the value of the structure functions at x/d 0 = 1, are observed to increase linearly with [Fig. 8(a)]. Figure 8(b) shows the evolution of the large-scale roughness exponents ζ +

x and ζ + z as a function of . As expected, when tends toward 0, the values obtained in the no-porosity case (Table III) are recovered. They increase with , up to values of ζ + z ∼ ζ + x ∼ 0.48 for = 20%. No significant anisotropy is observed in the roughness amplitude [Fig. 8(a)]. Conversely, the scaling anisotropy, i.e., the difference between ζ +

x and ζ + z , observed for = 0 is recovered. However, it decreases as increases.

To quantify the scaling anisotropy evolution, we plot the ratio ζ + z /ζ + x as a function of in Fig. 9. It decreases from ζ + z /ζ + x = 1.3 ± 0.05 to 1 as goes from 0 to 10%. Above, isotropic scaling is obtained.

IV. DISCUSSION

A. Comparison with other experimental results

The roughness value ζ -measured at small scales is observed to be independent of the bead size, the direction measurement, and the porosity. It is compatible with the ζ = 0.8 value reported in the literature for a variety of materials [START_REF] Bouchaud | Scaling properties of cracks[END_REF][START_REF] Bouchaud | Fractal dimension of fractured surfaces: A universal value?[END_REF][START_REF] Måløy | Experimental measurements of the roughness of brittle cracks[END_REF]. Here, ζ -is observed at a scale smaller than the bead size. This small-scale regime is not, strictly speaking, a self-affine regime: The ζ -value can be interpreted as the measured roughness exponent of a piecewise linear profile analyzed at scales smaller than that of its segments (see [START_REF] Lechenault | Effects of finite probe size on self-affine roughness measurements[END_REF] for a related analysis).

Regarding the roughness properties at large scales, it is interesting to compare our results with what has been reported for ceramics obtained by sintering oxide glass beads [START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF]. The roughness amplitude [defined as √ S(δr = d 0 )] is found to increase linearly with as observed in Ref. [START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF] for glass ceramics. Conversely, our data converge toward a finite amplitude as → 0, while in glass ceramics the amplitude was vanishing. The main difference is that, in Ref. [START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF], decreasing induces a transition from an inter-to an intragranular propagation mode, which is not the case here.

The values of the roughness exponents are also, within the error bars, similar to the values of ζ 0.4-0.5 reported in these glass ceramics [START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF] and in sandstone [START_REF] Plouraboué | Experimental study of fracture surface roughness on rocks with crack velocity[END_REF][START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone[END_REF]. Our material is hence more representative of conventional heterogeneous brittle materials such as ceramics and some rocks. As a bonus, the smaller error bars and the larger range of explored have evidenced the increase in ζ with . They have also revealed scaling anisotropy in the absence of any porosity, which decreases as increases and disappears for 10%.

B. Comparison with theoretical results

Since the early 1990s, a great deal of theoretical work has focused on the spatial distribution of the fracture surface roughness in brittle heterogeneous materials. They can be classified into two categories: (i) elastic string models, which consider the crack front as an elastic line propagating through randomly distributed obstacles [START_REF] Bouchaud | Models of fractal cracks[END_REF][START_REF] Schmittbuhl | Interfacial crack pinning: Effect of nonlocal interactions[END_REF][START_REF] Larralde | The shape of slowly growing cracks[END_REF][START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF][START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF][START_REF] Katzav | Fracture surfaces of heterogeneous materials: A 2D solvable model[END_REF]; and (ii) random lattice models, which model the material by a network of fuses, springs, or beams with randomly distributed breakdown thresholds [START_REF] Hansen | Roughness of crack interfaces[END_REF][START_REF] Raisanen | Fracture of three-dimensional fuse networks with quenched disorder[END_REF][START_REF] Nukala | Fracture roughness in three-dimensional beam lattice systems[END_REF][START_REF] De Arcangelis | Scaling laws in fracture[END_REF]. The main difference between the predictions of these two classes of models is that models of type i naturally lead to anisotropic surfaces, where the direction of front propagation plays a specific role, while models of type ii lead to isotropic surfaces [START_REF] Ponson | Two-dimensional scaling properties of experimental fracture surfaces[END_REF].

In essence, elastic line models address situations of nominally brittle fracture, with nonporous materials. Hence, they are the relevant theoretical framework in which to discuss the observations reported in Sec. III B. In particular, elastic string approaches lead to (a) the length scale of the microstructural disorder being the single relevant length scale for the structure function and (b) anisotropic surfaces, where the direction of front propagation plays a specific role.

Both these predictions are in agreement with the experimental observations (Figs. 6 and8).

Conversely, the values of the roughness exponents, ζ + x ≈ 0.27 and ζ + z ≈ 0.35, observed here are clearly distinct from the values predicted theoretically. In particular, the most refined models [START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF][START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF] attempted to derive rigorously the equation of motion of the elastic string from linear elastic fracture mechanics. They yield either logarithmic scaling [START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF] or ζ x = 0.50 ± 0.05 and ζ z ≈ 0.385 ± 0.05 [START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF], according to the disorder introduced in the equation. The values reported here are clearly different, which indicates that some key ingredients are missing in the models. In particular, no model takes into account the T -stress influence. The T stress here has been determined [START_REF] Hattali | Evaluation using digital image correlation and finite element method of stress intensity factors and t-stress in wedge splitting test[END_REF] both by finite-element simulations and digital image correlations [START_REF] Mathieu | Identification of a crack propagation law by digital image correlation[END_REF]. As shown in Fig. 10, the T stress is positive, so that the crack path should be unstable toward any perturbation of its rectilinear path [START_REF] Cotterell | Slightly curved or kinked cracks[END_REF]. In other words, additional correlations are anticipated to result from a positive T stress, and hence a larger value for the roughness exponent measured along x. If this scenario is correct, it may partly explain the departure from logarithmic behavior.

As mentioned above, in elastic line models, the dimensionless curves S x /d 2 0 vs x/d 0 and S z /d 2 0 vs z/d 0 are anticipated to be universal. But we observed a scaling function of (see Fig. 7). This indicates that elastic line approaches stop being relevant as the porosity becomes finite. This is thought to be due to the presence of microcracks forming at the pores, causing the fracture propagation mode to shift gradually from nominally brittle to quasibrittle. The anisotropy decrease with is consistent with this scenario: As microcracking develops, the fracture surface stops being the trace left by a single propagating line in an otherwise intact material but results from the coalescence of multiple microcracks instead. Random lattice models include such processes and may be the relevant framework in which to describe this quasibrittle regime. In this context, it is interesting to note that the roughness exponents ζ + x = ζ + z ≈ 0.45 measured here on isotropic fracture surfaces, i.e., above = 10%, are consistent with the values ζ loc = 0.42, ζ loc = 0.5, and ζ loc = 0.48 observed, respectively, in three-dimensional random fuse [START_REF] Raisanen | Fracture of three-dimensional fuse networks with quenched disorder[END_REF], spring [START_REF] Parisi | Roughness of fracture surfaces[END_REF], and beam [START_REF] Nukala | Fracture roughness in three-dimensional beam lattice systems[END_REF] networks.

In the above scenario, porosity provides a tunable parameter to go gradually from nominally brittle fracture, where the disorder effect is distortion of the front propagation, to quasibrittle fracture, where the fracture surfaces emerge from a percolating path throughout the microcrack cloud. At present, the theoretical descriptions of these two situations belong to two distinct realms: elastic line models for the former and random lattice models for the latter. Unifying these two frameworks represents an important challenge for future investigations (see [START_REF] Gjerden | Universality classes in constrained crack growth[END_REF] for a recent thrust in this direction).

V. CONCLUSION

In this paper, we have used homemade cemented grain materials to investigate the influence of the microstructure size and porosity on the fracture surfaces. The following observations have been made:

(a) The roughness displays scale-invariant morphological features that depend on both the grain size and the porosity.

(b) At zero porosity, the structure functions measured along a direction and normalized by the grain size collapse into a single master curve.

(c) Scaling anisotropy is well pronounced for = 0 but decreases with and vanishes for > 10%.

( These observations are compatible with previous studies considering other sintered glass bead materials [START_REF] Ponson | Low self-affine exponents of fractured glass ceramics surfaces[END_REF] or sandstone [START_REF] Plouraboué | Experimental study of fracture surface roughness on rocks with crack velocity[END_REF][START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone[END_REF]. But in these previous studies, the precision was not sufficient to detect any clear variation with . However, the values of the roughness exponents, ζ +

x ≈ 0.27 and ζ + z ≈ 0.35, observed at low porosity are clearly distinct from the values predicted theoretically: Additional investigations are necessary to take into account some missing ingredients like the T stress. Conversely, for > 10%, the roughness exponents ζ + x = ζ + z ≈ 0.45 are consistent with the values of ζ loc between 0.4 and 0.5 obtained with random networks [START_REF] Raisanen | Fracture of three-dimensional fuse networks with quenched disorder[END_REF][START_REF] Nukala | Fracture roughness in three-dimensional beam lattice systems[END_REF][START_REF] Parisi | Roughness of fracture surfaces[END_REF]. Modulating has, thus, permitted us to modify quantitatively the value of the roughness exponents. Note that, in conventional interface growth problems, the roughness exponent value is characteristic of a universality class. Their continuous evolution with observed herein may pose a rather severe test for current and future competing models of heterogeneous fracture. The T -stress component acting on the fracture is controlled by the macroscopic shape of the porous samples; future wors will consider the effect of the T stress on the roughness.

FIG. 1 .

 1 FIG. 1. (Color online) (a) Geometry and dimensions of the wedge-splitting samples. (b) Evolution sketch of the temperature of the oven, of the sample, and of the force applied during the sintering protocol (see [38] for details). (c) Temperature values at the different stages of the sintering protocol for the different bead sizes.

FIG. 2 .

 2 FIG. 2. (Color online) (a) Sketch of the experimental setup. (b) Typical force-vs-displacement curve.

Figure 3

 3 Figure 3 presents microscope visualizations of typical fracture surfaces, for different bead diameters, porosities, and wedge velocities. The crack propagates from left to right.Figures3(a)-3(c) correspond to a null porosity and to samples broken at the same wedge velocity but with different bead diameters. The particle surfaces are observed to be completely deformed, resulting in sharp edges and corners, in agreement with the zero value of the porosity. Also, the surfaces are found to be made of smooth facets. Such a faceted morphology is characteristic of a brittle intergranular fracture, with a crack growing along the grain-grain interfaces. Note that the facets are slightly more elongated along x than along z [visible in Fig.3(b)]. This anisotropy is due to the sintering process. In fact, the compression was applied along z, or in a granular packing, the two components σ xx and σ zz of the stress tensor are classically related[START_REF] Janssen | Versuche ber getreidedruck in silozellen[END_REF] via σ xx = Kσ zz , where the Janssen constant K is less than 1, typically around 0.6-0.8. During sintering, the beads, hence, contract more along z than along x. Figure3(d), corresponding to a lower wedge speed, V wdg = 1.6 nm/s, presents a qualitatively different morphology: The facets are observed to be blurred, with a multitude of small-scale fragments. They betray an intragranular propagation mode, with a crack propagating throughout the grains, without necessarily following the graingrain interfaces.

FIG. 3 .

 3 FIG. 3. (Color online) Microscope visualization of wedge-split fractured surfaces of the sintered-PS-bead samples. The bead diameter d 0 , wedge speed V wdg , sintering load F s , and corresponding porosity are as follows. (a) d 0 = 21 μm, V wdg = 16 nm s -1 , F s = 8 t, = 0%. (b) d 0 = 228 μm, V wdg = 16 nm s -1 , F s = 8 t, = 0%. (c) d 0 = 583 μm, V wdg = 16 nm s -1 , F s = 8 t, = 0%. (d) d 0 = 583 μm, V wdg = 1.6 nm s -1 , F s = 8 t, = 0%. (e) d 0 = 42 μm, V wdg = 16 nm s -1 , F s = 0.1 t, = 15%. (f) d 0 = 42 μm, V wdg = 16 nm s -1 , F s = 1 t, = 1.45%.

FIG. 4 .

 4 FIG. 4. Experimental height profiles along the crack propagation direction of the same samples as in Fig. 3. At the right, the bead shape is given as a reference, using the same abscissa and ordinate scales as in the plots. Since these scales are different, the bead appears distorted. (a) Experiment 2: d 0 = 21 μm, V wdg = 16 nm s -1 , = 0%. (b) Experiment 13: d 0 = 228 μm, V wdg = 16 nm s -1 , = 0%. (c) Experiment 16: d 0 = 583 μm, V wdg = 16 nm s -1 , = 0%. (d) Experiment 14: d 0 = 583 μm, V wdg = 1.6 nm s -1 , = 0%. (e) Experiment 21: d 0 = 42 μm, V wdg = 16 nm s -1 , = 15%. (f) Experiment 3: d 0 = 42 μm, V wdg = 16 nm s -1 , = 1%

FIG. 5 .

 5 FIG. 5. Structure functions S x of intergranular experiments, along the propagation direction at = 0%. Each symbol corresponds to a given value of d 0 . Inset: Same distributions made dimensionless by d 2 0 .

FIG. 6 .

 6 FIG. 6. (Color online) Structure functions normalized by the optimal lengths x and z . Measurements performed along (a) x (propagation direction) and (b) z (front crack). Both axes are logarithmic in the main panels and semilogarithmic in the insets. Straight lines in the main panels are power-law fits with exponents reported in TableIII.

Fig. 6 .

 6 FIG. 7. Dimensionless structure functions along the propagation direction for different porosities .

FIG. 8 .

 8 FIG. 8. (a) Normalized roughness amplitudes in the x and z directions √ S(d 0 )/d 0 versus porosity . (b) Large-scale roughness exponents ζ + vs . Symbols match the experiments in Fig. 7. Crosses correspond to the averaged exponent values in TableIII.

FIG. 10

 10 FIG. 10. (Color online) T stress as a function of the crack tip position x tip obtained by classical finite-element simulations and digital image correlation.

TABLE I .

 I Values of experimental parameters associated with the different fracture surfaces analyzed here: bead diameter d 0 prior to sintering, sintering force F s , wedge velocity V wdg , and porosity .

				Experiment No.
		1	2	3	4	5	6
	d 0 (μm)	21	21	42	42	42	81
	F S ( k N )	4	8	1	8	8	2
	V wdg (nm s -1 )	16	16	16	16	16	16
	( % )	1	0	1	0	0	1
				Experiment No.
		7	8	9	1 0	1 1	1 2
	d 0 (μm)	81	81	81	81	228	228
	F S ( k N )	4	8	8	8	2	4
	V wdg (nm s -1 )	16	16	16	16	16	16
	( % )	1	0	0	0	1	0
				Experiment No.
		13	14	15	16	17
	d 0 (μm)	228	583	583	583	583
	F s ( k N )	8	8	8	8	8
	V wdg (nm s -1 )	16	1.6	16	160	1600
	( % )	0	0	0	0	0
				Experiment No.
		18	19	20	21	22	23
	d 0 (μm)	21	81	42	42	42	42
	F S (kN)	0.1	0.1	0.1	0.2	0.6	0.8
	V wdg (nm s -1 )	16	16	16	16	16	16
	(%)	11.65	18.7	15.1	10.4	≈3	≈2

TABLE III .

 III Roughness exponents obtained by linear fits performed on the curves in

  d) The values of the roughness exponents increase with , from ζ + x ∼ 0.27 and ζ + z ∼ 0.35 at = 0% to ζ + x ∼ ζ + z ∼ 0.48 at = 20%.

ACKNOWLEDGMENTS

We acknowledge support by the ANR Mephystar (Grant No. ANR-09-SYSC-006-01). We acknowledge P. Viel for the use of the Bruker profilometer and T. Bernard for technical support. We thank R. Pidoux, L. Auffray, and T. Bernard for realizing the samples.


Antoine Campagne

Basile Gallet

Frederic Moisy

Pierre-Philippe Cortet

Frédéric Moisy

Disentangling inertial waves from eddy turbulence in a forced rotating-turbulence experiment

Keywords: number(s): 47.27.-i, 47.32.Ef, 47.35.-i

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

The energy content of turbulence is usually characterized by the energy distribution among spatial scales, either in physical or in Fourier space. For rotating, stratified, or magnetohydrodynamic turbulence [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF], waves can propagate and coexist with classical eddies and coherent structures, which advocates for a spatiotemporal description of such flows. While temporal fluctuations are usually slaved to the spatial ones via sweeping effects in classical turbulence [START_REF] Chen | [END_REF]3], they are expected to be governed by the dispersion relation of the waves for time scales much smaller than the eddy turnover time. The latter regime is the subject of wave-turbulence theory, in which the assumption of weak nonlinear coupling between waves allows one to predict scaling laws for the spatial energy spectrum [START_REF] Zakharov | Wave Turbulence[END_REF][START_REF] Nazarenko | Wave Turbulence[END_REF].

It is a matter of debate whether wave-turbulence theory (also known as weak-turbulence theory) is a good candidate to describe rotating turbulence in the rapidly rotating limit. Solutions to the linearized rotating Euler equation can be decomposed into inertial waves, which satisfy the anisotropic dispersion relation

σ (k) = 2 |k | |k| , (1)
where is the rotation rate and k the component of the wave vector k along the rotation axis (referred to as the vertical axis by convention) [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF]. Accordingly, only fluid motions at frequencies σ smaller than the Coriolis frequency 2 correspond to wave propagation. Fluid motions of weak amplitude and slowly varying in time (σ 2 ) can be described in terms of waves with nearly horizontal wave vectors: They tend to be two-dimensional three-component (2D3C), invariant along the rotation axis, a result known as the Taylor-Proudman theorem.

The trend towards two-dimensionality is a landmark in rotating turbulence, observed in both experiments and numerical simulations [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF][START_REF] Hopfinger | [END_REF][8][9][10][11]. It originates from the modification of the nonlinear interactions by the Coriolis force, which yields preferential energy transfers towards modes with almost horizontal wave vectors [12][13][14][15]. In the frequency domain, this corresponds to the generation of slow dynamics compared to the characteristic frequency at which energy is supplied to the system. These anisotropic energy transfers can be accounted for in terms of resonant and nearly resonant triadic interactions of inertial waves [16][17][18][19]. A major feature of rotating turbulence is the emergence of inverse energy transfers in the horizontal plane [13,[19][20][21][22][23][24][25][26]. Inverse transfers between 3D fast wave modes and the 2D slow vortex mode, mediated by near-resonant triadic interactions, are allowed at finite Rossby number only [22,23,[27][28][29]. The 2D mode is therefore fed either from the coupling with the 3D modes at finite Rossby number or from direct energy input by the forcing. One naturally expects the energy within this 2D mode to undergo an inverse energy cascade, similar to that of classical (nonrotating) 2D turbulence [30,31].

Such coexistence between 2D and 3D flows is relevant to most experiments and numerical simulations and cannot be accounted for by wave-turbulence theory, which describes the direct energy cascade arising from resonant triadic interactions of 3D wave modes only [32,33]. This theory therefore provides only a partial description of rotating turbulence in realistic systems and careful experimental and numerical studies remain necessary to assess its range of validity.

Laboratory experiments differ from most numerical and theoretical studies by the presence of rigid horizontal boundaries, where the rotating flow achieves no-slip conditions through Ekman layers [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. In a laminar Ekman layer, the balance between the viscous and Coriolis forces leads to a boundary layer thickness δ Ek √ ν/ . The belief is that, provided the experimental tank is deep enough, the bulk turbulent flow away from the top and bottom boundaries should resemble the one obtained in the ideal 3D periodic or stress-free domains considered in most numerical and theoretical studies. Closer to the horizontal walls, the boundary layers induce Ekman friction that is not taken into account by most numerical studies.

In the laboratory, the energy dissipation of rotating turbulent flows originates from three main contributions: bulk viscous dissipation of 3D flow structures, bulk viscous dissipation of quasi-2D flow structures (somewhat similar to the bulk energy dissipation of 2D turbulence), and dissipation through Ekman friction on the horizontal boundaries.

In spite of the importance of the 2D mode in most geophysical and laboratory flows, the 3D fluctuations still play a crucial role in the dynamics of rotating turbulence, because they are much more efficient at dissipating energy. This key feature is illustrated in Fig. 2, from data obtained in the present experiment (setup sketched in Fig. 1; see Sec. III for details): We decompose the turbulent velocity field into a vertically averaged 2D flow and a vertically dependent 3D remainder and show the corresponding energies and energy dissipation rates as a function of global rotation. For maximum rotation, although the 3D component contains a small fraction of the total kinetic energy, its dissipation rate is as large as that of the vertically invariant 2D component. Moreover, both of these dissipations are larger than an estimate of the frictional losses due to laminar Ekman layers (see Sec. III). An accurate description of the 3D structures of the flow is therefore essential to characterize the energy fluxes in rotating turbulence at moderate Rossby number.

A primary goal in this direction is to determine the range of scales and frequencies for which 3D fluctuations follow the inertial-wave dispersion relation. This requires a full spatiotemporal analysis, which is very demanding in general for wave-turbulence systems: The accessible range of scales is usually limited in experiments, whereas long integration times are prohibitive in numerical simulations. The case of rotating turbulence is particularly delicate because of the specific form of the dispersion relation (1): The frequency is not related to the wave number, as in conventional isotropic wave systems such as surface waves [START_REF] Berhanu | [END_REF] or elastic waves [36,37], but to the wave-vector orientation only.

The recent studies of Clark di Leoni et al. [38] and Yarom and Sharon [39] constitute important steps forward in this respect. Using numerical simulation of rotating turbulence forced at large scale, Clark di Leoni et al. [38] observe a clear concentration of energy along the dispersion relation of inertial waves and provide a detailed analysis of the various time scales of the system. They observe a wave-dominated regime at large scale and a sweeping-dominated regime at small scale (see also Ref. [40]). In the experiment of Yarom and Sharon [39] the forcing consists of a random set of sources and sinks at the bottom of a rotating water tank. They measure 3D2C velocity fields using a scanning particle image velocimetry (PIV) technique and observe also good agreement with the inertial-wave dispersion relation. In both Refs. [38] and [39], the inertial waves are observed at scales smaller than the injection scale, suggesting that they are fed by forward energy transfers, which is consistent with the predictions of wave-turbulence theory.

The aim of the present paper is to further analyze experimentally the range of spatiotemporal scales at which inertial waves can be detected in rotating turbulence. Stationary rotating turbulence is produced by a set of vortex dipole generators that continuously inject turbulent fluctuations towards the center of a rotating water tank where measurements are A laser sheet illuminates a vertical slice through a horizontal glass lid covering the fluid. Two-dimensional three-component velocity measurements are performed using stereoscopic particle image velocimetry in a vertical square domain of size x × z = 14 × 14 cm 2 , shown as a dashed square in (a). performed. We showed in Ref. [26] that this configuration generates a double energy cascade at large rotation rate: an inverse cascade of horizontal energy and a direct cascade of vertical energy, which behaves approximately as a passive scalar advected by the horizontal flow. Here we perform a detailed spatiotemporal analysis using two-point spatial correlations of the temporal Fourier modes computed from time-resolved 2D3C velocity fields measured by stereoscopic PIV in a vertical plane. We observe that, at large scales and frequencies, the spatiotemporal anisotropy of the energy distribution is well described by the dispersion relation of inertial waves, whereas smaller scales are dominated by the sweeping of the waves by the energetic large-scale flow.

II. EXPERIMENTAL SETUP

The experimental setup, sketched in Fig. 1, is similar to the one described in Refs. [26,41] and is only briefly described here. It consists of a 125 × 125 × 65 cm 3 glass tank, filled with H = 50 cm of water and mounted on a 2-m-diam rotating platform that rotates at a rate in the range 0.21-1.68 rad s -1 (2-16 rpm). Turbulence is produced in the rotating frame by a set of ten vertical vortex dipole generators organized as a circular arena of 85 cm diameter around the center of the water tank. This forcing device was initially designed to generate turbulence in stratified fluids and is described in detail in Refs. [42,43]. Each generator consists of a pair of vertical flaps, 60 cm high and L f = 10 cm long, alternatively closing rapidly and opening slowly in a cyclic motion of period T 0 = 2π/σ 0 = 8.5 s. The closing stage is achieved with the flaps rotating at an angular velocity σ f = 0.092 rad s -1 and a random phase shift is applied between the generators.

In the laminar regime, a single pair of flaps generates vortex dipoles with core vorticity ω f . Additional PIV measurements close to a vortex dipole generator indicate that this core vorticity is governed by the vorticity in the viscous boundary layer of the flap ω f ∼ σ f L f /δ, where σ f L f is the flap velocity and δ the viscous boundary layer thickness. In the present experiment, the vortex dipoles are unstable and the closing of the flaps therefore produces small-scale 3D turbulent fluctuations that are advected towards the center of the arena by the remaining large-scale dipolar structure.

The turbulent Reynolds number, computed from the rms velocity and the horizontal integral scale, is about 400 in the center of the flow and the turbulent Rossby number covers the range 0.30-0.07 for = 2-16 rpm [26]. We measure the three components of the velocity field u = u x e x + u y e y + u z e z (with e z oriented vertically, along the rotation axis) in a vertical square domain of size x × z = 14 × 14 cm 2 located at the center of the circular arena at mid-depth, using a stereoscopic PIV system [START_REF] Davis | Anna-Vandenhoeck-Ring 19[END_REF][START_REF] Moisy | PIVMat toolbox for MATLAB[END_REF] embarked on the rotating platform. These 2D3C velocity fields are sampled on a grid of 80 × 80 vectors with a spatial resolution of 1.75 mm. Two acquisition sets are recorded for each rotation rate : one set of 10 000 fields at 0.35 Hz and one set of 1000 fields at 1.5 Hz. The combination of these two time series results in a temporal spectral range of three decades.

III. 2D VS 3D FLOW COMPONENTS

In the present experiment, energy is primarily injected in the 2D mode (vertically invariant), but the instabilities in the vicinity of the flaps rapidly feed 3D fluctuations that are advected in the central region. Energy transfers between the 2D and 3D flow components, which vanish in the weak-turbulence limit (Ro → 0), are allowed in our system because of the moderate value of the Rossby number. It is therefore of interest to quantify the energy contained in the 2D and 3D components of the flow. We estimate the vertically averaged 2D flow as the average of the velocity field over the vertical extent z of the PIV field, and the remaining z-dependent 3D flow as u 3D = uu 2D . We compute the energy per unit mass of these two flow components as u 2 2D /2 and u 2 3D /2, with the overline denoting the temporal average and angular brackets the spatial average over the PIV field. They are plotted in Fig. 2(a) as a function of the rotation rate . Because of the limited height of the PIV field, the 2D flow estimated from Eq. ( 2) unavoidably contains 3D fluctuations associated with vertical scales larger than z, so the measured 2D energy may overestimate the true one.

u 2D = 1 z z 0 u(x,z)dz, (2)
Figure 2(a) shows that without rotation the 2D and 3D components of the flow have comparable energy. With rotation, the 2D energy increases with , following approximately the power law 2/3 [41], whereas the 3D energy remains approximately constant and represents only 5% of the total energy at the largest rotation rate. Although most of the energy is contained in the 2D flow component for = 0, a significant fraction of the dissipation still arises from the 3D fluctuations. Assuming axisymmetry, we compute an estimate of the energy dissipation rate = ν (∂u i /∂x j ) 2 (

This dissipation rate, computed for both u 2D and u 3D , is shown in Fig. 2(b). Since the derivatives are obtained from finite differences at the smallest resolved scale, the computed dissipation underestimates the true one (the PIV resolution is 1.75 mm while the Kolmogorov scale is of order of 0.6 mm [41]). However, we expect the measured evolution of with to reflect the true one. We first compare these bulk energy dissipation rates to an estimate of frictional losses due to laminar Ekman layers Ek ν

U 2 ⊥rms δ 2 Ek δ Ek H = √ ν U 2 ⊥rms H
, where U ⊥rms is the root-mean-square horizontal velocity. This estimate ranges from 8 × 10 -9 m 2 s -3 for = 2 rpm to 1 × 10 -7 m 2 s -3 for = 16 rpm; it is smaller than the bulk energy dissipation of both the 2D and 3D parts of the turbulent flow, by a factor of 10 for slow rotation and 4 for rapid rotation. A detailed experimental characterization of these Ekman layers would be necessary to validate the assumption of laminar layers, but it is beyond the scope of the present study.

We now compare the bulk energy dissipation rates in the 2D and 3D parts of the turbulent flow. Remarkably, while the 3D fluctuations represent a small fraction of the total energy, they account for a large fraction of the dissipation at all rotation rates. It is therefore of interest to investigate these 3D modes and to determine to what extent they can be described in terms of inertial waves.

IV. TEMPORAL ANALYSIS

We now focus on the temporal dynamics of the velocity field, which we characterize through the energy distribution of turbulent fluctuations as a function of angular frequency σ . This temporal energy spectrum is defined as

E(σ ) = 4π T | ũi (x,σ )| 2 , ( 4
)
where

ũi (x,σ ) = 1 2π T 0 u i (x,t)e -iσ t dt (5)
is the temporal Fourier transform of the velocity field u i (x,t) (with i = x,y,z), T the acquisition duration, and the angular brackets the spatial average. The normalization is such that

u 2 i = ∞ 0 E(σ )
dσ , with the overline the temporal average. We use the standard Welch method [START_REF] Welch | [END_REF] to improve the statistical convergence of the power spectrum.

We plot E(σ ) for each rotation rate in Fig. 3. For = 0 we observe a global increase with of the energy at all frequencies, consistently with the behavior of the overall energy in Fig. 2(a). These spectra for = 0 strongly differ from the nonrotating spectrum, with relatively much more energy at low frequency in the rotating case: Global rotation induces slow dynamics. A step towards a description of the flow anisotropy in the frequency domain can be provided by further decomposing the power spectrum density (4) as This decomposition highlights the frequency-dependent componentality of turbulence, i.e., the distribution of energy among the different velocity components, which is related to the polarization anisotropy [12,15,16]. This is not to be confused with the frequency-dependent dimensionality of turbulence, which compares the vertical and horizontal characteristic scales at a given frequency (described in Sec. V).

E(σ ) = E (σ ) + E ⊥ (σ ), (6)
The temporal spectra E (σ ) and E ⊥ (σ )/2 are shown in Fig. 4 for all rotation rates. Without rotation, energy is nearly equally distributed among the three velocity components (i.e., E ⊥ 2E ). There is actually a slight overrepresentation of horizontal energy, a consequence of the forcing device geometry that preferentially injects energy in horizontal motions. As the rotation rate increases, the high frequencies remain nearly isotropic (isocomponent), whereas the low frequencies become gradually anisotropic, with E (σ ) nearly flat and E ⊥ (σ ) approaching a power law close to σ -5/3 . This anisotropy is related to the fact that, as increases, the decorrelation frequency (i.e., the frequency below which the spectrum becomes flat) gets significantly smaller for the horizontal velocity (σ dec⊥ = 0.04 ± 0.01) than for the vertical velocity (σ dec , increasing from 0.05 to 0.40 rad s -1 for from 2 to 16 rpm). For the largest rotation rate ( = 16 rpm), there is a clear range of frequencies over which the horizontal spectrum E ⊥ (σ ) follows a σ -5/3 power law. This range gets narrower for decreasing rotation rate . The decorrelation frequency of the vertical velocity appears to scale as σ dec * = σ dec /2 = 0.25 ± 0.05, which becomes evident when plotting the ratio 2E /E ⊥ as a function of the normalized frequency σ * = σ/2 (Fig. 5). In this figure, for σ * > 1, for which no inertial waves can exist, energy is nearly equally distributed among the velocity components (2E /E ⊥ 1). The frequency range over which a power law σ 5/3 * is approached is bounded by the two decorrelation frequencies: on the left by σ dec⊥ and on the right by σ dec . Interestingly, we also observe a small frequency domain σ dec * < σ * 0.6 over which energy in the vertical component is slightly larger than in each horizontal component. Such a slight overrepresentation of the vertical velocity is compatible with the componential anisotropy of an assembly of linear inertial waves: A single plane inertial wave has a componential anisotropy 2E /E ⊥ = 2(1σ 2 * )/(1 + σ 2 * ) (shown as a dotted line in Fig. 5), with a larger rms velocity along the vertical than along any horizontal direction for σ * < 1/ √ 3 0.6. As can be seen in Fig. 2, for the experiment under rapid rotation, nearly 90% of the energy is contained in the 2D vortex mode. As discussed quantitatively in Campagne et al. [26], this strong 2D nature of the flow drives an inverse cascade of energy for the horizontal velocity and a direct cascade of energy for the vertical velocity. The horizontal velocity consequently exhibits slow dynamics, while the vertical velocity fluctuations are found at higher frequencies (Fig. 4). This behavior is consistent with the usual phenomenology of rapidly rotating turbulence:

The flow becomes approximately 2D at low frequency and the vertical velocity behaves as a passive scalar, which is stretched and folded by the horizontal velocity. This produces thin vertical layers swept by the horizontal flow, yielding rapidly changing time series of the vertical velocity. In a similar fashion, the σ -5/3 power law of E ⊥ (σ ) could originate from the stochastic sweeping by the large-scale horizontal flow of a k -5/3 ⊥ spatial spectrum, reminiscent of the inverse energy cascade of 2D turbulence.

V. SPATIOTEMPORAL ANALYSIS

A. Spatiotemporal correlations

We now turn to a combined spatiotemporal analysis of the PIV time series, focusing on the signature of inertial waves in terms of dimensional anisotropy. This signature is sought here in terms of characteristic horizontal and vertical scales of the turbulent structures as a function of their frequency. For this, we define the frequency-dependent two-point spatial correlation of the temporal Fourier transform of the velocity field

R(r,σ ) = 2π T ũi (x,σ ) ũ * i (x + r,σ ) + c.c. , ( 7
)
with the asterisk denoting the complex conjugate (here again Welch's method is used to improve convergence). Instead of the spectra considered in Refs. [38,39], we compute spatial correlations, because the former are sensitive to finite-size effects arising from the PIV field being of limited extent compared to the largest flow structures. The correlation [START_REF] Hopfinger | [END_REF] probes the energy distribution among vector separations r for each frequency σ . Its angular average provides an estimate for the cumulative energy from scale r = |r| to r = ∞ for turbulent motions of frequency σ . The single-point limit of this correlation is the temporal energy spectrum (4), i.e., E(σ ) = R(r = 0,σ ).

Maps of the normalized correlation R(r,σ )/E(σ ) are plotted in Fig. 6 for = 0 and 16 rpm at four selected angular frequencies σ . In the nonrotating case, the iso-R lines are approximately circular at all scales and frequencies, indicating the overall isotropy of turbulence. The strongly peaked correlation that develops around r = 0 as σ increases indicates that rapid turbulent fluctuations are found at small scales only. In the rotating case, the iso-R lines evolve gradually from quasivertical at small frequency ("cigar" anisotropy) to more horizontal for σ * ∼ 1 ("pancake" anisotropy), with σ * = σ/2 the normalized frequency. The cigar anisotropy observed at σ * 1 is consistent with the 2D3C vertical invariance predicted by the Taylor-Proudman theorem for vanishing frequency: It corresponds to the zero-frequency limit of Eq. ( 1) for a nearly horizontal wave vector. Similarly, the tendency towards pancake anisotropy, observed for σ * → 1 and sufficiently large r ⊥ , is consistent with the nearly vertical wave-vector limit of Eq. [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF].

A natural way to characterize the frequency-dependent anisotropy would be to compute integral scales along and normal to the rotation axis at each frequency. Here we consider a finer approach, which also takes into account the scale dependence of this anisotropy: For each frequency σ and horizontal scale r ⊥ , we identify the vertical scale (r ⊥ ,σ ) at which the correlation along the vertical axis is equal to the one at r ⊥ e ⊥ , i.e., such that R(r = e ,σ ) = R(r = r ⊥ e ⊥ ,σ ). In practice, we compute as the vertical semiaxis obtained from the fit of the iso-R line defined by R(r,σ ) = R(r ⊥ e ⊥ ,σ ) with an ellipse of prescribed horizontal semiaxis r ⊥ . This method allows us to filter out the noise in the iso-R lines at small R. It also extends the analysis to values of larger than the PIV field height ( z = 140 mm), which is useful at small σ * for nearly vertically invariant R. Poor fits defined by a correlation coefficient less than 0.9 or such that is larger than 2 z are discarded. We finally define the scale-and frequency-dependent anisotropy factor as

A(r ⊥ ,σ ) = r ⊥ (r ⊥ ,σ ) . ( 8
)
It is equal to 1 for isotropic turbulence, to 0 for vertically invariant (2D3C) turbulence, and to ∞ for horizontally invariant (1D2C) turbulence.

If the anisotropy of the two-point correlation R at frequency σ * 1 is governed by linear inertial waves, we expect A to be independent of the scale and to be set by the dispersion relation [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF]. A simple estimate, assuming that the wave vector k is in the vertical measurement plane and identifying on dimensional grounds ∼ k -1 and r ⊥ ∼ k -1 ⊥ , yields a frequency-dependent anisotropy factor A IW (σ * ) (σ -2 * -1) -1/2 . Considering now an assembly of inertial waves with an axisymmetric wave-vector distribution, an analytic computation (given in the Appendix) leads to a very similar result

A IW (σ * ) = 2 σ -2 * -1 . ( 9
)
In Figs. 7(a)-7(c) we compare the anisotropy factor A measured in the four rotating experiments to the inertial-wave prediction (9) at three horizontal scales r ⊥ = 9, 18, and 50 mm. For all scales and σ * 1, A is an increasing function of σ * , confirming that slow fluctuations are more vertically elongated than fast fluctuations. We find that the inertial-wave prediction (9) provides a good description of the data at large horizontal scales and large rotation rate. For such large scales (r ⊥ 50 mm), the anisotropy factor is no longer accessible for σ * < 0.1 because it corresponds to much larger than the height of the PIV field. On the other hand, at smaller horizontal scale the prediction (9) fails, with small frequencies more isotropic than predicted by the inertial-wave argument.

Because of the moderate Reynolds and Rossby numbers of the present experiment, two effects may be considered to explain why large scales follow the inertial-wave prediction whereas small scales do not: viscous damping and sweeping of small scales by the velocity at larger scales. Viscosity introduces an imaginary term iν|k| 2 in the dispersion relation [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF] without modifying its real part. Waves (σ,k) such that |k| r -1 ν are therefore damped, with r ν = √ ν/σ a viscous cutoff. This viscous cutoff is of the order of 1 to 10 mm for the normalized frequencies σ * = σ/2 in the range [10 -2 ,1] considered in Fig. 5. However, since viscous damping affects the wave amplitude without modifying the wave-vector components, its should not affect the anisotropy. We therefore focus in the following on the sweeping effect.

B. Sweeping effect

Sweeping corresponds to the advection of the waves by the large-scale flow, which leads to a modification of their apparent frequency. An inertial wave propagating in a time-independent uniform flow U has a Doppler-shifted frequency

σ = σ i + k • U, (10)
where σ i is the intrinsic frequency given by ( 1) and σ is the frequency at which the wave is detected in the frame of the rotating tank. In our experiment, the energetic large-scale 2D flow may be thought of locally as a uniform sweeping flow U that evolves slowly in time, inducing a scrambling of the waves' spatiotemporal signature. An order of magnitude of the typical Doppler shift can be estimated by k ⊥ U ⊥rms , where U ⊥rms is the root-mean-square horizontal velocity. A key difference between Eqs. ( 10) and ( 1) is that the frequency σ now depends on the magnitude of k, with small-scale waves more affected by sweeping.

For an ensemble of inertial waves with axisymmetric wavevector statistics, the intrinsic frequency σ i can be related to the anisotropy through Eq. ( 9). Substituting the corresponding expression into (10) and estimating the Doppler-shift term on dimensional grounds, we obtain

σ 2 √ 1 + 2A -2 + C U ⊥rms r ⊥ , ( 11
)
where C is a constant of order unity. This indicates that the parameter

N = 2 r ⊥ U ⊥rms 1 + 2/A 2 (12)
should be a unique function of the sweeping parameter S = U ⊥rms /σ r ⊥ . Here, N corresponds approximately to the intrinsic frequency of inertial waves rescaled by the advective time r ⊥ /U ⊥rms , whereas S is the observed period of the waves, rescaled by the same advective time. Figure 8 confirms this picture: The data for different values of , r ⊥ , and σ collapse onto a master curve N = f (S). This collapse indicates that sweeping is indeed responsible for the departure from the inertial-wave prediction at small frequencies and/or small scales. The expected asymptotic behavior for a small sweeping parameter is N 1/S, which corresponds to the prediction (9) for an axisymmetric ensemble of nonswept inertial waves. The data are in quantitative agreement with this small-S prediction, shown as a dashed line in Fig. 8. For large S, Eq. (11) indicates that N should asymptotically approach a constant value N C, which again is compatible with the data.

The master curve in Fig. 8 has the following simple interpretation: High-frequency or large-scale waves are hardly affected by sweeping. The Doppler-shift term is negligible compared to their intrinsic frequency and their location in a space-time energy distribution is given by the dispersion relation [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF]. This is the low-S behavior in Fig. 8. By contrast, when focusing on low frequencies σ or small scales in the frame of the tank, one measures the inertial waves with intrinsic frequency σ i = σ , but one also detects many waves with σ i = σ that are Doppler shifted back to frequency σ by the advective term in (10). The anisotropy measured at low frequency σ therefore results from strongly swept inertial waves with various intrinsic frequencies and the information from the dispersion relation [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF] is lost in the space-time correlation. The limit S 1 corresponds to frequencies σ that are much lower than the inverse advective time. In this σ → 0 limit, one detects mostly waves with σ i σ that are Doppler shifted by the horizontal flow in such a way that they are almost steady in the frame of the tank: This is a σ -independent regime that corresponds to the large-S plateau in Fig. 8.

VI. CONCLUSION

In the present experiment, the anisotropy of the turbulent energy distribution at a given spatiotemporal scale (r ⊥ ,σ ) is well described by the inertial-wave dispersion relation at high frequency and/or large scale only. The smaller-scale waves are subject to intense sweeping by large-scale turbulent motions contained predominantly in the 2D vortex mode. This conclusion is compatible with the numerical findings of Clark di Leoni et al. [38], who also identify the sweeping time scale as the relevant decorrelation time at small scale.

Such sweeping by the 2D mode has strong implications for wave-turbulence theories. Indeed, most waves do not follow the inertial-wave dispersion relation and the assumptions of weak-turbulence theory break down even at the linear stage in wave amplitude: Instead of the dispersion relation [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF], the linear problem consists in determining the evolution of waves embedded in a turbulent 2D flow. This is a formidable task in general because the 2D flow is space and time dependent: In the discussion of our data, we simplified the problem by assuming that the 2D flow is at much larger scales and slower frequencies than the waves, therefore including it as a simple Doppler-shift term in the dispersion relation.

We conclude with a discussion on the dimensionality of the forcing. In the present experiment, the flow is driven by vertically invariant flaps: Such a quasi-2D forcing device enhances two dimensionalization and the resulting sweeping of the 3D flow structures. Nevertheless, accumulation of energy in the 2D mode is a robust feature of rotating turbulence, which takes place for arbitrary forcing geometry, even if the forcing does not input energy directly into the 2D mode. A careful and extensive numerical study of this issue has been recently reported for the fully 3D Taylor-Green forcing [47]: For rapid global rotation and low viscosity, energy accumulates in the 2D mode until the Rossby number based on the velocity of this 2D flow is of order unity. If these findings are confirmed, the sweeping of the most energetic 3D structures would be an inevitable outcome of this accumulation of 2D energy.

FIG. 1 .

 1 FIG. 1. (Color online) Experimental setup: (a) side view and (b) top view. An arena of ten pairs of flaps forces a turbulent flow in the central region of a water tank mounted on a rotating turntable.A laser sheet illuminates a vertical slice through a horizontal glass lid covering the fluid. Two-dimensional three-component velocity measurements are performed using stereoscopic particle image velocimetry in a vertical square domain of size x × z = 14 × 14 cm 2 , shown as a dashed square in (a).

3 )FIG. 2 .

 32 FIG. 2. (Color online) (a) Energy and (b) energy dissipation rate per unit mass for the 2D and 3D modes as a function of the rotation rate . In both figures, the first data points (shown with arrows, at arbitrary abscissas) correspond to the nonrotating case = 0.

FIG. 3 .

 3 FIG. 3. (Color online) Temporal energy spectrum E(σ ) as a function of the angular frequency σ for different rotation rates . Here σ 0 indicates the frequency of the opening-and-closing cycle of the flaps. The Coriolis frequency 2 is highlighted with closed symbols.

  with E (σ ) = 4π | ũz (x,σ )| 2 /T the spectrum of the vertical velocity and E ⊥ (σ ) = 4π | ũx (x,σ )| 2 /T + 4π | ũy (x,σ )| 2 /T the spectrum of the horizontal velocity.

FIG. 4 .

 4 FIG. 4. (Color online) Temporal energy spectra of the vertical E [light gray (blue)] and horizontal E ⊥ /2 (black) velocity components for different rotation rates . For visibility, there is a vertical shift by a factor of 10 3 between couples of curves at different . The dashed lines show power laws σ -5/3 .

FIG. 5 .

 5 FIG. 5. Componential anisotropy ratio as a function of the normalized frequency σ * = σ/2 . Isotropy is indicated by the horizontal dashed line at 2E /E ⊥ = 1. The dotted line indicates the prediction for a plane inertial wave, i.e., 2E /E ⊥ = 2(1σ 2 * )/(1 + σ 2 * ).

σ

  FIG. 6. (Color online) Maps of the normalized two-point correlation R(r,σ )/E(σ ) in the vertical plane (r ⊥ ,r ) for = 0 (left) and = 16 rpm (right), at four frequencies σ = 0.1,0.7,2.4,3.3 rad s -1 . In the rotating case, the corresponding normalized frequencies are σ * = σ/2 = 0.03,0.20,0.71,0.98. Isocontour lines for R(r,σ )/E(σ ) < 0.2 are noisy and are not shown.

FIG. 7 .

 7 FIG. 7. Anisotropy factor A (8) as a function of the normalized frequency σ * = σ/2 at different rotation rates (same symbols as in Fig. 3), for three horizontal scales (a) r ⊥ = 9, (b) r ⊥ = 18, and (c) r ⊥ = 50 mm. The solid line represents the inviscid inertial-wave prediction A IW (9).

1 NS

 1 FIG. 8. Rescaled intrinsic frequency N (12) as a function of the sweeping parameter S. The symbols indicate the different rotation rates and are the same as in Fig. 3. The dashed line N ∼ 1/S shows the low-S prediction for (nonswept) ensembles of inertial waves [Eq. (12)]. For S 1 one sees the signature of the dispersion relation (1), while for S 1 one detects swept inertial waves.

  from the six terms of the velocity gradient tensor accessible in the 3C2D measurements,

	ν 2	∂u x ∂x		2	+ 2	∂u y ∂x	2	+ 2	∂x ∂u z	2
	+	∂u x ∂z	2	+	∂u y ∂z	2	+	∂u z ∂z

2

.

ACKNOWLEDGMENTS

We acknowledge P. Augier, P. Billant, and J.-M. Chomaz for kindly providing the flap apparatus and A. Aubertin, L. Auffray, C. Borget, and R. Pidoux for their experimental help. This work was supported by the ANR Grant No. 2011-BS04-006-01 "ONLITUR". This work is supported by "Investissements d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM). F.M. acknowledges the Institut Universitaire de France.

APPENDIX: ANISOTROPY FACTOR FOR A STATISTICALLY AXISYMMETRIC DISTRIBUTION OF INERTIAL WAVES

We compute the anisotropy factor A for an ensemble of independent plane inertial waves, with axisymmetric wavevector statistics. The temporal Fourier transform of the velocity field reads

where a(k,σ ) is the space-time Fourier amplitude of the velocity field at wave number k and frequency σ . The twopoint velocity correlation at frequency σ (7) can be written

where the angular brackets is the space average and r is a separation vector inside the PIV plane. Introducing spherical coordinates with vertical polar axis, we denote by ϕ the azimuthal angle between k and the vertical PIV plane. The argument of the cosine becomes

Let us first consider an ensemble of inertial waves having the same wave number |k| = k. For a given reduced frequency σ * = σ/2 , the dispersion relation determines the ratio |k |/k ⊥ and because of statistical axisymmetry |a(k,σ )| 2 is independent of ϕ. The spatial correlation at frequency σ becomes

where J 0 is the Bessel function of the first kind and G(k,σ ) is a prefactor proportional to the squared amplitude of the waves at wave number k and frequency σ . For a given frequency σ and horizontal scale r ⊥ , the vertical scale (r ⊥ ,σ ) is determined from the isolines of R in the (r ⊥ ,r ) plane. An isoline of R starting on the horizontal axis at r ⊥ intersects the vertical axis at r = (r ⊥ ,σ ). According to expression (A4), such isolines connecting the two axes exist provided the argument of the Bessel function is smaller than the first zero of this function, i.e.,

where J 0 (C 0 ) = 0 (C 0 2.40). Equating expression (A4) computed for (r ⊥ ,r = 0) and for (r

and an anisotropy factor

This anisotropy factor depends very weakly on kr ⊥ : It is minimum for low kr ⊥ , where Taylor expansion for kr ⊥ 1 leads to

whereas it is maximum for the maximum value of r ⊥ allowed by (A5), which gives

Because the numerical prefactors in (A8) and (A9) differ by less than 10%, we can say that the anisotropy factor of this statistically axisymmetric distribution of inertial waves is given by expression (A8) within 10% accuracy. Because this anisotropy factor is almost independent of k, we finally expect it to be approximately given by expression (A8) for a realistic superposition of inertial waves with different wave numbers k.


D Monnet

J Ninin

L Jaulin
email: luc.jaulin@ensta-bretagne.fr

Computing an Inner and an Outer Approximation of the Viability Kernel *

Keywords:

The viability kernel corresponds to the set of all state vectors of a controlled dynamic system that are viable, i.e., such that there exists an input such that the system will not enter inside a forbidden zone. In this paper, we propose a method which computes an inner and an outer approximation of the viability kernel in a guaranteed way. Our method is based on interval analysis and uses the notions of V-viability and capture basin. We illustrate our approach on the car on the hill problem. A software package has been developed to solve any 2D-problem.

Introduction

Safety verification of controlled systems has been approached with different tools such as viability theory [START_REF] Aubin | Viability Theory[END_REF], reachability analysis, or the concept of barrier function (see e.g., [START_REF] Bouissou | Computation of parametric barrier functions for dynamical systems using interval analysis[END_REF]). Safety is often expressed as a set of constraints in which the system must stay. For example, safety verification problems may consist of ensuring a safe configuration during the landing [START_REF] Bayen | Aircraft Autolander Safety Analysis Through Optimal Control-Based Reach Set Computation[END_REF] or the take off [START_REF] Seube | Aircraft take-off in windshear: A viability approach[END_REF] of an airplane, or collision avoidance [START_REF] Desilles | Collision analysis for an UAV[END_REF] with other aircraft. When we have an input to the system, we may want to find a controller which guarantees that the system is safe for all conditions. The set of all state vectors such that a system can stay (for the right input) within a set of constraints is called the viability kernel.

The computation of the viability kernel has been addressed by several methods. For low-dimensional problems and non-linear systems, most methods are based on gridding the state space [START_REF] Deffuant | Approximating Viability Kernels With Support Vector Machines[END_REF][START_REF] Saint-Pierre | Approximation of the viability kernel[END_REF], but due to the finite precision of the computer, gridding methods cannot be considered as guaranteed. For high-dimensional systems, the viability kernel can be approximated using methods based on Lagrangian methods [START_REF] Maidens | Lagrangian methods for approximating the viability kernel in high-dimensional systems[END_REF] for linear systems, or invariant sets [START_REF] She | Brief paper: Computing an invariance kernel with target by computing Lyapunov-like functions[END_REF] for polynomial systems. The viability kernel can also be approached using tools from the reachability analysis. Reachability analysis consists of computing the set of all state vectors that the dynamic system can reach from a given initial state. This problem has been considered with Hamilton-Jacobi equations [START_REF] Mitchell | Validating a Hamilton-Jacobi Approximation to Hybrid System Reachable Sets[END_REF], and was most recently approached with interval analysis [START_REF] Rego | Determination of inner and outer bounds of reachable sets through subpavings[END_REF]. If the target set belongs to a safe region, the viability kernel can be approximated by computing the set of state vectors from which the system can reach the target [START_REF] Kaynama | Computing the Viability Kernel Using Maximal Reachable Sets[END_REF]. Now, for our problem of computing the viability kernel, we cannot assume a priori that such a target set, included in the viability kernel, is available. Therefore, a subset of the viability kernel must be first computed in a reliable way.

In this paper, we consider low-dimensional and non-linear dynamic systems, and we propose to compute a guaranteed approximation of the viability kernel with interval analysis tools [START_REF] Jaulin | Applied Analysis[END_REF]. Instead of manipulating a set of points as gridding methods do, interval methods consider connected sets of state spaces. Using interval techniques, computations are reliable and guaranteed, which allows us to compute an inner and an outer approximation of the viability kernel.

This paper is organized as follows. Section 2 contains the notation and gives the definition of the viability kernel and the capture basin. In Section 3, we present the concept of V -viability, which will allow us to find subsets of the viability kernel. Section 4 provides theoretical results to compute the capture basin and Section 5 gives an algorithm to get an inner and an outer approximation of the viability kernel. In Section 6, we illustrate our approach with a two-dimensional example. Section 7 concludes this paper.

Notation and Definitions

In this paper, an interval is a continuous set, denoted as follows: [x] = [x, x], x ≤ x, with x the lower bound and x the upper bound. A vector of intervals [x] = ([x1], . . . , [xq]) is commonly called a box. A dynamic system S is defined by the following differential equation :

ẋ(t) = f (x(t), u(t)), u(t) ∈ U, (1)
where:

• x(t) ∈ R n is an evolution of the state variables,

• U is a compact subset of R m which represents the set of admissible values for the control,

• u ∈ U = {u : R + → U} is the control vector,

• f : R n ×U → R m is the evolution function of S. We assume that f is continuous, locally Lipschitzian and bounded in R n × U.

Let ϕ be the flow map of S, i.e., with the initial condition x0 = x(0) and a control function u(t), the system S reaches the state ϕ(t, x0, u) at time t.

Let us define the viability kernel and the capture basin of S.

Definition 2.1 Let S be a dynamic system and let K ⊆ R n be a compact set. The viability kernel of K under S is the set ViabS (K) of initial states x ∈ K from which at least one evolution does not leave K for all t ≥ 0. We have

ViabS (K) = {x ∈ K | ∃u ∈ U, ∀t ≥ 0, ϕ(t, x, u) ∈ K}.
We now define the capture set (also called robust invariant set), which is closely related to the concept of viability kernel.

Definition 2.2 Let T ⊂ K be a target. The capture basin of T viable in K under S is the set CaptS (K, T) of initial states x ∈ K from which at least one evolution of S in K reaches the target T in a finite time

CaptS (K, T) = {x0 ∈ K | ∃t ≥ 0, ∃u ∈ U , ϕ(t, x0, u) ∈ T, ϕ([0, t], x0, u) ⊆ K. . where ϕ([t1, t2], x0, u) = {x ∈ R n | ∃t ∈ [t1, t2], x = ϕ(t, x0, u)}. ( 2
)
3 V -Viability Definition 3.1 Consider a dynamic system S and a differentiable function

V : R n → R. S is said to be V -viable if ∀x ∈ R n such that V (x) = 0, ∃u ∈ U, ∇V (x) • f (x, u) < 0.
Figure 1 illustrates a V -viable system. For each point of the curve V (x) = 0 at least one potential evolution of S points strictly inward into the gray set, which represents V -1 (R -). This is the geometrical interpretation of ∇V (x) • f (x, u) < 0. The definition of V -viability implies the resolution of a quantified constraint satisfaction problem that can be solved with cylindrical decomposition [START_REF] Caviness | Quantifier Elimination and Cylindrical Algebraic Decomposition[END_REF] if V is polynomial or with methods based on interval analysis [START_REF] Hladík | Efficient solution of a class of quantified constraints with quantifier prefix Exists-Forall[END_REF] in the general case.

Theorem 3.1 Consider K a subset of R n . If S is V -viable and V -1 (R -) is a compact subset of K, then V -1 (R -) ⊆ ViabS (K).
Proof: Consider x ∈ K and the subset L of K defined by

L = {x ∈ K|V (x) ≤ 0} = V -1 (R -)
The contingent cone T L (x) to L at x ∈ L, V (x) = 0 can be defined with the gradient of V . According to [1, pp. 123], T L (x) is defined by

T L (x) = {y ∈ K| V (x) = 0 =⇒ ∇V (x), y ≤ 0}
Here we recall the definition of a viability domain [2, pp. 84]:

We shall say that a subset L ⊂ K is a viability domain of K if and only if

∀x ∈ L, f (x, U) ∩ T L (x) = ∅
As S is V -viable, we have:

∀x such as V (x) = 0, ∃u ∈ U, V (x), f (x, u) < 0 =⇒ f (x, U) ∩ T L (x) = ∅
Therefore, L is a viability domain. Moreover, according to the Local viability theorem [2, pp. 91], we have:

x0 ∈ L =⇒ ∃u ∈ U such as ∀t ≥ 0, ϕ(x0, t, u) ∈ K.
Thus, by definition, L ⊆ ViabS (K).

Theorem 3.1 can be used to prove that a set is inside the viability kernel of K. Figure 1 illustrates this by the gray set, which belongs to ViabS (K). To find a viable subset of the viability kernel, a function V must be chosen so that S has a high probability to be V -viable. Such a function can be found using Lyapunov theory [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF][START_REF] Slotine | Applied Nonlinear Control[END_REF]. Some Lyapunov methods can approximate attraction domains of S [START_REF] Delanoue | An algorithm for computing a neighborhood included in the attraction domain of an asymptotically stable point[END_REF][START_REF] Ratschan | Providing a basin of attraction to a target region of polynomial systems by computation of Lyapunov-like functions[END_REF]. If S is Vviable and the second condition of Theorem 3.1 is checked, V -1 (R -) is proved to be a subset of ViabS (K).

Capture Basin of Viable Sets

In Section 3, a method to compute subsets of the viability kernel from the notion of V -viability is provided. These subsets represent a first approximation of ViabS (K). In order to obtain a better approximation, these viable subsets are enlarged with their capture basins.

The following theorem states that the capture basin of a subset of ViabS (K) is also a subset of ViabS (K). Theorem 4.1 Let S be a dynamic system, K be a closed subset of the state space of S, and T be a subset of K. If T ⊂ ViabS (K), then CaptS (K, T) ⊂ ViabS (K).

Algorithm to Compute Viab S (K)

In this section, we suppose that we have found a set E ⊆ ViabS (K) using techniques of Section 3. We present an algorithm that computes Vinner, a guaranteed inner approximation of CaptS (K, E), and H, a guaranteed approximation of the complement of ViabS (K) in K. This algorithm is an extension of the one presented in [START_REF] Lhommeau | Capture basin approximation using interval analysis[END_REF] to compute the viability kernel. We choose to subdivide the state space into boxes, so that we can easily use Proposition 4.2. Figure 2a shows the initial problem. Its representation with a subpaving is displayed in Figure 2b. for all

[x i ] in S do 3:
Choose u ∈ U.

4: if [ϕ]([0, t], [x i ], u) ⊆ K and [ϕ](t, [x i ], u) ⊆ V inner then 5: V inner := V inner ∪ [x i ], S := S \ [x i ],
6:

else if [ϕ](t, [x i ], U) ⊆ ((R n \ K) ∪ H) then 7: H := H ∪ [x i ], S := S \ [x i ],
8:

end if 9:
end for 10:

Bisect boxes of S. 11: end while 12: return V inner and H. Compared to the algorithm of [START_REF] Lhommeau | Capture basin approximation using interval analysis[END_REF], our version contains technical and computational improvements to accelerate the convergence:

• The enclosure of [ϕ] is computed using the technique of [START_REF] Chapoutot | Validated explicit and implicit Runge-Kutta methods[END_REF], which improves the quality of the bounds and increases the chance to include a box in Vinner or H.

• The data structure of Vinner and H is carefully implemented. A regular paver, represented by a binary tree, is used to perform fast intersections and set unions .

• In Algorithm 1, the subboxes composing S are divided only if it is not possible to prove that they belong to Vinner or H. In the algorithm presented in [START_REF] Lhommeau | Capture basin approximation using interval analysis[END_REF], S is entirely dividied into small boxes with the same size, then it tries to prove if they belong to Vinner or H. Algorithm 1 divides boxes only if needed; this allows it to include large boxes directly in Vinner or H.

These points aim to break down the exponential complexity of this algorithm.

6 Application to the Car on the Hill Problem

We illustrate our approach with the car on the hill problem. In this application, we want the car to stay on a landscape represented by the parametric function

g : s → -1.1 1.2 cos(1.2s) + 1.2 1.1 cos(1.1s) 2 ,
where s ∈ [0, 12] denotes the longitudinal variable. The function g is plotted on Figure 4. The acceleration of the car can be controlled within a limited range. We also consider a friction force that slows down the car. The dynamic system S associated with our problem is described by ẋ1(t) = x2(t), ẋ2(t) = -9.81 sin( ġ(x1(t))) -0.7x2(t) + u(t), where x1 = s represents the position of the car on the landscape, x2 = ṡ its velocity, and u ∈ [-2, 2] is the control function. The viability problem is formulated as follows:

we want to keep x1 between values 0 and 12 with a limit on the velocity of the car between -6 and 6. Thus, the set of constraints K is defined by

K = {x | x1 ∈ [0, 12], x2 ∈ [-6, 6]}.
We limit the precision on s and ṡ to 0.1. To compute ViabS (K) we first apply the method of Section 3 to find a subset E of ViabS (K). To do so, S is linearized around an equilibrium point of S. Then, a quadratic Lapunov function is computed for the linearized system. A function V is created from the Lyapunov function so that S is V -viable. This procedure is repeated for several equilibrium points. Figure 5 illustrates the result obtained. The elliptic shapes of subsets of ViabS (K) in Figure 5 come from the quadratic expressions of the Lyapunov functions. Moreover, the centers of these viable subsets correspond to a position of the car at the tops or bottoms of the hills with a null velocity, where the car is in an equilibrium state.

The computation took 10 seconds on a 2.5 GHz Intel Core i5-2450M processor with 6 Gb RAM. Next, we apply Algorithm 1 to compute an inner and outer approximation of ViabS (K). The result is shown in Figure 6. The CPU time is 1h 23min. The darkgray set corresponds to Vinner and is the largest set that we could prove to be viable in a guaranteed way. Nothing could be proved for the boxes of the white set. The light gray set is proved to be outside the viability kernel. The union of the white set and the dark-gray set corresponds to Vouter. Figures 5 and6 have been obtained with a solver available at http://www.ensta-bretagne.fr/monnet/Viabibex/.

Conclusion

This paper proposes an interval method to approximate the viability kernel of a nonlinear dynamic system. Our approach merges the concepts of V-viability and reachability analysis. An inner approximation of the viability kernel is computed on an infinite time horizon. Interval analysis tools ensure the reliability and provide a guarantee for all the results. The non-linear constraints on the evolution function enable application of our method to a large variety of problems. Numerical results for the car on the hill problem demonstrate the feasibility of our approach. Moreover, our method can be generalized to an n-dimensional problem.

Figure 1 :

 1 Figure 1: The vector field is displayed for several controls along the level curves of V .

Figure 2 : 1

 21 Figure 2: Subpaving representing the initial inner approximation for ViabS (K)

Figure 3 :

 3 Figure 3: Illustration of the method for characterizing ViabS (K)

Figure 4 :

 4 Figure 4: Car on the hill problem. We aimed to avoid the car falling off the cliffs at s = 0 and s = 12.

Figure 5 :

 5 Figure 5: Subsets of ViabS (K) are found with method of Section 3

Figure 6 :

 6 Figure 6: Approximations of the viability kernel.

Proof: Consider x0 ∈ CaptS (K, T). From the definition of CaptS (K, T), we have ∃t1 ≥ 0, ∃u1 ∈ U , ϕ(t1, x0, u1) ∈ T and ϕ([0, t1], x0, u1) ⊂ K.

Then, ∀t > 0, ϕ(t, x0, u) ∈ K. So, according to the definition of ViabS (K), we have x0 ∈ ViabS (K).

Proof: This proof is trivial since it comes from the definition of ViabS (K). Proposition 4.1 may be used to prove that a state belongs to CaptS (K, E) or does not belong to ViabS (K).

Proof: The proof of (i) and (ii) are the same of those of Proposition 4.1 considering sets of states variable instead of unique point.

Methods using interval arithmetic are able to compute a guaranteed enclosure of the flow map ϕ from an initial box [x0] of the state space, see [START_REF] Chapoutot | Validated explicit and implicit Runge-Kutta methods[END_REF][START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] for details. Using these techniques, it is possible to compute

• an enclosure [ϕ](t, [x0], U) of all evolutions ϕ(t, x, u) with all initial states x ∈ [x0] and all the possible control functions u ∈ U .

Such enclosures are generally overestimated. But if the size of the box [x0] is small enough, this estimation can be accurate.