W Hamouche
email: walid.hamouche@dalembert.upmc.fr

C Maurini

A Vincenti
email: angela.vincenti@upmc.fr

S Vidoli

Basic criteria to design and produce multistable shells

Keywords: Shallow shells, Shape morphing, Bistability, Tristability

A shell can have multiple stable equilibria either if its initial curvature is sufficiently high or if a suitably strong pre-stress is applied. Under the hypotheses of a thin and shallow shell, we derive closed form results for the critical values of curvatures and pre-stresses leading to bistability and tristability. These analytical expressions allow to easily provide guidelines to build shells with different stability properties.

Introduction

The highly non-linear behaviour of shell structures along with anisotropic material properties and prestresses can induce interesting phenomena, such as multistability, i.e. the existence of multiple distinct configurations of stable equilibria. Similar complex non-linear behaviours recently raise a growing scientific interest; they are the object of a large number of research works ranging from the analysis of biological systems, such as the study of multistable natural structures [START_REF] Forterre | How the venus flytrap snaps[END_REF][START_REF] Koehl | How kelp produce blade shapes suited to different flow regimes: a new wrinkle[END_REF] or the problem of the growth of elastic membranes, see for instance [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Klein | Shaping of elastic sheets by prescription of non-euclidean metrics[END_REF][START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF], to the analysis and design of morphing structures, i.e. structures that can largely change their shape in order to adapt to different functioning regimes [START_REF] Portela | Analysis of morphing, multi stable structures actuated by piezoelectric patches[END_REF][START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF][START_REF] Norman | Novel multistable corrugated structures[END_REF][START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF][START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF][START_REF] Seffen | Prestressed morphing bistable and neutrally stable shells[END_REF], and more generally to the study of slender structures undergoing large displacements, see for instance [START_REF] Seffen | Heating of a uniform wafer disk[END_REF][START_REF] Modes | Gaussian curvature from flat elastica sheets[END_REF].

The study of the non-linear response of shell structures is generally perceived as a complex problem of structural mechanics [START_REF] Calladine | Theory of shell structures[END_REF] not prone to a simple analytical treatment. However, as pointed out in [START_REF] Pirrera | On the thermally induced bistability of composite cylindrical shells for morphing structures[END_REF], it is extremely difficult to extrapolate informations on the global non-linear behaviour of shells only from the results of direct finite-element simulations. In fact, the numerical analysis of fully non-linear shells is a difficult task and only parametric studies can be conducted for limited combinations of geometry, material properties and loading conditions, see for instance [START_REF] Gigliotti | Loss of bifurcation and multiple shapes of thin [0/90] unsymmetric composite plates subject to thermal stress[END_REF].

The use of simplified analytical models is of paramount importance, because it can provide a global picture of the qualitative behaviour of non-linear shells. This may be obtained through reduced order models assuming restrictive hypotheses on the admissible deformation modes and on the order of the geometrical nonlinearities. The application of the Föppl-von Karman (FvK) kinematics has represented an important contribution to the study of shallow shells [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Klein | Shaping of elastic sheets by prescription of non-euclidean metrics[END_REF][START_REF] Giomi | Multi-stability of free spontaneously curved anisotropic strips[END_REF], and specifically a deeper insight of the stability of anisotropic shells was possible thanks to the use, within FvK models, of the syimplifying assumption of Uniform Curvature (UC) [START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF][START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF][START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF][START_REF] Seffen | Prestressed morphing bistable and neutrally stable shells[END_REF][START_REF] Guest | A zero-stiffness elastic shell structure[END_REF]. Despite its simplicity, such assumption allowed to foresee and experimentally demonstrate the existence of three stable equilibria for orthotropic shells [START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF][START_REF] Coburn | Tristability of an orthotropic doubly curved shell[END_REF].

Although the accepted domain of validity of the FvK model is limited to plates and shallow shells, many of the works cited above seem to indicate that, at least for shells loaded only by inelastic deformations, the FvK can furnish results in good agreement with experiments for an impressive large range of magnitude of curvatures (both initial geometric curvatures, as well as elastic or inelastic curvatures induced by external loads), even when the model is further simplified by the assumption of uniform curvatures [START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF][START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF][START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF]. In a recent interesting work, [START_REF] Lewicka | The foppl-von karman equations for plates with incompatible strains[END_REF] rigorously show by asymptotic expansion and Gamma-convergence that FvK is indeed the pertinent model to include the dominant geometrically non-linear effects in plates or shallow shells with inelastic deformations. However a better numerical validation and theoretical investigation of the range of validity of the results obtained by uniform curvature FvK model seem to be necessary. First steps in this direction were presented in [START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF] and [START_REF] Pirrera | On the thermally induced bistability of composite cylindrical shells for morphing structures[END_REF]. In [START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF] the simplified UC FvK model was applied to the study of the non-linear behaviour of orthotropic plates undergoing large inelastic curvatures, from thermo-elastic or piezo-elastic origin, and the validity of the analytical results was confirmed by means of fully non-linear finite-element simulations. In [START_REF] Pirrera | On the thermally induced bistability of composite cylindrical shells for morphing structures[END_REF] the results obtained with polynomial approximations of FvK and higher order shell models were compared to non-linear finite element simulations.

In the present work, we intend to investigate the applicability of the uniform curvature assumption in order to model and to design multistable shells with high initial and inelastic curvatures. In Section 2 we describe the two basic bifurcations which lead to the multistable behaviour of shells. For the sake of synthesis, we illustrate these concepts in the case of spherical initial and inelastic curvatures. Section 3 illustrates the uniform curvature model and summarizes the useful results for design purposes. Based on the simple "recipes" given in Section 4, in Section 5 we present experimental prototypes for each of the three basic relevant cases of multistability. The analytical results and the simplified models are checked against fully non-linear numerical simulations based on a commercial FE package (Abaqus). Section 6 draws the conclusions. In particular we discuss agreements and failures of the analytical models with respect to the numerical simulations, thus assessing the range of validity of the Uniform Curvature Föppl-von Karman model for the description of the highly non-linear behaviour of shell structures.

Basic phenomena

We study how the stable equilibrium shapes of freestanding shells depend on two distinct classes of control parameters: (i) the curvatures of the shell in its initial stress-free configuration and, (ii) the inelastic curvatures induced by thermal, plastic or hygroscopic loads.

We consider here the simple case in which the natural and inelastic curvatures are uniform in space and spherical (i.e. equal in all directions). Under these hypotheses, they are completely characterised by two independent scalar parameters: the magnitudes c 0 and c i of the initial and inelastic curvatures, respectively. To simplify the analytical treatment, we will further assume that also the current equilibrium configuration has uniform (but generally not spherical) curvatures.

Here and henceforth, we will denote by K x , K y , K xy the three independent components of the symmetric tensor characterising the current curvature in this simplified Uniform Curvature (UC) model. Despite its simplicity, this setting allows us to faithfully describe the rich non-linear behaviour of the shells under the effect of the control parameters (c 0 , c i ) and to design structures with multiple (up to three) stable equilibrium shapes or a continuous set of configurations at the same energy level (neutrally stable shells). The results of the simple UC model are checked against more refined models, as a model with Quadratic Curvature (QC, see Section 3) and fully non-linear shell finite-element models with a large number of degrees of freedom.

Figures 1, 2 and 3 show bifurcation diagrams describing the effect of the natural and inelastic curvatures. Details on the procedure that allows to draw these graphics are given in Section 3. As three, qualitatively different, stability regions arise, in Section 4 we sketch a procedure to design shells within each one of them. Section 5 describes the actual production of three prototypes corresponding to points B, C and T in Figure 3. In the rest of this Section, we comment on the key qualitative phenomena highlighted in Figures 1232.1 Shells with initial curvature c 0 and vanishing inelastic curvature (c i = 0) Shells with a (spherical) initial curvature c 0 and vanishing inelastic curvature are characterised by a natural stress-free configuration with

K x = K y = c 0 , K xy = 0. (1)
For c 0 = 0 the natural shape of the shell is a spherical sector of radius 1/c 0 . For c 0 = 0 the shell degenerates to a flat plate. Because of geometrical non-linearities, stable equilibrium configurations different from the natural configuration above may exist. Figure 1 shows the emergence of alternative stable equilibria when the absolute value of the natural curvature c 0 exceeds a critical threshold. Shells which are sufficiently curved in their initial configuration possess an additional (not stress-free) stable everted configuration. black points refer to Abaqus FE results, solid black and brown respectively to the UC and QC models, whilst light gray refers to unstable equilibria.

2.2 Shells with inelastic curvature c i and vanishing initial curvature (c 0 = 0)

The inelastic curvature affects the constitutive relation between the bending moments and the (geometric) curvatures of the shell. In a linear constitutive model, the bending energy is quadratic with respect to the misfit between the current curvature K and the spherical inelastic curvature c i . The elastic energy density due to bending vanishes only when K is spherical and equal to c i . The plate would like to conform in a spherical shape to minimize the bending energy. However, because of the geometrical compatibility, an initially flat plate can require a significant amount of extensional energy to deform as a sphere. It would prefer to deform in cylindrical shapes. This latter effect, related to geometrical nonlinearities, becomes predominant in the large curvature regime. Indeed, when the absolute value of the inelastic curvature c i exceeds a critical threshold, the spherical deformation mode bifurcates toward an almost cylindrical one, as shown in Figure 2. In perfectly isotropic circular discs, the structure has not a preferred bending axis. All the configurations obtained by rotating the axis of the cylinder would be neutrally a stable equilibrium, producing a neutrally stable shell [START_REF] Seffen | Prestressed morphing bistable and neutrally stable shells[END_REF]. However, small imperfections in the shape of the disc, in the constitutive properties, or in the inelastic curvature would modify this idealized behaviour. In practical situations, the disc is bistable in this regime, with two stable equilibria, where the shell is bended in two orthogonal directions (see Figure 2). In conclusion, inelastic curvatures can lead to shells with multiple equilibrium configurations. black points refer to Abaqus FE results, solid black and brown lines (superposed) are obtained with the UC and QC models, whilst light gray refers to unstable equilibria.

Combined effects of initial and inelastic curvatures

The stability diagram is sensibly enriched when combining the effects of initial and inelastic curvatures. Only one half of the plane, namely c i > 0, is shown as for c i < 0 the effect of (c 0 , -c i ) is equivalent to that of (-c 0 , c i ). As to be expected, the two phenomena shown in Figures 1 and2 superpose to produce a richer scenario. For instance choosing the parameters c 0 and c i within the darker gray region produces a shell having both the everted and the bifurcated shapes as stable equilibria. However note that the tristability region shown here is different from the case discussed in [START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF] and [START_REF] Coburn | Tristability of an orthotropic doubly curved shell[END_REF]. In those references the tristable behaviour of an orthotropic shell was induced using only a natural curvature field with K x suitably larger than K y . Here the tristable nature of some shells is due to the contemporary occurrence of natural and inelastic curvatures.

Discrete reduced models and analytical results

In this Section we briefly recall the procedure to deduce, from the Föppl-von Kármán shallow shell equations, reduced models with few degrees of freedom and the main analytical results arising from similar models. For further details the reader is addressed mainly to [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF].

In particular we focus on the results arising from the Uniform Curvature (UC) reduced model that proves to be the most useful one for design purposes.

Föppl-von Kármán assumptions

The Föppl-von Kármán (FvK) assumptions are often used to simplify the fully non-linear Koiter shell model to the case of thin shallow shells. The non-linear coupling between membrane deformations and transverse displacements is still accounted for and, actually, it represents the key to describe many interesting phenomena as the multistable behaviour of such structures. Two displacements fields are taken as primary unknowns: with reference to an initially flat configuration S ⊂ IR 2 these are the in-plane, say u = (u x , u y ), and the transverse or normal, say w, displacement fields 1 . The stable equilibria of the FvK model are found as the local minimizers of the total energy:

U(u, w) = S Ā (e -f ) • (e -f ) 2 + D (k -h) • (k -h) 2 d S, (2)
where e and k are respectively the in-plane distortion and the curvature given e.g. by [START_REF] Mansfield | The Bending & Stretching of Plates[END_REF]:

e := {e x , e y , 2e xy } = =            ∂ux ∂x + 1 2 ∂w ∂x 2 , ∂uy ∂y + 1 2 ∂w ∂y 2 , ∂ux ∂y + ∂uy ∂x + ∂w ∂x ∂w ∂y            , (3)
k := {k x , k y , 2k xy } = ∂ 2 w ∂x 2 , ∂ 2 w ∂y 2 , 2 ∂w ∂x∂y . ( 4
)
1 Sometime the membrane problem, aimed at determining u, is translated in terms of an Airy stress function, say Φ.

Note that the in-plane distortion e includes the secondorder geometric contributions due to transverse displacements, which is the source of the non-linear behaviour; correspondingly the compatibility between e and k requires:

curl curl e := ∂ 2 e y ∂x 2 + ∂ 2 e x ∂y 2 -2 ∂ 2 e xy ∂x∂y = k x k y -k 2 xy =: det k, ( 5
)
where det k means the shell Gaussian curvature and curl is the standard curl operator in S. In Eq. ( 2) Ā and D are respectively the 3 × 3 matrices representing the membrane stiffness and the bending stiffness. Moreover, still in Eq. ( 2), the fields f and h are decomposed into two separate contributions:

f = f 0 + f i , h = h 0 + h i : curl curl f 0 = det h 0 . ( 6
)
The contributions f i and h i describe inelastic membrane and bending deformations as those associated to thermal, plastic, hygroscopic or piezoelectric effects, while f 0 and h 0 give the shape of the initial configuration.

Reduction under the assumption of uniform curvature

A necessary condition to minimize (2) is to solve an elliptic problem on the two-dimensional domain S for the membrane displacement fields (u x , u y ). Noting that this problem is linear and completely defined once assigned the transverse displacement w, an efficient procedure to approximate the FVK energy with few discrete degrees of freedom was proposed in [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF].

In particular, assuming the transverse displacement in the form

w(x, y) = kx x 2 2 + ky y 2 2 + kxy x y (7)
leads to a particularly simple reduced model labeled as UC. This name refers to the fact that ( 7) is tantamount to assume the curvature (4) to be uniform, k = k := { kx , kx 2 kxy }. In [START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF] was proven that, under the hypothesis (7) the minimization of ( 2) is reduced to the minimization of

U (K) = 1 2 D (K -H 0 -H i ) • (K -H 0 -H i ) + 1 2 (det K -det H 0 ) 2 , ( 8
)
in terms of the three state parameters

K = {K x , K y , 2K xy }.
Here det K = K x K y -K 2 xy , and

D := D/ D11 , D := 1 L 2 S D d S, (9)
is an averaged and normalized bending stiffness matrix. The curvatures K := R k, H 0 := R h 0 and H i := R h i are expressed in dimensionless form using the characteristic radius

R = 12 ψ L 2 t , t = 12 D11 / Ā11 . (10)
In ( 9) and ( 10) L 2 is the area of the planform of the shell and t its equivalent thickness.

The scalar number ψ which determines the curvature scaling is of fundamental importance as it measures the ratio between bending and membrane energy. It is function only of the planform shape Ω and of the normalized membrane stiffness A := Ā/ Ā11 . Let us introduce dimensionless coordinates X = x/L and Y = y/L and a corresponding scaled domain Ω. Choosing L so that the area of Ω is 1, the number ψ is computed as

ψ := Ω A -1 Σ • Σ, Σ = ∂ 2 Φ ∂Y 2 , ∂ 2 Φ ∂X 2 , - ∂ 2 Φ ∂X∂Y , (11)
where the Airy stress function Φ satisfies the elliptic problem

∆∆Φ = 1 in Ω, Φ = 0, ∂Φ/∂n = 0 on ∂Ω. ( 12
)

Estimation of the parameter ψ

Clearly the evaluation of the integral [START_REF] Lewicka | The foppl-von karman equations for plates with incompatible strains[END_REF] and the solution of problem [START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF] are the crucial, and possibly more difficult, steps in the reduction procedure. However it turns out that ψ can be estimated with good accuracy by simply applying a correction coefficient to the analytical expression [START_REF] Modes | Gaussian curvature from flat elastica sheets[END_REF] available for elliptical planforms.

For Ω an elliptical shape with axes ratio a/b and a membranally orthotropic material:

A =   1 ν 0 ν β 0 0 0 γ   , γ = ρ(1 -ν 2 /β), (13)
the problem ( 12) is solved in closed form, see [START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF], and ψ is evaluated to be:

ψ ell = β -ν 2 24π 2 (3(a/b) 2 β + (β/ρ -2ν) + 3(b/a) 2 ) . ( 14
)
For isotropic materials (β = 1, 2ρ = 1/(1 + ν)) and circular shapes (a/b = 1) this expression reduces to the simple:

ψ ic = 1 -ν 2 192 π 2 5.3 * 10 -4 (1 -ν 2 ). (15)
In Figure 4, we have plotted as a solid black curve the ratio ψ ell /ψ ic in the isotropic case in order to elucidate the dependence of ψ on the planform aspect ratio a/b: it is evident that for nearly circular planforms, the membrane energy has the maximal relevance with respect to the bending energy.

Ψ i Ψ ic
Fig. 4 The ratio ψ i /ψ ic for isotropic material as function of the planform aspect ratio a/b: exact solution for ellipses (black), FE solution (points) and the approximation (16) (dashed) for rectangular planform.

For rectangular shapes we claim that multiplying the closed form solution ( 14) by a correction factor is sufficient to get a quite accurate estimation; in particular we state that

ψ rect 0.738 * ψ ell . ( 16
)
As a matter of fact, in Figure 4, we have also reported as black points the values ψ rect /ψ ic , computed numerically solving [START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF] on rectangles with different aspect ratios and the function 0.738 * ψ ell /ψ ic as a gray dashed curve. One can check that the error associated to the estimation ( 16) is quite small and specifically is below 5% for a/b ∈ (0.5, 2). The number 0.738 has been numerically obtained as the ratio between the values of ψ for a circle and a square of isotropic material, but similar error bounds still hold for non isotropic materials.

More complex ansatz

Displacement ansatzs more complicated than (7) can be introduced to satisfy the boundary conditions for the bending moments and improve the accuracy of the UC solution. In particular we will refer in Section 6 to some results related to the Quadratic Curvature (QC) model arising from the ansatz

w(x, y) = kx x 2 2 + ky y 2 2 + kxy x y+ c 1 x 2 x 2 -a 2 + c 2 y 2 y 2 -b 2 . ( 17
)
This leads to a quadratically varying curvature field able to satisfy on average the boundary conditions for the bending moment in a rectangular domain. Again the interested reader is addressed to [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF] for details.

The analysis of the multistable behaviour of shells based on the QC model produces results which are very close to the ones obtained by the reduced UC model (see Figures 123).

Validation of UC model by comparison with Finite-Element (FE) simulations

In order to validate the analytical results for the UC model, we performed fully non-linear FE simulations: all the FE results were obtained with the Abaqus S4R element in a general non-linear setting by means of a Full-Newton scheme. FE results are plotted in Fig- ures 1-3 (details on how the stability boundaries are traced can be found in the Appendix) and show very good agreement with the analytical predictions: all the boundaries of the monostability region fall within a 3% error.

The FE model is not bound to the assumptions of shallow shells, as it is the case in the Föppl-von Kármán model, or its UC and QC approximations. Yet the overall stability scenario (Figure 3), as well as the main effects leading to multistable behaviour of shells (Figures 12), are confirmed by FE simulations even for large curvatures. Further comments on the differencies between analytical and FE results can be found in Section 6.

Based on these remarks, in the following we present design rules for multistable shells based on the UC model. These rules or recipes are applied in order to build examples of multistable shells (see Section 5).

Equilibria and stability properties

Aimed at giving simple design directives, we derive closed form solutions regarding the existence and stability of the UC equilibria. We focus on the case of orthotropic (particularly square-symmetric/isotropic) shells.

For square-symmetric shells the normalized stiffness matrix D can always be put in the form:

D :=   1 µ 0 µ 1 0 0 0 γ   , -1 < µ < 1, γ > 0, (18)
where µ := ν/ √ β and γ := ρ(1 -ν 2 /β)/ √ β are the only two constitutive parameters. Here ν is the in-plane Poisson's ratio, β := E y /E x the ratio between Young moduli in the orthotropy directions and ρ the shear modulus 2 . Isotropic shells correspond to β = 1 and γ = (1 -µ)/2.

Assuming ( 18) and a spherical initial natural curvature H 0 represented by the only parameter c 0 such that:

H 0x = H 0y = c 0 and H 0xy = 0 (21)
the stationarity equations (∂U/∂K = 0) give, after some algebra, the following equilibrium equations:

(K x -K y ) detK -c 2 0 + µ -1 = 0, (detK + µ + 1)∆ -c 2 0 ∆ -2(c i + c 0 ) (µ + 1) = 0, K xy detK -c 2 0 -2γ = 0, (22)
where

∆ := 4detK + 4K 2 xy + (K x -K y ) 2 .
Stability of equilibria is related to the sign of the Hessian ∂ 2 U/∂K 2 , and the stability margin is then represented by the condition:

∂ 2 U ∂K 2 = 0. ( 23
)
By solving [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF] 1 and ( 22) 3 , we see that three different types of equilibria are possible which are briefly discussed in what follows. In each of the cases discussed hereafter, the values of K x and K y at equilibrium are obtained by solving (22) 2 .

1. K x = K y and K xy = 0. These are "untwisted" configurations and they are spherical when β = 1. These equilibria can be either one or three. If we combine the stationarity condition [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF] 2 with the stability margin [START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF], in this case one can find the expression f I :

f I = 4c 6 0 -54c 0 c i (µ + 1) 2 -(µ + 1) 2 27c 2 i + 4µ + 4 + -12c 4 0 (µ + 1) -15c 2 0 (µ + 1) 2 ,
2 As shown in [START_REF] Fernandes | Multiparameter actuation for shape control of bistable composite plates[END_REF], a differential scaling of the curvature components can always be introduced

Kx → Kx β, Ky → Ky, Kxy → Kxy 4 β, (19)
Hx → Hx β, Hy → Hy, Hxy → Hxy 4 β,

which reduces the generic orthotropic material to the analysis of a square-symmetric material with β = 1.

(24)

it is easily seen that in the region f I < 0 there is only one stable equilibrium. Instead when f I > 0 there are three equilibria only two of which are stable, while on the set f I = 0 one of them looses stability as a turning point. Hence f I > 0 is the gray region in Figure 1 where everted configurations are stable and curve f I = 0 represents the boundary between monostable and bistable regions. 2. detK = c 2 0 +1-µ > 0 and K xy = 0. These equilibria correspond to "untwisted" configurations with positive assigned Gaussian curvature. If we combine the stationarity condition [START_REF] Vidoli | Discrete approximations of the Foppl-von Karman shell model: From coarse to more refined models[END_REF] 2 with the stability margin [START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF], in this case one can find the expression f II :

f II = (µ (c i + c 0 ) + c i -c 0 ) (µ (c i + c 0 ) + c i + 3c 0 )+4(µ-1), (25)
it is easily seen that these equilibria are stable when f II > 0. They originate from type I equilibria through a pitchfork bifurcation on the set f II = 0. Hence f II > 0 is the gray region in Figure 2 which allows for two bifurcated configurations, and the curve f II = 0 represents the boundary between monostable and bistable regions. 3. K x = K y and detK = c 2 0 + 2 γ. These equilibria correspond to "twisted" configurations with assigned Gaussian curvature. However it easily checked that these equilibria can be stable only if 0 < γ < (1µ)/2 ≤ 1 i.e. for non-isotropic shells with low shear rigidity.

If γ = (1 -µ)/2, as, for instance, in isotropic shells, then

detK = c 2 0 + 1 -µ ≡ c 2 0 + 2 γ, (26)
and two equilibrium equations are contemporarily solved by this value of the Gaussian curvature. Therefore when γ = (1 -µ)/2 the twisting curvature K xy and the difference (K x -K y ) can be chosen arbitrarily provided:

4 K 2 xy + (K x -K y ) 2 = (c i + c 0 ) 2 (µ + 1) 2 -4 (c 2 0 + 1 -µ). (27)
We then have a whole one-dimensional set of configurations, (27), sharing the same level of energy. The corresponding stability margin is vanishing and similar shells actually behave as zero-stiffness structures [START_REF] Seffen | Prestressed morphing bistable and neutrally stable shells[END_REF][START_REF] Guest | A zero-stiffness elastic shell structure[END_REF]. We notice that these equilibria can be seen as a degeneration of the case II described above and that a similar result was recently found for cylindrical shells in [START_REF] Guest | A zero-stiffness elastic shell structure[END_REF].

Useful analytical results for design purposes

Materials verifying the condition γ < (1 -µ)/2 are rare. We focus here on the case of the two basic classes described in Section 2 without considering twisting instability. In the important limit case of isotropic shells, twist becomes a degenerate zero-stiffness deformation mode. Having found the analytical boundaries, f I = 0 and f II = 0, of the monostability, bistability and tristability regions, one can derive precise closed form indications for the design of multistable shells. We focus in particular on the set c 0 > 0 and c i > 0 which is shown in Figure 5. With reference to this figure we first define the two critical values

c * 0 = 2 1 + µ, c * i = 2 √ 1 -µ 1 + µ , (28)
which are obtained intersecting f I = 0 with c i = 0 and f II = 0 with c 0 = 0 respectively. The value c * 0 gives the critical stress-free curvature over which a shell is bistable (i.e has an everted shape) in absence of inelastic curvatures; the value c * i gives the critical inelastic curvature over which a shell is bistable (i.e bifurcates its equilibrium) in absence of natural curvatures.

The two sets f I = 0 and f II = 0 intersect at the apex of the tristability region, namely in the point

c 0 = c * * 0 := 7 + µ, c i = c * * i := 4 √ 2 1 + µ -7 + µ. (29)
In the same point we compute the tangents to the two stability boundaries (gray vectors in Figure 5) in order to be able to move inside the tristability region. In particular we find that any point in the form

{c * * 0 , c * * i } + λ {cos φ, sin φ}, λ > 0, (30)
with

φ = arctan 3 √ 7 + µ + p √ 7 + µ - √ 2(1 + µ) √ 2 p (1 + µ) (31)
and 0 < p < 1, lies inside the tristability region. For p 1 the point (30) stays close to the lower stability boundary, while for p 1 is closer to the upper boundary. The point T, for instance, corresponds to p 0.24 and λ 3.4. Finally in Figure 6 we examine the dependance of the introduced critical values for c 0 and c i as functions of the material parameter µ. Whilst for c * 0 and c * * 0 we observe a moderate dependence, both c * i and c * * i tends to zero when µ → 1. This fact is of great practical importance as even small inelastic loads can bifurcate the equilibrium in materials having the parameter µ close to 1. Recalling the definition of µ, the condition µ → 1 is obtained when

ν 2 → E y /E x , (32)
a remarkable material property already found and discussed in [START_REF] Vidoli | Tristability of thin orthotropic shells with uniform initial curvature[END_REF] or [START_REF] Giomi | Multi-stability of free spontaneously curved anisotropic strips[END_REF].

Design recipes

Based on the results obtained above, we detail, in the following subsections, three recipes to actually build multistable shells in each of the regions shown in Figure 3:

(B) Bistable shells with natural and everted configurations (f I > 0, c 0 > 0); (C) Bistable shells with bifurcated quasi-cylindrical configurations (f II > 0, c i > 0); (T) Tristable shells with bifurcated quasi-cilindrical configurations and everted configuration (f

I > 0, f II > 0, c 0 > 0, c i > 0).
Preliminarly to every recipe there is the choice of both the planform shape Ω and the material. Concerning the shape, the best results are obtained with Ω a circle because it minimizes the critical values of curvatures for multistability. In fact, the dimensional values of curvatures are inversely proportional to the radius R (see [START_REF] Koehl | How kelp produce blade shapes suited to different flow regimes: a new wrinkle[END_REF]), which on its turn is proportional to √ ψ (the number ψ is defined in ( 11)): thus, the bigger is ψ, the bigger is R and the lower the corresponding dimensional curvatures. It can be seen from the results in Figure 4 that the maximum value for ψ are obtained for the circular shape. Concerning the material, in the orthotropic case we suppose given the Poisson's ratio ν, the Young modulus E x , the Young moduli ratio β = E y /E x and the shear modulus ρ. The best results are obtained for materials having sufficiently high shear modulus 2ρ ≥ 1 + ν (in order to eliminate twisting intabilities: see explanations in Section 3) and satisfying (32), in order to maximize the stability margins.

Steps common to all recipes

(I) Depending on the chosen shape Ω and on the material compute ψ:

using ( 14) if the shape is elliptical, using ( 16) if the shape is rectangular; either numerically solving [START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF] and [START_REF] Lewicka | The foppl-von karman equations for plates with incompatible strains[END_REF] or approximating ψ as in the previous cases if the shape is neither elliptical nor rectangular. (II) Choose the area of the planform L 2 and the thickness t. Note that the ratio L 2 /t is a scaling factor for the characteristic radius; hence, the larger is the number L 2 /t, the smaller curvature will be needed for multistability. (III) Compute the characteristic radius R in (10).

Bistable shells with natural and everted configurations

(IV) Compute c * 0 in (28). (V) Build a shell having, in a stress-free configuration, curvatures

h0x = β h * 0 , h0y = h * 0 , h0xy = 0, h * 0 ≥ 1.05•c * 0 /R. (33)
As, with respect to refined FE simulations, the UC model typically overestimates within a 3% error the critical natural curvature, the coefficient 1.05 in (33) guarantees the existence of a stable everted configuration. Clearly the value in (33) 3 represents the minimum curvature, but larger values are also possible. Note, however, that the energy needed for the shell eversion increases rapidly with the distance from the critical value, since it scales as

(det h0 -(h * 0 ) 2 ) 2 .
4.3 Bistable shells with two bifurcated configurations (IV) Compute c * i in (28). (V) Build a flat shell and apply to it inelastic curvatures

hix = β h * i , hiy = h * i , hixy = 0, h * i ≥ c * i /R. ( 34
)
Since, with respect to refined FE simulation, the UC model estimates very well the critical inelastic curvature, condition (34) should guarantee the existence of two bifurcated configurations without any additional coefficient. Clearly the value (34) 3 represents the minimum inelastic curvature to induce bifurcation of the equilibrium. Larger values are possible and needed to increase the difference between the bifurcated shapes; however, in these shells the energy, needed to move from one equilibrium to the other, increases less rapidly with the distance from the minimal critical value as it scales as ( hi -h * i ) 2 .

4.4 Tristable shells with two bifurcated and one everted configurations (IV) Compute c * * 0 and c * * i in (29). (V) Build a shell with stress-free curvatures

h0x = β h * * 0 , h0y = h * * 0 , h0xy = 0, (35)
and apply to it inelastic curvatures

hix = β h * * i , hiy = h * * i , hixy = 0, (36)
with

{h * * 0 , h * * i } = {c * * 0 , c * * i }/R + λ {cos φ, sin φ}/R (37)
and φ given in (31). A word of caution is needed in order to properly choose the parameters λ and p. Since the UC model underestimates the apex of tristable region by 14% in norm, the condition λ > 0 should be substituted by a safer λ > 1. Moreover the angle φ (controlled by 0 < p < 1) is better chosen to stay closer at the bottom boundary of the tristability region; safer values for p are within 0.05 < p < 0.5. This last choice allows to avoid local instabilities that the UC model is not able to describe and which are discussed in Section 6.

Sample prototypes

We present an experimental example for each of the three cases examined in the previous Section and corresponding to points B, T and C in Figs. 3 and5.

In particular the prototype corresponding to point B in Figure 5 is half a tennis ball, while the prototypes corresponding to points C and T are made of silicone. These two families of materials are parametrized by µ = 0.5, which is the only parameter affecting the shape of the stability diagram of Figure 3; hence, this diagram is valid for all cases. We remark that this value of Poisson's ratio is the highest possible for isotropic materials, thus the one which is closest to the condition (32).

Bistable spherical shell: half tennis-ball

A tennis ball is a spherical shell essentially made of rubber. One half of such a ball constitutes a nice example of bistable shell with both natural and everted configurations as shown in Figure 7. Modelling it as a shell with circular planform of radius r 0.03 m, thickness t 0.003 m, and using an isotropic rubber material with β = 1, ν = µ 0.5 and γ 0.25, we calculate:

ψ = 3.96 * 10 -4 , L = √ 2 πr 0.075 m, L/t 25, R = 0.13 m, c * 0 = 2.45, h * 0 = c * 0 /R = 18
.84 m -1 . Since β = 1 and the stress-free curvature of half tennis ball, namely h0 = 1/r = 33.33 m -1 , is sensibly larger than the threshold h * 0 , the everted configuration is stable, Figure 7. Clearly half a tennis ball is not a shallow shell; yet the UC model predicts the correct stability properties.

Bistable bifurcated shell and tristable shell: silicones

Whilst we presented the very simple case of half a tennis ball to illustrate the case of the bistable shell in B, the points C and T in Figure 5 do require the presence of inelastic curvatures: fabricating shells with large initial and inelastic curvatures is not a trivial task from the experimental point of view. The solution we adopted here is to use the silicones Zhermack Elite double 32 (referred as Green) and Zhermack Elite double 8 (referred as Pink). These soft and isotropic materials are easily moldable and can exhibit controlled inelastic deformations if suitably prepared; they are delivered in the form of liquid polymers which, once mixed with a catalyst, undergo a rapid polymerisation reaction.

The Green and Pink silicones have similar mechanical properties, namely ν 0.5, β = 1, ρ = 1/3 (leading to µ = 0.5 and γ = 0.25), but are characterized by different Young moduli (E Green = 1, 364 MPa, E Pink = 0, 186 MPa) and by different hygroscopic properties. We use this circumstance to produce shells with controlled inelastic effects: indeed we observed that the polymerisation of a layer of Pink silicone, which is laid over a layer of already polymerized Green silicone, produces a bilayered medium with inelastic curvature (an effect which is not documented in the technical material specifications). Moreover, this effect is also obtained when the top layer is made of a mixture of Pink and Green silicone, and the resulting inelastic curvature can be controlled by tuning the relative percentage of Pink and Green silicone used in the mixture. In particular, we suppose the bottom layer being made entirely of Green silicone and we characterize the mixture of the upper layer by its percentage ρ of Pink silicone (ρ ∈ [0, 100], the percentage of Green silicone being then (100 -ρ)). We measured the resulting inelastic curvature h i experimentally on 2mm-thick two-layered strips corresponding to different values of ρ: results are shown in Figure 8. We note that in such assemblies also a small percentage of Pink silicone in the upper layer suffices to produce an inelastic curvature higher than the critical threshold h * i . For the material properties of these silicones and a circular planform the preliminary steps described in Section 4. Entering with this last value in Figure 8b dictates a Pink percentage ρ in the upper layer higher than 5%; we choose in particular ρ 40% which corresponds to the inelastic curvatures hix = hiy 13 m -1 and to point C in Figure 3. The resulting bifurcated configurations are shown in Figure 9. The choice made for λ and p allows to have the same inelastic curvature of the previous case and, therefore, it is obtained as before with a mixture of 40% Pink and 60% Green in the upper layer. But now the bilayer plate must not be flat and must have curvature radius 1/h * * 0 5.5 cm equal in all directions since β = 1. To this end, using a 3D printer, we have produced the mold shown in Figure 11, and the polymerization phases for both the Green and Pink-Green layers has taken place inside this mold. The three stable shapes of the resulting tristable shell are shown in Figure 10.

Concluding remarks

We have reported on modelling, numerical approximation and design recipes of simple multistable shells. Using soft isotropic materials, we have built sample prototypes for three relevant cases. The following conclusions can be drawn:

1. All the FE results have been obtained with fully non-linear simulations of shells in Abaqus (S4R element, Full-Newton scheme: see the Appendix for details). Even if the Föppl-von Kármán model, or its UC and QC approximations, are based on the assumptions of thin (t L) shallow (k 1/L) shells, yet the UC and QC models are able to give an accurate description of the stability scenario (Figure 3) even for large curvatures. 2. The UC model gives faithful indications of the stability, associated to rather large natural and inelastic curvatures, despite the simplicity of the Uniform Curvature assumption. In particular the predictions on all the boundaries of the monostability region fall within a 3% error. We note that this remains true even when the actual (or FE computed) shapes violate the UC assumption. Figure 12a shows, for instance, the UC shapes (yellow) of a square tristable shell as compared with the ones predicted by the QC model or by FE simulations (brown). For the everted configuration, Figure 12b shows the variations of the FE curvature fields with respect to their averages. 3. The QC model, based on a quadratic ansatz of the curvature field, is slightly more precise than the UC model for predicting the number of stable equilibria (see, for instance, the lower boundary of the tristability region in Figure 3) and the actual equilibrium shapes (Figure 12a). However, in a possible tradeoff between these QC precision enhancements and the simple closed-form solutions allowed by the UC model, the latter should, in our opinion, be preferred because it is simpler but it suffices to well describe the multistability phenomena both qualitatetively and quantitatively. 4. Still with reference to Figure 3 the more pronounced difference between the UC (or QC) predictions and the FE results lies in the upper stability boundaries of the tristable regions. In Figure 13a we have plotted the average shell curvatures K x and K y as pre-dicted by UC (gray), QC (brown) and FE (black) for c 0 = 4.4. Clearly in the FE simulation the everted configuration looses stability for c i 4.6, i.e. much earlier than the UC and QC predictions. Plotting, in Figure 13b, the FE spatial distribution of the component K xy as c i 4.6 is approached, we understand that the stability is lost for a localization of the twisting curvature; a bifurcation mode that evidently could not be contemplated under the simple hypotheses of a uniform or quadratically varying curvature.

Appendix: Identification of stability boundaries for multistable shells by FE simulations

In Figure 3, we resume the prediction of the stability boundaries for multistable shells: results obtained by applying the UC and QC models (solid and dashed lines, respectively) can be compared to the ones computed by refined Finite Element approximations (marked as black dots). In particular, we performed FE simulations with Abaqus CAE using the S4R element, which is a shell element with four nodes and five degrees of freedom per node with reduced integration. Geometrical non linearities are taken into account in a Full-Newton scheme. In order to validate the analytical results as well as to locate the experimental response of the prototypes described in Section 5, results of simulations correspond to the same material system as in experiments, that is 2mm-thick bilayered shells made of Green-(Green+Pink) silicone material. The non-dimensional diagram of Figure 3 is valid for both silicone and rubber materials since they have the same Poisson's ratio (µ = 0.5) which is the only parameter affecting its boundaries. The planform of the shell is circular, but simulations based on square planforms were also performed and gave similar results. The diameter was initially fixed at D = 10 cm in order to reproduce the experimental results. Simulations were also performed with a diameter of D = 20 cm in order to apply sufficiently high initial curvatures c 0 (points with c 0 > 5 in Figure 3).

In the analytical models, the shell is free and no boundary condition is applied; thus, in the FE simulations we only clamped the central node in order to avoid rigid body motions.

In order to reproduce the stability diagram of Figure 3, the two parameters to be varied in the FE study are the natural curvature c 0 and the inelastic curvature c i (we remind that c 0 and c i are adimensional curvatures, according to definitions given in Section 3). Different levels of natural curvature c 0 simply correspond to different paraboloid shell surfaces c0 R (x 2 + y 2 ), where R is the characteristic radius defined by equation [START_REF] Koehl | How kelp produce blade shapes suited to different flow regimes: a new wrinkle[END_REF]: changing c 0 required to build a new mesh for the initial stress-free configuration, the number of elements of the mesh varying from 1456 (9114 d.o.f.) for a disk (c 0 = 0) to 1828 (11346 d.o.f.) for the deepest shell considered (c 0 = 9). Different levels of inelastic curvature c i were induced by applying an equivalent uniform temperature change ∆T to the shell (different expansion coefficients are defined for the two layers, thus inducing an inelastic curvature of the shell: values of the ficticious expansion coefficients are identified from the experimental behaviour of the silicone material system, shown in Figure 8).

From an operational point of view, we fixed a discrete set of c 0 values (see dots on Figure 3). For each value of c 0 , we built the corresponding mesh (stress-free configuration) of the shell and we found the corresponding everted configuration at c i = 0, if it existed, by applying imposed displacements onto the shell: clearly there is a range of values of c 0 (c 0 ∈ [-c * 0 , c * 0 ]) where the everted configuration cannot be found. Then, starting from each natural configuration and from each everted configuration, we have progressively increased the inelastic curvature value c i (in practice, the temperature change ∆T ) until the stability was lost, thus marking the critical values of inelastic curvature in order to obtain bifurcation or loss of equilibrium.

In order to find the different points on the stability diagram (see Figure 14), we followed the procedures described hereafter. 3. Dashed gray boundary (c 0 > 0): start from a spherical shell with initial curvature c 0 > 0 and apply a positive inelastic curvature c i until finding the value corresponding to the pitchfork bifurcation of the natural configuration into two quasi-cylindrical shapes.

4. Dashed black boundary (c 0 > 0): first find the everted configuration, then apply a positive inelastic curvature c i until it looses stability (above the dashed black line the everted configuration is no longer stable: tristability is attested above the dashed gray line and below the dashed black line).

5. Solid black boundary (c 0 < 0): start from a curved shell with initial negative c 0 curvature, and apply a positive inelastic curvature c i until it looses stability.

6. Solid gray boundary (c 0 < 0): apply a positive inelastic curvature c i onto an everted configuration corresponding to intial negative natural curvature c 0 and increase c i until the everted configuration bifurcates into two quasi-cylindrical shapes.

Fig. 1

 1 Fig. 1 Average curvatures at the equilibrium when varying the initial curvature c 0 with vanishing inelastic curvature c i = 0:

Fig. 2

 2 Fig. 2 Average curvatures at the equilibrium when varying the inelastic curvature c i with vanishing initial curvature (c 0 = 0):

Fig- ure 3

 3 shows the regions of monostability (white), bistability (light gray) and tristability (darker gray) in the plane (c 0 , c i ), obtained by combining curves from Figures 1 and 2.

Fig. 3

 3 Fig. 3 Stability diagram in the half-plane c 0 , c i > 0. The gray level refers to the number of stable equilibria: 1 (white), 2 (light gray) or 3 (darker gray). Black points report the Abaqus FE results, dashed and solid curves report respectively the stability boundaries found with the UC and QC reduced models.

Fig. 5 Figure 3 .

 53 Fig. 5 Zoom for c 0 > 0, c i > 0 of the stability diagram in Figure 3. Definitions of the limits c * 0 , c * i , c * * 0 and c * * i .

Fig. 6

 6 Fig. 6 Critical values of the curvatures as functions of µ: c * 0

Fig. 7

 7 Fig. 7 Natural (a) and everted (b) configurations for half a tennis ball.

1 )Fig. 8

 18 Fig.8Inelastic curvature h i of bi-layered Pink-Green silicone strips induced by the percentage ρ of Pink material in the Pink-Green mixture of the upper layer. The critical value h * i , needed in the sample case discussed, and the inelastic curvatures of points T and C are reported: these last correspond to a 40% Pink-60% Green mixture.

  1 give: (I) ψ = 3.96 * 10 -4 . (II) The planform radius and thickness are chosen respectively 5 cm and 2 mm. Hence L 2 /t = 3.93 m. (III) R = 0.271 m. To build the bistable shell with quasi-cylindrical configurations, we follow the steps in Section 4.3 to get: (IV) c * i = 0.943. (V) h * i > 3.48 m -1 .

Fig. 9

 9 Fig. 9 Bifurcated shapes of the silicone shell corresponding to point C in Figure3. The principal directions of curvature are marked on the shell in red and blue.

Finally, to

  build the tristable shell we follow the steps in Section 4.4 to get: (IV) c * * 0 = 2.74 and c * * i = 1.03 (V) The choice λ = 3.4 and p = 0.24 leads to h * * 0 18.1 m -1 and h * * i 13 m -1

Fig. 10

 10 Fig. 10 Shapes of the silicone shell corresponding to point T in Figure3. The principal directions of curvature are marked on the shell in red and blue.

Fig. 11

 11 Fig. 11 Mold with radius of curvature 5.5 cm used to build the tristable shell in Figure 10.

Fig. 12 a

 12 Fig. 12 a) Upper row: Shapes predicted by UC (yellow) and QC (brown) models for a square shell within the tristable region c 0 > 0, c i > 0; actually, the FE predictions are almost confused with the QC ones and are not reported exactly. b) Lower row: Variations of the curvatures computed with FE with respect to their spatial averages: from the left components Kx -Kx, Ky -Ky and Kxy -Kxy in the everted configuration of a square shell.

Fig. 13 a

 13 Fig. 13 a) Comparison of the average shell curvatures for c 0 = 4.4 and c i increasing: UC (gray), QC (brown), FE (black). b) The Kxy field when c i reaches the critical value making the everted configuration unstable as computed in Abaqus.

1 .

 1 Intersection between solid gray and dashed gray lines (black point, c * i ): start from a flat disk and increase the inelastic curvature c i until finding the critical value c * i inducing bistability through a pitchfork bifurcation (evolution of curvatures is shown in Figure 2.a).2. Intersection between dashed black line and the hor-izontal axis c i = 0 (gray point, c * 0 ): start from a spherical shallow shell and impose displacements onto the shell in order to evert the curvature; repeat the simulations by increasing the initial curvature c 0 until the existence of the everted configuration is established at c 0 = c * 0 (see results on Figure1.a).

Fig. 14

 14 Fig. 14 Multistability boundaries for all materials parametrized by µ = 0.5.

Acknowledgements The authors acknowledge the financial support of Project ANR-13-JS09-0009 (Agence Nationale de la Recherche, 2014).


Alexandre Radjesvarane

Alain Damlamian

Tatsien Li

François Murat

dedicated to Professor Luc Tartar

HAL is


Claire Lestringant

Corrado Maurini

Arnaud Lazarus

Basile Audoly

Buckling of an elastic ridge: competition between wrinkles and creases

Keywords: numbers: 46.32.+x, 46.70.De, 83.80.Va

come

In extended systems, elastic instabilities generally produce smooth patterns having a well-defined wavelength. There are numerous examples involving an elastic beam [1] or a thin film [2][3][4] on an elastic foundation, a bulk elastic material with inhomogeneous elastic properties [START_REF] Lee | Thouless[END_REF], or rod-like solids with large incompatible strain [6][7][8], with applications ranging from morphogenesis [1] to the active control of surface properties [9]. An important exception to this rule is when the bifurcation problem has no intrinsic length-scale as happens for a compressed hyperelastic block, a problem considered by Biot [10,[START_REF] Biot | Mechanics of Incremental Deformations[END_REF]: a continuum of linear modes appear simultaneously at the bifurcation threshold with all possible wavelengths. This free-surface instability has been characterized numerically and experimentally only recently, and was found to be subcritical, localized and non-linear in essence [START_REF] Hong | [END_REF][13][14][15][16][17][18][19][20]. In spite of recent progress [21,22], there is no simple and systematic theoretical argument that explains why and in which circumstances localized creasing patterns are to be observed, nor whether scaleinvariance is a sufficient condition for localization.

Here, we analyze a variant of Biot's compressed elastic block, in which we replace the half-space geometry by a prism. The ridge angle φ brings in an additional parameter. Experimentally, we find buckling patterns reminiscent of creasing when the prism is flat enough (φ close to 180 • ), consistent with prior work [START_REF] Hong | [END_REF][13][14][15][16][17][18][19][20]. For acute enough ridge angles, however, a smooth buckling mode develops near the ridge, with a well-defined wavelength. We carry out a linear stability analysis of the compressed hyperelastic prism and investigate the competition between smooth and localizing buckling modes.

In our experiments, we use an isosceles-triangular prism made of a silicon elastomer (Ecoflex). This elas- tomer is nearly incompressible with a Young's modulus E p ≈ 0.06 ± 0.02MPa. Its lower face is bonded to a parallepipedic silicone block made of vinylpolysiloxane whose Young's modulus is ∼ 20 times larger, E b ≈ 1.3±0.05MPa. Both the prism and the base are obtained by casting liquid polymer into molds made of PMMA obtained by laser-cutting. We stretch the base to a length L 0 prior to glueing the prism onto it, see Fig. 1. By bringing the ends of the base closer to one another we induce a compressive axial strain = L0-L L0 in the prism, that depends on the current length L < L 0 of the base. At a critical value of the strain c an instability is observed which is localized along the free ridge of the prism, opposite to the base. Using different molds we repeat the buckling experiment for different ridge angles φ in the range 20 • to 120 • . The height h of the triangular prism is chosen at least ten times smaller than L 0 , so we can ignore finite-length effects in the analysis.

For φ smaller than a critical value φ ≈ 90 • we observe a smooth, extended buckling mode whereby the ridge bends out of the plane of symmetry of the prism, see Fig. 2a; we will refer to this as antisymmetric wrinkling (AW). For a given angle φ ≤ φ the wavelength λ scales close to linearly with the height of the prism h for the range of heights tested in the experiments, see the inset of Fig. 2.

For φ ≥ φ the buckling mode is entirely different, see Fig. 2c: localized creases are initiated at the ridge. As the strain is increased beyond c , more creases are formed and they spread along the lateral faces toward the base. The gap between successive creases does not appear to be regular. This buckling mode will be referred to as surface creasing (SC).

Overall, c increases steadily with the ridge angle φ until it reaches a plateau at c ≈ 0.42 for φ = φ , where the nature of the buckling mode changes, see Fig. 2d. This value is lower than the critical Biot strain Biot = 0.55 calculated by Biot [10,[START_REF] Biot | Mechanics of Incremental Deformations[END_REF] for the surface instability, as explained below.

We set up a bifurcation analysis, with the aim of characterizing the instabilities and of explaining the competition between the localized and extended buckling modes. The system is modeled as an infinitely long prism with triangular cross-section D made of an hyperelastic material. Its elastic energy density is denoted by

W 3D (E) where E = 1 2 (F T • F -1) is the strain tensor, F = ∂(x,y,z)
∂(X,Y,Z) the transformation gradient, (X, Y, Z) are the coordinates in reference configuration with Z aligned with the prism axis and Y along the axis of symmetry of the triangular cross-section D, and (x, y, z) are the coordinates in deformed configuration. The expression of W 3D (E) reflects the choice of a material law; we use a Gent model, as described in the supplementary material, with a choice of material parameters that makes this constitutive law practically equivalent to an incompressible neo-Hookean model. Working in the framework of finite elasticity, we denote by ϕ(X, Y, Z) = (x, y, z) -(X, Y, Z) the displacement. The non-linear equilibrium is obtained by the principle of virtual work as

∀ ϕ(X, Y, Z), L0 0 D Σ : F T • F dX dY dZ = 0, (1)
where Σ = ∂W 3D ∂E denotes the stress and F = ∂ ϕ ∂(X,Y,Z) the virtual increment of deformation gradient. As the average strain is imposed by the base, we consider only ad- missible virtual displacements ϕ whose incremental axial strain is zero on average. Taking advantage of the fact that the buckling patterns are localized near the ridge in the experiments, we simplify the boundary conditions at the interface with the base which we replace by a free boundary.

The unbuckled solution is in a state of homogeneous 'simple' compression as described by ϕ 0 = η( ) X e x + Y e y -Z e z . Here η( ) captures the dilation of the cross-section by Poisson's effect and is found from the constitutive law by solving ∂W 3D ∂η , η( ) = 0. We consider a small perturbation ϕ 1 , to this invariant solution ϕ 1 = ξ x (X, Y ) e x + ξ y (X, Y ) e y + i ξ z (X, Y ) e z e i q Z in the form of a pure Fourier mode with wavenumber q. The virtual displacement ϕ is sought in a similar form. Upon linearization and discretization using the finite-element method, the equation of equilibrium (1) takes the form

∀ ξ, ξ • K + q C + q 2 M • ξ 1 = 0, ( 2
)
where ξ 1 = (ξ x , ξ y , ξ z ) and ξ = ξ x , ξ y , ξ z are two vectors collecting the Fourier amplitudes of the real and virtual nodal displacements on the cross-section D. The Fourier analysis thus yields a 2-d eigenvalue problem in which the third dimension enters through the wavenumber q only. For details of this 2-d formulation and of its implementation, see Supplemental Material (Supp. Mat.) and [8]. To discretize and solve the eigenvalue problem, we make use of the finite-element library FEniCS [26] and of the SLEPc library [27]. Equation ( 2) is invariant when a homothety is applied to both the solution and the wavelength 2π/q; in addition, the domain D is scale-invariant near the tip (ridge). As a result, an infinite number of modes that are homothetic one to another appear concurrently at the critical strain c . These modes are localized near the ridge and are associated with all possible wavenumbers: there is no selection of the wavenumber in this scale-invariant linear bifurcation analysis, see Supp. Mat. for details. By contrast, the critical strain c and the shape of the buckling mode (up to a dilation) are selected as a function of φ.

As the unbuckled configuration is mirror-symmetric with respect to the (yz) plane, the buckling modes can be either symmetric or anti-symmetric. When φ is smaller than ≈ 105 • , the first critical buckling mode predicted by the FEM analysis is an anti-symmetric wrinkling mode (AW), see Fig. 4a. It involves lateral undulations of the ridge, see Fig. 3b, similar to the buckling mode seen in the experiments. The corresponding critical strain c is plotted in Fig. 3a (disks) and compared to experimental results (open circles): c (φ) is in good agreement with the experiments, and increases with φ.

In the limit of an acute ridge angle, φ → 0, the prism can be modeled as a thin, infinitely long plate whose thickness t(y) varies linearly with the distance to the ridge, t = φ |h -y|. For the unbuckled solution, the midsurface of the plate is contained in the (yz) plane and has an axial pre-stress σ 0 = E . When linearized about this solution, the Föppl-Von Kàrmàn equations for elastic plates yield, see for instance [START_REF] Audoly | Elasticity and geometry: from hair curls to the nonlinear response of shells[END_REF], (m αβ ) ,αβ + t(y)σ 0 w ,zz = 0 plate, E is Young's modulus and ν is Poisson's ratio. A comma in subscript denotes a partial derivative and Greek symbols are restricted to in-plane directions, α, β ∈ {y, z}. We use Einstein's convention for implicit summation on repeated indices.

In the plate model, we consider perturbations that are harmonic in the axial direction and rescale the vertical coordinate using the wavelength, w(y) = w(q y) e i q z . When expressed in terms of w and q, the boundary value problem (3) and the associated boundary conditions depend on the two dimensionless parameters ν and

σ 0 = 12 (1-ν 2 ) σ 0 E φ 2
. A numerical solution based on a shooting method yields the critical value σ 0 c (ν), see Supp. Mat. for details. The corresponding critical strain is

c = σ 0 c (ν)
12 (1-ν 2 ) φ 2 . For our particular 3-d constitutive law, ν = 0.45 and we obtain c (0.45) ≈ 0.35 φ 2 ; the dependence on Poisson's ratio is mild, c (0.50) ≈ 0.33 φ 2 in the incompressible case. This prediction has no adjustable parameter and is plotted in Fig. 3a (red dashed line): it agrees asymptotically with the finite element analysis for φ → 0. Note that c ∼ φ 2 is small when φ → 0, which is consistent with the linear elastic behavior assumed in the plate model.

Symmetric wrinkling modes (SW) are also found in the numerical bifurcation analysis. They extend on the adjacent faces on both sides of the ridge and involve an undulation of the ridge in the plane of symmetry (yz), see Fig. 4b. The strain at which the first symmetric mode appears is ≈ 0.55, a value which hardly depends on the ridge angle φ, see Fig. 4a. This value is consistent with the critical Biot strain Biot = 0.55 corresponding to the existence of a marginally stable surface mode in a pre-stressed neo-Hookean half-space [10,[START_REF] Biot | Mechanics of Incremental Deformations[END_REF]. This is consistent with the fact that the SW mode is localized just beneath the faces of prism, see Fig. 4b. When φ reaches ≈ 105 • , the critical strain c of the AW mode becomes larger than Biot = 0.55: the numerical analysis then predicts that the first buckling mode switches from an antisymmetric mode (AW) to a symmetric (SW) mode, see Fig. 4a.

This linear analysis therefore predicts a symmetric wrinkling mode (SW) which is smooth and sinusoidal, in apparent contradiction with the localized pattern observed in the experiments. When this mode becomes unstable, all wavelengths appear concurrently: it is known that the non-linear coupling between the different wavelengths gives rise to a creasing instability through a sub-critical bifurcation [14,21]. The buckling strain for the creasing instability in a neo-Hookean half-plane crease ≈ 0.44 is therefore lower than that predicted by the linear analysis Biot ≈ 0.55, see [14,17]. Extrapolating to our problem, this suggests that our SW modes are subcritical as well, and that the critical strain Biot predicted by the linear analysis needs to be corrected: the value crease has been included in Fig. 4a and indeed corresponds to the plateau observed in the experiments, see Fig. 4. Accordingly, the critical ridge angle φ can be found by equating the critical strain for antisymmetric modes c (φ) with the creasing strain crease : this yields φ = 88 • , see Fig. 4, which accurately matches the experimental value φ ≈ 90 • .

Our linear stability analysis correctly captures the dependence of the critical strain on the ridge angle, c (φ), as well as the shape of the antisymmetric mode. In our scale-free formulation, there is no selection of the wavelength. To account for the wavelength of the antisym-metric mode, one would need to consider additional ingredients in the analysis, such as subtle nonlinear effects and/or small-scale regularization. By contrast with the antisymmetric mode, the symmetric mode predicted by the linear stability analysis is not observed, as it gives rise to creasing by a subcritical bifurcation. Combining our linear analysis with the non-linear threshold for creasing, we have explained the critical value of the ridge angle φ ≈ 90 • at which the pattern changes. Interestingly, close to φ , the system displays a mix of the two behaviors: creases superimposed onto the smooth antisymmetric mode are shown in Fig. 2b, probably resulting from the non-linear interaction between the symmetric and antisymmetric modes.

The creasing localization has been explained in earlier work by non-linear coupling of the buckling modes. A remarkable finding of our experiments is that our system features both localized creases and a smooth extended buckling pattern: the coupling between modes of different wavelengths is effective for the symmetric mode (leading to creases) but it is not effective for the antisymmetric mode, surprisingly. Therefore, scale-invariance in not a sufficient condition for localization and the exact conditions in which modes of different wavelengths can cooperate remain to be elucidated: the compressed hyperelastic prism provides a workbench for future non-linear analyses of creasing.

FIG. 1 :

 1 FIG.1: Sketch of the experimental setup: to set the prism in axial compression, we first strech the substrate, then glue the prism to the substrate while keeping the substrate in tension, and finally release the substrate. This induces buckling of the prism: extended wrinkling (top row, φ = 40 • ) and localized creasing (bottom row, φ = 120 • ) are observed, depending on the value of φ. Insets: experimental pictures.

FIG. 2 :

 2 FIG. 2: Experimental results. (a-c) Top-views in the (x, z)plane. Each set of pictures is for increasing axial strain . A: φ = 40 • , h = 10 mm and L0 = 100 mm, B: φ = 90 • , h = 5 mm and L0 = 100 mm, and C: φ = 120 • , h = 10 mm and L0 = 100 mm. Black arrows highlight the creases visible on the faces. (d) Critical strain c(φ). The thin, hand-drawn curves reveal the trend of the experimental data points. Inset: rescaled wavelength of the extended mode λ/h for φ = 30 • .

( 3 )FIG. 3 :

 33 FIG. 3: Anti-symmetric wrinkling (AW). (a) Phase diagram c(φ) from experiments (open circles), simulations (disks) and analytical model (dashed curve). (b) Numerical buckling mode for φ = 40 • , h = 5 mm, shown with an arbitrary amplitude. The two colormaps show the amplitude of the lateral displacement (left) and of the incremental hoop strain E 1 θθ

FIG. 4 :

 4 FIG. 4: (a) Full bifurcation diagram, comparing the modes predicted by the linear bifurcation analysis (AW and SW), Biot's threshold, the non-linear creasing threshold, and experiments. (b) Numerical linear buckling mode for φ = 120 • , h = 0.5. Sketch of the deformed prism superimposed with a colormap of the amplitude of the vertical displacement ξy. Colormap of the amplitude of the hoop strain E 1 θθ . (c) Sketch of the experimental surface creasing mode (SC) for φ = 120 • . (d) Sketch of the sub-critical bifurcation curve A( ) for creasing.

Acknowledgments: We thank A. El Ouardy for his contribution to the experiments. CM acknowledges the financial support of Project ANR-13-JS09-0009 (Agence Nationale de la Recherche, 2014).