Claire Lestringant

Basile Audoly

Elastic rods with incompatible strain: macroscopic versus microscopic buckling

Keywords: B. Beams and columns, C. Stability and bifurcation, C. Asymptotic analysis

We consider the buckling of a long prismatic elastic solid under the combined effect of a pre-stress that is inhomogeneous in the cross-section, and of a prescribed displacement of its endpoints. A linear bifurcation analysis is carried out using different structural models (namely a double beam, a rectangular thin plate, and a hyper-elastic prismatic solid in 3-d): it yields the buckling mode and the wavenumber q c that are first encountered when the end-to-end displacement is progressively decreased with fixed pre-stress. For all three structural models, we find a transition from a long-wavelength (q c = 0) to a short-wavelength first buckling mode (q c = 0) when the inhomogeneous pre-stress is increased past a critical value. A method for calculating the critical inhomogeneous pre-stress is proposed based on a small-wavenumber expansion of the buckling mode. Overall, our findings explain the formation of multiple perversions in elastomer strips, as well as the large variations in the number of perversions as a function of pre-stress and cross-sectional geometry, as reported by Liu et al. (

Introduction

By pre-stressing a thin or slender elastic body, it is possible to prevent the propagation of cracks associated with tensile stress, as in the classical examples of tempered glass [START_REF] Aben | Photoelasticity of Glass[END_REF] and pre-stressed concrete [START_REF] Freyssinet | Exposé sur l'idée de précontrainte[END_REF]. When it is compressive, the pre-stress can lead to instabilities: the analysis of structures with incompatible strain, arising e.g. from growth or thermal effects, has applications ranging from morphogenesis [START_REF] Savin | On the growth and form of the gut[END_REF][START_REF] Osterfield | Three-dimensional epithelial morphogenesis in the developing Drosophila egg[END_REF] and the delamination of thin films [START_REF] Moon | Directed assembly of fluidic networks by buckle delamination of films on patterned substrates[END_REF] to the rolling and levelling of thin metal sheets [START_REF] Tomita | Buckling behavior in thin sheet metal subjected to nonuniform membrane-type deformation[END_REF][START_REF] Komori | Analysis of cross and vertical buckling in sheet metal rolling[END_REF][START_REF] Fischer | Buckling phenomena related to rolling and levelling of sheet metal[END_REF]Abdelkhalek et al., 2015). These instabilities fall in two classes: they are either macroscopic, meaning that the wavelength of the buckling mode is set by the size L of the structure, or microscopic, when the wavelength is instead set by the dimension h L of the cross-section. In recent experiments, a pre-stressed elastomer strip has been shown to display both types of behaviors [START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF]: depending on the geometry and of the pre-strain profile, wavelengths ranging from the microscopic scale h to the macroscopic scale L have been reported. The present work aims at explaining these unusually large variations of the buckling wavelength which, as we will see, are typical of slender elastic bodies subject to large and inhomogeneous pre-stress.

In the experiments of [START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF], two long elastomer strips I and II are used. The strip II is first first subjected to a pre-stretch p > 1, see figure 1.1(a-b); then the strips are glued and held by terminal forces, which are then released progressively, see figure 1.1(c-d). An instability takes place during unloading as the stress eventually becomes compressive in strip I. In the post-buckled regime, the shape of the strip is made up of pieces of helices having alternate chiralities [START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF], that are connected by localized defects termed 'perversions' in prior work [START_REF] Goriely | Spontaneous helix hand reversal and tendril perversion in climbing plants[END_REF][START_REF] Mcmillen | Tendril perversion in intrinsically curved rods[END_REF]. A striking feature of the experiment, which has not been explained to date, is that the number of perversions varies 2014). (a-c) Preparation of the bistrip with inhomogeneous pre-stress in the cross-section; (d) the terminal tensile load is decreased and the bistrip buckles. (e) Incipient buckling mode (top) and post-buckled solution (bottom), for macroscopic (left) and microscopic (right) buckling. The position of the perversions in the post-buckled solutions is determined by the crests and valleys of the incipient buckling pattern (star symbols): a macroscopic buckling mode yields at most one perversion in the entire strip (left) while a microscopic buckling mode yields two perversions per wavelength (right).

greatly as a function of the amount of pre-stretch p and of the aspect ratio of the cross-section, as documented by [START_REF] Liu | Structural transition from helices to hemihelices[END_REF]. For small pre-stretch p and a stubby cross-section (i.e. a large aspect-ratio h/v, where h is the thickness and v the width, see figure), there is at most one perversion over the entire length L of the strip and the instability is macroscopic; conversely, for large pre-stretch and a flat cross-section (small h/v), a large number of perversions is obtained, their spacing being comparable to the cross-sectional dimension h: the instability is then microscopic. Even though perversions appear progressively in the post-buckled range, we observe that the selection of the number of perversions is entirely determined by the wavelength of the incipient buckling mode, when the planar solution first becomes unstable. Indeed, the experimental results sketched in figure 1.1(e) show that the perversions form at the points where the deflection of the incipient buckling mode is maximal; upon releasing the endpoints further, the distance between successive perversions remains close to half the wavelength of the incipient buckling mode, in material coordinates. Accordingly, we address in this paper the wavelength selection of the first buckling mode, as sketched in figure 1.1(c-d). Our goal is to account for the large variations of the wavelength, from the macroscopic scale ∼ L down to the microscopic scale ∼ h depending on the amount of pre-stretch and of the aspect-ratio h/v. In their paper, [START_REF] Liu | Structural transition from helices to hemihelices[END_REF] successfully compare their experiments with finite-element simulations; we complement their approach by deriving analytical models based on the same set of mathematical equations as their simulations, thereby providing an interpretation of their results.

In view of the slender geometry of the bistrip, it is tempting to model it as a thin elastic rod. Thin rod models, however, do not correctly account for the wavelength selection in this particular system: they systematically favor long-wavelength instabilities and fail to predict the microscopic buckling mode seen in the experiments. This has not been fully appreciated, and a discussion of the selection of the number of perversions has remained elusive so far in the literature. [START_REF] Goriely | Spontaneous helix hand reversal and tendril perversion in climbing plants[END_REF]; [START_REF] Mcmillen | Tendril perversion in intrinsically curved rods[END_REF] have described an isolated perversion mathematically, by deriving a localized non-linear solution to the equilibrium equations for a naturally curved thin elastic rod. Using the same model of a naturally curved thin elastic rod, [START_REF] Domokos | Multiple helical perversions of finite, intrinsically curved rods[END_REF] have identified numerical solutions comprising multiple perversions but they report that these solutions live on a branch of equilibria that is not connected to the fundamental (unbuckled) branch: this confirms that thin rod models cannot account for the experimental sequence sketched in the right-hand side column of figure 1.1(e), where the origin of multiple perversions can be traced back to a microscopic buckling mode. Recently, [START_REF] Liu | Structural transition from helices to hemihelices[END_REF] have suggested that a microscopic wavelength can be selected by inertial effects using a thin rod model, but inertia is associated with time scales much shorter than the inverse loading rate of the experiments, and is therefore likely irrelevant. As we show in this paper, the classical Euler-Bernoulli rod model cannot account for the wavelength selection in the bistrip, as the microscopic bucking mode deforms the cross-sections in a way that is incompatible with the kinematic assumptions underlying the Euler-Bernoulli rod model.

Having ruled out the possibility to use an Euler-Bernoulli rod model, one might hope that an alternative, more sophisticated 1-d model may be identified in the vast literature concerned with the justification of rod and plate models from 3-d elasticity-see for instance the seminal works of Le [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Friesecke | A hierarchy of plate models derived from nonlinear elasticity by Gammaconvergence[END_REF] in this area. This is not the case: the experiments of [START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF] involve a finite pre-stress with significant variations across the cross-section, while available convergence results are limited to the case of a weak pre-stress (i.e. which vanishes for h → 0) and/or to the case of a pre-stress that is homogeneous through the thickness. The bistrip experiment therefore challenges known results on dimension reduction for elastic rods. The main results available in the literature can be summarized as follows. For thin elastic bodies, the convergence of the solutions of 3-d elasticity has been established towards variants of the classical plate, shell and membrane models in the limit h → 0. In the plate or shell models thus obtained, the pre-stress modifies the reference metric and/or the reference curvature [START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: A constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF][START_REF] Efrati | Elastic theory of unconstrained non-euclidean plates[END_REF][START_REF] Lewicka | The föppl-von kármán equations for plates with incompatible strains[END_REF][START_REF] Bhattacharya | Plates with incompatible prestrain[END_REF]; these thin plate or shell models have been applied to a variety of problem involving growth [START_REF] Marder | Theory of edges of leaves[END_REF][START_REF] Audoly | Self-similar structures near boundaries in strained systems[END_REF][START_REF] Klein | Shaping of Elastic Sheets by Prescription of Non-Euclidean Metrics[END_REF][START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Liang | Growth, geometry, and mechanics of a blooming lily[END_REF]. For slender elastic bodies, a convergence result has been established by [START_REF] Cicalese | On local and global minimizers of prestrained thin elastic rods[END_REF], whereby the equilibrium solutions for a long prismatic 3-d hyper-elastic solid are shown to converge asymptotically to the solutions of a thin elastic rod model with natural curvature: this result, which holds under the assumption of weak pre-stress, brings us back to the Euler-Bernoulli rod model, and therefore cannot account for microscopic buckling. In the absence of an applicable 1-d, we shall resort to 3-d finite elasticity theory to analyze the perversions in the bistrip experiment, see §5 of this paper.

This paper builds on important previous works. Euler buckling has been approached by [START_REF] Scherzinger | Asymptotic analysis of stability for prismatic solids under axial loads[END_REF] based on the 3-d bifurcation analysis of a prismatic hyper-elastic solid with an arbitrary cross-section: this instability is macroscopic, and the critical load was found to be asymptotically consistent with the Euler-Bernoulli rod model in the limit of a small aspect-ratio h/L; here, we extend this analysis to include the effect of pre-stress and arrive at a different conclusion. The pre-stressed strut on an elastic foundation [START_REF] Timoshenko | Theory of elastic stability, 2nd Edition[END_REF] is a textbook example of a microscopic instability, of which a number of variants and extensions have been discussed [START_REF] Lee | Surface instability of an elastic half space with material properties varying with depth[END_REF][START_REF] Savin | On the growth and form of the gut[END_REF]; here, we revisit this example, with the first strip I and the second strip II playing the roles of the strut and of the foundation, respectively: unlike a normal elastic foundation, the second strip II is not anchored to a fixed base and can be carried away by the strut, thereby opening the possibility of a long-wavelength instability as well, as we shall show. In the context of the rolling of thin metal sheets [START_REF] Fischer | Buckling phenomena related to rolling and levelling of sheet metal[END_REF]Abdelkhalek et al., 2015), buckling modes that are sometimes macroscopic and sometimes microscopic have appeared in numerical simulations of pre-stressed elastic plates [START_REF] Rammerstorfer | Buckling of free infinite strips under residual stresses and global tension[END_REF][START_REF] Kpogan | Simulation numérique de la planéité des tôles métalliques formées par laminage[END_REF]; here, we propose a systematic analysis of the macroscopic-to-microscopic transition, beyond the limited framework of elastic plates.

In this paper, the buckling of a pre-stressed prismatic solid is analyzed as follows. We carry out a linear bifurcation analysis of the cylindrically invariant solution as a mean to predict the distance between perversions, see figure 1.1(e): we seek the first critical wavenumber q c when the average imposed stretch λ is progressively decreased. We consider an infinitely long prismatic solid, L → ∞, with a cross-section having fixed dimensions h × v: in this framework, a macroscopic buckling mode corresponds to a vanishing first critical wavenumber, q c = 0, while a microscopic one corresponds to a non-zero wavenumber, q c = 0. The pre-stress distribution in the solid is assumed to be inhomogeneous in the plane (xy) of the cross-section, but independent of the axial coordinate z: with suitable assumptions on the material symmetry (see below), this allows us to use Fourier analysis with respect to z. The pre-stress is finite, i.e. it is not assumed to depend in any way on a small parameter-incidentally, note that the aspect-ratio parameter v/L is zero (and not a small parameter) as we work with L = ∞. We use different structural models for the solid, all of which allow the cross-section to deform: as mentioned earlier, this is a necessary condition for the existence of microscopic buckling modes.

Three models are considered in this paper. Starting at a qualitative level, a double-beam model is ). The mean imposed stretch λ is progressively decreased (meaning that the diagram has to be read from right to left) and the homogeneous solution bifurcates along the solid curve. (a) Case of a microscopic buckling instability: the system first bifurcates at λ = λc and a wavenumber qc = 0 is selected. We denote by λ (0) the stretch such that a linear mode with an infinite wavelength is available near the fundamental branch: λ (0) < λc for microscopic buckling. (b) Case of a macroscopic buckling instability: a zero wavenumber is selected, qc = 0, and λ (0) = λc.

considered in §3; its mathematical formulation is relatively straightforward and it successfully captures the salient features of the bistrip experiment, including the macroscopic-to-microscopic transition. Gradually increasing the level of the detail (and of technicality), we analyze a non-linear plate with inhomogeneous pre-stress in §4. Finally, we turn in §5 to a bifurcation analysis of a full 3-d, hyper-elastic model. In addition, an expansion method inspired by Koiter's weakly non-linear expansion is proposed to characterize the critical pre-stress the macroscopic-to-microscopic transition. This method is first derived in an abstract setting in §2, and then applied to the three structural models successively. The sections of the paper are largely independent. In the sections §3, §4 and §5, different model requiring different levels of sophistication are used to derive results that all agree qualitatively. Possible reading paths can be suggested as follows. In a first reading, one can focus on the double-beam model in §3, skipping the discussion of the expansion method in §3.6. The plate model in §4 offers an intermediate option: it presents all the main ideas of the papers in a mostly analytical setting. The plate model can also serve as a first illustration of the expansion method: the interested reader is advised to read §2 and then §4.4. The most complete results are presented in §5, where we make use of finite elasticity and of finite-element simulations.

Small-wavenumber expansion of the bifurcation modes

We consider a long prismatic solid having inhomogeneous pre-stress in its cross-section. The displacement of its endpoints is prescribed and the a mean imposed stretch is denoted by λ. In this section we propose a generic bifurcation analysis of the configuration that is invariant along the axis z, towards a buckled configuration characterized by an axial wavenumber q. The first critical wavenumber encountered when the imposed stretch λ is decreased is denoted by q c . As we consider an infinitely long solid, L = ∞, the bifurcation mode is macroscopic (long-wavelength) when q c = 0, and microscopic (short-wavelength) when q c > 0. The bifurcation occurs along the so-called 'curve of marginal stability' that is traced out in the λ, q) plane, see figure 2.1. This curve is obtained by solving the equations for adjacent equilibria near the cylindrically invariant solution.

Near the point λ (0) , 0 where it intersects the λ-axis, the curve of marginal stability can be fitted by a parabola as λ = λ (0) + q 2 λ (1) (dotted curve in the figure). In this section we propose an expansion method that allows one to calculate the coefficient λ (1) . The sign of λ (1) determines the direction in which the curve of marginal stability is curved near λ (0) , 0 : as illustrated in figure 2.1, a long-wavelength instability corresponds in general to λ (1) < 0, while a short-wavelength instability corresponds to λ (1) > 0. By allowing one to calculate λ (1) , the proposed expansion method can therefore predict whether the bifurcation mode is microscopic or macroscopic.

In this section, the expansion method is established in a general and formal setting. The method will be applied later to specific structural models, namely a double-beam ( §3.6), a plate ( §4.4), and a 3-d elasticity model ( §5.7). The forthcoming section is technical in some places. In a first reading, it is sufficient to read §2.1 (geometry), §2.3 (canonical form of the eigenvalue problem), §2.4 (form of the expansion) and §2.11 (summary of steps required to calculate λ (1) ).

Geometry, invariant solutions

We use a cylindrically invariant configuration of the solid as the reference configuration. Let (x, y, z) be a set of Lagrangian variables, such that the axis z is parallel to the axis of the solid in this reference configuration. We use (x, y, z) as Lagrangian variables to track the subsequent deformation of the solid, with (x, y, z) denoting the position of a material point in the reference configuration. In its reference configuration, the solid is pre-stressed. The pre-stress may be inhomogeneous in the cross-section (x, y) but is independent of z. More specifically, we assume that both the pre-stress distribution and the material properties of the solid are invariant both by a translation along the axis z (z → z + τ ) and by reflections about the plane of the cross-section (z → -z).

The solid is deformed by prescribing the displacement of its remote ends, which impose a mean stretch λ. Let ϕ(z) denote the collection of degrees of freedom that characterize the configuration of the cross-section with Lagrangian coordinate z. In a finite-element discretization, for instance, ϕ(z) is the vector collecting the degrees of freedom of the nodes that belong to the cross-section z.

For any value of the stretch λ there exists an equilibrium solution which is invariant both by translation along z and by reflection about a cross-section. This invariant solution is of the form ϕ(z) = ϕ λ 0 : the degrees of freedom depends on the mean imposed stretch λ and on the cross-sectional coordinate, but not on the axial coordinate z.

Initial form of the polynomial eigenvalue problem

We analyze the bifurcations that can take place when the mean imposed stretch λ is decreased: such bifurcations give rise to another family of solutions that are not invariant in the axial direction. The bifurcation analysis starts by introducing perturbations that depend harmonically on the axial coordinate z, as imposed by the invariance of the base solution: in complex notation,

ϕ(z) = ϕ λ 0 + ξ 1 e i q z + • • • (2.1)
Here, q is the wavenumber and ξ 1 is a vector collecting the complex amplitudes associated with each one of the degrees of freedom on a particular cross-section. As we shall show in the different illustrations, ξ 1 can be taken to be a real vector as this amounts to impose the phase of the buckling mode.

Next, the perturbation (2.1) is inserted into the equations of equilibrium, and the latter are linearized with respect to the amplitude ξ 1 . The form of the equations of equilibrium depend on which particular structural model is used, and they will be specified later. We work with the weak form of the equilibrium equations, as obtained by the principle of virtual work. For an infinitely long solid, the only virtual motions that produce non-zero virtual work with the particular family of perturbations (2.1) are those that depend harmonically on the axial variable z as well, and have the same wavenumber q: we consider virtual motions of the form ϕ(z) = ξ e i q z . In weak form, the linearized equilibrium equations take the generic form:

∀ ξ, ξ • d 0 λ • ξ 1 + q ξ • d 1 λ • ξ 1 + • • • + q n ξ • d n λ • ξ 1 = 0, (2.2)
as we will check later. This defines an eigenvalue problem having the wavenumber q as the eigenvalue and the mode shape ξ 1 as the eigenvector. As we have used Fourier analysis, the coordinate z has been replaced by the wavenumber q: the eigenvalue problem (2.2) is formulated on a the cross-section, which is a 2-d domain.

The order n and the bilinear operators d 0 λ , . . . , d n λ depend on the particular structural model chosen (beam, plate, or 3D elasticity): they will be calculated explicitly in the applications in §3-5. In a context of finite elasticity, these operators may depend on the mean imposed stretch λ through the base solution ϕ λ 0 . 5 They also depend on the distribution of inhomogeneous pre-stretch, although this is implicit in our notation. These operators satisfy symmetry properties, as a consequence of the invariance of the solid by reflection about cross-sections (z → -z): if, for instance, all the degrees of freedom are displacements in the plane of the cross-section (and are therefore unaffected by the reflection symmetry), the invariance q → -q dictates that the eigenvalue problem (2.2) contains only even powers of q. The detailed symmetry properties of the operators will be discussed later in §3-5.

For simple structural models such as beams ( §3), there is a finite number of degrees of freedom in each cross-section: then, ξ 1 is truly a vector, the operators d j λ can be represented as matrices and the virtual vector ξ can be readily eliminated from (2.2). For more complex structural models ( §4-5), however, a cross-section has infinitely many degrees of freedom, and ξ 1 is typically a collection of functions defined over the cross-sections. In the latter case, the dot products in (2.2) are meant to involve integrals over the cross-section, and the elimination of the virtual displacement ξ is not always desirable: this is why we retain the virtual displacement ξ in (2.2).

By repeatedly solving the eigenvalue problem (2.2) for different values of λ, one can plot the curves of marginal stability in the (λ, q) plane for a given distribution of pre-stress. This approach will be referred to as the direct method later on. The curves of marginal stability obtained in this way (see §3-5) are typically similar to those sketched in figure 2.1.

Solving polynomial eigenvalue problems is a numerically intensive task. In the rest of this section, we propose an expansion method that yields directly the 'initial curvature' λ (1) of the curve of marginal stability: this allows one to predict whether the first bifurcation mode is macroscopic or microscopic, without the need for an extensive eigenvalue analysis.

Canonical form of the eigenvalue problem

As mentioned above, the order n of the polynomial eigenvalue problems (2.2) depends on which structural model is used to represent the solid body. To unify the rest of the presentation, we start by observing that for, any value of n, equation (2.2) can be rewritten as a quadratic eigenvalue problem containing no linear term in q:

∀ Ξ, Ξ • A λ • Ξ + q 2 Ξ • B λ • Ξ = 0, (2.3a)
where A λ and B λ are symmetric operators, and the null space of A λ is independent of λ (these important properties are recapitulated in §2.5 below). In this equation, we have introduced a generalized eigenvector Ξ, by concatenation (direct sum) of the original eigenmode ξ 1 and of vectors of the form q k ξ 1 for some well-chosen powers k: in block-vector notation,

Ξ =   ξ 1 q k ξ 1 • • •   (2.3b)
In (2.3a), the virtual motion Ξ is an arbitrary vector having the same length as Ξ: it is not necessarily of the special form (2.3b).

To avoid a overly formal proof, we will admit the canonical form (2.3) of the eigenvalue problem, without attempting to derive it from its original form (2.2): this will be done later for each of the particular structural model, see §3-5. Doing so, we will obtain an explicit construction of the generalized eigenvector Ξ, as well as an explicit definition of the operators A λ and B λ appearing in (2.3a).

Expansion

At this point, we have outlined the derivation of the eigenvalue problem (2.3a) governing the linear bifurcation analysis. We can proceed to identify solutions in the vicinity of the bifurcation point λ, q = λ (0) , q = 0 depicted in figure 2.1.

To to so, we use q as an expansion parameter and seek a solution in the form

Ξ = Ξ (0) + q 2 Ξ (1) + q 4 Ξ (2) + • • • (2.4a) λ = λ (0) + q 2 λ (1) + q 4 λ (2) + • • • (2.4b)
The absence of odd powers of q in these expansions is imposed by the symmetry q → -q of (2.3a), which is itself a consequence of the symmetry of the homogeneous solution by a reflection about the plane (x, y).

We will need an expansion of the operators A λ and B λ as well. They are known in terms of the load λ; in view of (2.4b), their expansion can be calculated as

A λ = A (0) + q 2 λ (1) A (0) + • • • (2.5a) B λ = B (0) + q 2 λ (1) B (0) + • • • (2.5b)
where we use the notation

A (0) = A λ=λ (0) , B (0) = B λ=λ (0)
and

A (0) = dA λ dλ λ=λ (0) , A (0) = d 2 A λ dλ 2 λ=λ (0) , B (0) = dB λ dλ λ=λ (0) .
The notation A (0) will turn out to be useful later on.

Properties of the operators

As a preliminary step, one should check the properties of A λ and B λ announced below equation (2.3a):

A λ = A λ T , B λ = B λ T , ker A λ is independent of λ.
(2.6*)

In the absence of ambiguity, the null space will then be denoted as ker A, omitting the subscript λ.

For any null vector Ξ in ker A, the equality A λ • Ξ = 0 holds for any λ: by differentiating successively with respect to λ, this yields a useful property,

∀Ξ ∈ ker A, A (0) • Ξ = 0, A (0) • Ξ = 0, A (0) • Ξ = 0, • • • (2.7)
2.6. Solution at order 0 We are now ready to start the proper expansion procedure. Inserting the expansions (2.4) and (2.5) into the eigenvalue problem (2.3a) and reading off the result at order q 0 = 1, we find Ξ • A (0) • Ξ (0) = 0 for any Ξ. In different words, Ξ (0) is in the kernel of A λ , Ξ (0) ∈ ker A.

(2.8)

As we shall show later, ker A contains rigid-body modes of deformation, so we have just shown that the buckling modes living on the curve of marginal stability tend to rigid-body modes when the curve of marginal stability meets with the λ-axis, i.e. for q → 0.

Eigenproblem at order 2

Reading off the eigenvalue problem (2.3a) now to order q 2 , and noting that the term A (0) • Ξ (0) = 0 by equations (2.7) and (2.8), we obtain

∀ Ξ, Ξ • A (0) • Ξ (1) + Ξ • B (0) • Ξ (0) = 0.
(2.9*)

This equation is considered in the two following sections: a solvability condition is derived in §2.8, which provides an implicit equation for the stretch λ (0) , and the equation is solved in §2.9 for the linear bifurcation mode Ξ (1) .

Solvability condition at order 2

This linear equation for Ξ (1) involves the matrix A (0) which is singular by equation (2.8). Before attempting to solve it for Ξ (1) , one must therefore enforce a solvability condition (Fredholm alternative); this condition is derived by considering virtual motions Ξ belonging to ker A λ = ker A λ T which we denote by ker A (see 2.6): from (2.9), this yields

∀ Ξ ∈ ker A, Ξ • B (0) • Ξ (0) = 0.
(2.10) Let B * (λ) denote the restriction of the operator B λ to the subspace ker A. By equation (2.10), Ξ (0) is a null vector of B * (λ (0) ), implying that B * (λ) a singular operator when λ = λ (0) . In all the examples analyzed later, the dimension of the vector space ker A is small (at most 8 in our particular examples) even when the original eigenvalue problem is infinite-dimensional: it will be convenient to represent the operator B * (λ) as a square, symmetric matrix. The condition that B * (λ (0) ) is singular reads det B * (λ (0) ) = 0.

(2.11*)

This is an implicit equation for the load λ (0) at which the branch of marginal stability meets with the λ-axis in figure 2.1. Equation (2.10) can be rewritten as

Ξ (0) ∈ ker B * (λ (0) ), (2.12a*)
which is a stronger statement than (2.8) as ker B * (λ) ⊂ ker A.

Let n * ≥ 1 denote the dimension of ker B * (λ (0) ),

n * = dim ker B * (λ (0) ). (2.12b*)
As Ξ (0) lives in ker B * (λ (0) ) by (2.12a), it is defined by its n * coordinates in a basis of B * (λ (0) ). We will determine these coordinates later on in the expansion, up to a global scaling factor.

2.9. Solution at order 2 Having imposed the solvability condition for the equation (2.9), we can proceed to solve this equation. Let Ξ

(1) p denote a particular solution, which can typically be found by solving a linear algebra problem (discrete case, see §5.7), or by rewriting equation (2.9) in strong form and solving the resulting differential equations (continuous case, see §4.4)-here, by 'discrete' or 'continuous' cases, we mean that there are finitely many, or infinitely many degrees of freedom in the cross-section, respectively.

The general solution of (2.9) is then the sum of this particular solution and a solution of the homogeneous problem, Ξ (1) = Ξ (1) p + Ξ

(1) h , with Ξ

(1) h ∈ ker A.

Solvability condition at order 4

When expanded to order q 4 , equation (2.3a) yields a linear equation for Ξ (2) . As earlier, it involves the singular operator A (0) and a solvability condition must be enforced before one attemps to solve for Ξ (2) . The solvability condition is found again by taking the virtual motion Ξ in ker A: this cancels the term involving the unknown Ξ (2) , as well as other terms containing the dot products Ξ

• A (0) = 0 and Ξ • A (0) = 0, see (2.7). The result is ∀ Ξ ∈ ker A, Ξ • B (0) • Ξ (1) p + Ξ • B (0) • Ξ (1) h + λ (1) Ξ • B (0) • Ξ (0) = 0.
Consider the special case Ξ ∈ ker B * (λ (0) ) and recall ker B * (λ (0) ) ⊂ ker A: the second term in the left-hand side then cancels, as Ξ

• B (0) • Ξ (1) h = Ξ (1) h • B (0) • Ξ = Ξ (1) h • B * (λ (0) ) • Ξ = Ξ (1)
h • 0 = 0. The solvability condition takes the simple form

∀ Ξ ∈ ker B * (λ (0) ), Ξ • B (0) • Ξ (1) p + λ (1) Ξ • B (0) • Ξ (0) = 0.
(2.13*)

Equation (2.13) yields n * equations when Ξ is successively replaced by n * independent vectors in ker B * (λ (0) ). As we shall check later, these n * equations yield (i ) the value of λ (1) and (ii ) the value of the n * unknown coefficients in Ξ (0) , see equation (2.12a), up to a global scaling factor. In other words, the expansion method presented here allows one to calculate the 'initial curvature' λ (1) of the curve of marginal stability, as announced earlier, as well as the direction of the vector Ξ (0) . The magnitude of Ξ (0) remains undetermined, as the eigenvalue problem (2.3a) we started from is linear with respect to Ξ.

In the case n * = 1, for instance, ker B * (λ (0) ) is spanned by Ξ (0) : setting Ξ = Ξ (0) in equation (2.13), we find

λ (1) = - Ξ (0) • B (0) • Ξ (1) p Ξ (0) • B (0) • Ξ (0) (n * = 1). (2.14)
Note that this value of λ (1) is independent of both the arbitrary norm of Ξ (0) ∈ ker B * λ (0) -as Ξ

(1) p scales proportionally with this norm-, and of which particular solution Ξ

(1) p of (2.9) has been chosen. The case n * > 1 is treated similarly: an example is worked out in §5.7.

Summary

Among the equations listed above, only those marked by a star need to be considered for the purpose of applying the method: the other equations are required for the proof only. Practically, the expansion is carried out by the following sequence of steps:

1. given the operators A λ and B λ , check the properties (2.6*); 2. calculate the restriction B * of B λ to ker A and determine the critical load λ (0) by solving the implicit equation (2.11*); 3. determine the dimension n * of ker B * (λ (0) ) and introduce the n * unknown components of Ξ (0) in a basis of ker B * (λ (0) ), see (2.12a*-2.12b*); 4. find a particular solution Ξ

(1) p of the eigenvalue problem (2.9*) at order q 2 in terms of the n * components of Ξ (0) ; 5. using equation (2.13*), derive n * equations for the n * components of Ξ (0) and for λ (1) ; 6. solve these equations for the value of λ (1) and for the direction of the eigenmode Ξ (0) (its magnitude remains undetermined).

This general procedure is implemented in §3.6, §4.4 and §5.7: in each case we start by deriving the expressions of the operators A λ and B λ , which vary from one structural model to the other (double-beam, plate or 3-d prismatic solid).

2.12. Comment: similarities and differences with a Koiter expansion The expansion method proposed above has some similarities with the well-known Lyapunov-Schmidt-Koiter method (van der Heijden, 2008). Both methods deal with singular operators and therefore involve solvability conditions. Equation (2.11) is similar to the condition that the total stiffness operator is singular at the critical load; equation (2.12a) is similar to the condition that the linear mode is a null vector of the total stiffness operator; equation (2.12b) defines the multiplicity of the mode; equation (2.13) is similar to the amplitude equation, that sets the amplitude of the perturbation as a function of the increment of load.

There is an important difference, however: the Lyapunov-Schmidt-Koiter is a non-linear expansion of the amplitude as a function of the load, while in our method an expansion of the critical wavelength as a function of the load is sought based on linearized equations of equilibrium-those equation are still non-linear with respect to the wavenumber q.

A double-beam model

In this section, we analyze the bistrip experiment sketched in figure 1.1 based on a simple double-beam model, which captures the important features of the instability at a qualitative level, as we shall see. Using this double-beam model, we calculate the first critical load and the corresponding wavenumber that govern the bifurcation away from the homogeneous solution (direct method), and discuss in particular whether the bifurcation mode is macroscopic or microscopic. The double-beam model: δ is a mismatch strain, and an imposed mean strain (imposed by the displacement of the the remote endpoints). For the sake of legibility, an extensional pre-stretch is shown, + δ/2 > -δ/2 > 0; in reality the structure does not becomes unstable unless of the pre-stretch is compressive, -δ/2 < 0.

Formulation

In the double-beam model, each half of the bistrip used in the original experiment in figure 1.1 is represented by a different beam: the beams are denoted by I and II, consistent with the notation in figure 1.1: the beam II represents the strip that has been pre-stretched prior to gluing with I as sketched on figure 3.1, and a more tensile (less compressive) pre-stress is assigned to beam II than to beam I, see below. The two beams are coupled by an elastic layer of springs. In the original experiment, the instability is truly threedimensional. By contrast, the double-beam model proposes a simplified two-dimensional representation, whereby the centerlines of each half of the original bistrip have been projected onto a common plane (Oyz). The projection takes place along the direction x initially parallel to the width of the bistrip.

Let z denote the axial coordinate and y I (z) and y II (z) the deflection of each beam. According to the notation introduced in §2.1, the configuration of a cross-section of the double-beam is denoted by

ϕ(z) = y I (z), y II (z) .
The strain energy of the double-beam per unit length is

W (ϕ) = 1 L L 0 1 2   i=I,II C y i 2 + G + χ i δ 2 y i 2 + R (y 1 -y 2 ) 2   dz, (3.1)
where χ i = ±1 is a sign,

χ i = -1 for beam i = I, +1 for beam i = II.
In equation (3.1), C = E I is the bending modulus of the beams (where E is the Young's modulus and I the geometric moment of inertia of the cross-section), G = E A is their traction modulus (where A is the area of the cross-section of one beam), and R is the modulus of the elastic layer of springs connecting the two beams, see figure 3.1. For each beam, the first term C y i 2 /2 is the bending energy, and the second term is the energy associated with the pre-strain ± δ/2 corresponding to a pre-stress G ( ± δ/2)). The beams are subjected to different amount of pre-strain: the mean imposed pre-strain 1 is the result of the displacement 1 As the double-beam model is formulated in small-strain context, the mean imposed strain is a loading parameter; it replaces the mean imposed stretch λ used in the rest of the paper, = λ (p+1) 2

-1 ≈ λ -1 (when the pre-stretch is small: p ≈ 1).

applied at the remote endpoints, while the difference in pre-strain δ arises from the mismatch strain prior to gluing the two strips, see figure 3.1(b). The last term in equation 3.1 captures in a qualitative way the strain energy associated with the deformation of the cross-section in the original experimental set-up. For a general discussion of mechanical models obtained by assembling two elastic rods, including the extension to three dimensions and to general constitutive laws, see the work of [START_REF] Lessines | Morphoelastic rods. part II: Growing birods[END_REF].

As in Euler buckling, the pre-strain δ and the mean imposed strain can make the homogeneous solution ϕ(z) = (0, 0) unstable. The corresponding bifurcation problem is studied in the rest of this section. To ease the discussion, we first reformulate the model in terms of dimensionless quantities.

Dimensionless variables

In terms of the model parameters, one can define a typical mismatch strain δ † and a typical wavenumber

q † as δ † = √ C R G , q † = R C 1/4
, as well as the following dimensionless quantities: wavenumber q, axial coordinate z, mismatch strain δ, average strain , and deflection y i (z)

q = q q † , z = z q † , δ = δ δ † , = δ † , y i (z) = y i (z = z/q † ). (3.2)
With W = W/R denoting the rescaled strain energy per unit length and ϕ(z) = y 1 (z), y II (z) the configuration of the centerline, the energy (3.1) of the double-beam can be written in dimensionless form as

W (ϕ) = 1 L L 0 1 2 ϕ (z) • c • ϕ (z) + ϕ (z) • b • ϕ (z) + ϕ(z) • a • ϕ(z) dz, (3.3a)
where the dimensionless length is again infinite, L = q † L → ∞, and we have introduced the 2 × 2 symmetric matrices

a = 1 -1 -1 1 , b = -δ 2 0 0 + δ 2 , c = 1 0 0 1 . (3.3b)
Note that the operator b acting on the first gradient ϕ depends on both the mismatch strain δ and on the mean strain , but the dependence on δ is implicit in our notation.

Eigenvalue problem governing bifurcation

The homogeneous, unbuckled solution corresponds to ϕ 0 = (0, 0). We consider bifurcations from the homogeneous solution by seeking perturbed equilibria in the form y j (z) = 0 + y j 1 e i q z (j = I, II). To match up with the generic expansion ϕ(z) = ϕ 0 +ξ 1 e i q z introduced earlier in (2.1), we collect the complex amplitudes into a vector ξ 1 = y I 1 , y II 1 . It can be checked easily that the real part of y I 1 is uncoupled with the imaginary part of y II 1 and vice versa: as a result, one can assume without loss of generality that the complex amplitudes are real.

The linearized equations of equilibrium are obtained in weak form from (3.3a) as

∀ ξ, ξ • a • ξ 1 + q 2 ξ • b • ξ 1 + q 4 ξ • c • ξ 1 = 0. (3.4a)
This equation appears to be a special case of the general equation (2.2) derived earlier, when the operators are identified as d 0 λ = a, d 2 λ = b and d 4 λ = c, and d 1 λ = d 3 λ = 0 by symmetry (recall that the mean imposed stretch λ is replaced by the mean imposed strain in the current context of linear elasticity).

As there is a finite number of degrees of freedom, the strong form of the polynomial eigenvalue problem (3.4a) is found directly by eliminating ξ as

a + q 2 b + q 4 c • ξ 1 = 0. (3.4b)
This equation has both the eigenvector ξ 1 and the eigenvalue q as unknowns. * = √ 8 separates the regime of macroscopic buckling from the regime microscopic buckling. The dotted parabola is the result of the asymptotic analysis in equation (3.9) for δ = 5 and is checked to be consistent with the direct method (grey curve). The 'stable' and 'unstable' keywords indicate the expected region of stability and instability of the unbuckled solution. The critical values for the imposed average strain c and the associated wavenumber q c are denoted by black disks. (b) Wavenumber q c of the first unstable mode as functions of dimensionless mismatch strain δ, featuring the macroscopic to microscopic transition at the critical value of δ * . The dashed grey line represents the limit value of q c when δ → ∞; this corresponds to the classical model of a strut on a linear foundation. (c) The first unstable mode is an Euler-type of mode in the microscopic regime (y I 1 = y II 1 for δ < δ * ), and a strut-on-foundation-type of mode for large dimensionless mismatch strain

(|y II 1 | |y I 1 | when δ δ * ).

Direct solution: macroscopic versus microscopic buckling

In view of (3.4b), the bifurcation condition is det a + q 2 b + q 4 c = 0. Inserting the special form of the operators (3.3b), this yields

q 4 + q 2 + 1 2 - δ 2 q 4 4 + 1 = 0. (3.5)
The curves of marginal stability are obtained by solving this implicit equation in the plane ( , q) for different values of the dimensionless mismatch strain δ, see figure 3.2(a). The variable on the horizontal axis is the dimensionless average strain . The sequence of experimental steps depicted in figure 3.1 are as follows: the value of the dimensionless mismatch strain δ is fixed when the two rods are glued together; the bistrip is stretched (large positive , homogeneous solution is stable); the ends are then brought closer to one another (decreasing ) until a bifurcation is observed. This corresponds to the point c , q c denoted by a black disk in figure 3.2. In figure 3.2(a), all the curves of marginal stability meet with the horizontal axis at

(0) = 0, (3.6)
consistent with the fact that the buckling threshold of an Euler beam goes to zero as its length goes to infinity. Note that this (0) = 0 is different from the first critical strain c if the bifurcation mode is microscopic (see below).

As shown in figure 3.2, the double-beam model predicts a bifurcation from a macroscopic instability (q c = 0) at low values of the mismatch strain δ < δ * , to a microscopic instability (q c > 0) at large values of Equivalence with the Euler beam in the limit where the dimensionless mismatch strain is small, δ δ * : the curves of marginal stability in the plane ( , q) collapse onto the master curve corresponding to the Euler beam model, for which q 2 + = 0 (dashed curve). (b) Equivalence with the strut-on-a-foundation in the limit where the mismatch strain is large, δ → ∞: bifurcation curves in the plane ( -δ/2, q) collapse onto the master curve corresponding to the strut on a foundation, for which q 4 + ( -δ 2 ) q 2 + 1 = 0 (dashed curve).

The critical value of the dimensionless mismatch strain can be found by expanding (3.5) for small as q 4 q 4 + 2 (1 -δ 2 /8) ≈ 0, which leads to δ = √ 8.

(3.7)

Limiting cases

The double-beam model has two interesting limits.

In the absence of mismatch stress, δ = 0, or for stiff springs, R → ∞, a condition for the energy to be minimum is that y I (z) = y II (z): in the double-beam model, the two beams move together to avoid the energy penalty associated with the springs. The model is then equivalent to a single Euler beam with the pre-stress G . In dimensionless variables, this limits corresponds to δ = 0.

The other limit is when both → +∞ and δ → +∞ but with -δ/2 remaining finite: then beam II remains straight (y II (z) is constant) to avoid the penalty associated with the strongly stabilizing term G ( + δ/2) y II 2 /2; the classical model of a strut on a linear foundation is then recovered for the other beam I, with pre-stress G ( -δ/2). In dimensionless variables, this corresponds to both → ∞ and δ → ∞, with -δ/2 remaining finite.

The convergence of the solutions of the double-beam model is confirmed in figure 3.3, where the curves of marginal stability are shown to collapse onto a master curve corresponding to an Euler beam for small δ, and to a strut on a foundation for large δ.

The Euler beam model and the strut on an elastic foundation are the textbook examples of elastic continua featuring a macroscopic or a microscopic instability, respectively. Interestingly, the double beam model 'interpolates' between them using the dimensionless pre-strain parameter δ.

A first application of the small-wavenumber expansion

In this section, we recover the existence of a macroscopic-to-microscopic transition based on the general small-wavenumber expansion presented in §2 and summarized in §2.11. For a simple model such as the double beam, this approach turns out to be more complicated than the direct method carried out in §3.4. Still, this first example of application of the small-wavenumber expansion is instructive; it allows us to warm up before we move on to more accurate (and difficult) structural models-the real power of expansion method will then become fully apparent.

We need to first rewrite the quartic eigenvalue problem (3.4a) as a quadratic one to match with equation (2.3a). The trick is to consider a new eigenvector Ξ of dimension 2 + 2 = 4 obtained by the direct sum (concatenation) of the original eigenvector ξ 1 and of q 2 ξ 1 :

Ξ = ξ 1 q 2 ξ 1 =     y I 1 y II 1 q 2 y I 1 q 2 y II 1     .
(3.8a)

Using block-matrix notation, we define two 4 × 4 matrices by assembling the 2 × 2 operators a, b and c as follows,

A = a 0 0 -c , B = b c c 0 . (3.8b)
The quartic eigenvalue problem (3.4a) can then be rewritten as a quadratic one,

A • Ξ + q 2 B • Ξ = 0. (3.8c)
Indeed, one can check easily that to any set of eigenvector ξ 1 and eigenvalue q of the original problem (3.4a) corresponds an eigenvector Ξ and an eigenvalue q of (3.8c). The converse is also true. Let Ξ = Ξ 1 Ξ 2 and q denote an eigenpair of (3.8c): expanding the blocks in (3.8c) we find two equations, a

•Ξ 1 +q 2 (b •Ξ 1 +c•Ξ 2 ) = 0 and -c • Ξ 2 + q 2 c • Ξ 1 = 0, which indeed yield equation (3.4b) upon elimination of Ξ 2 = q 2 Ξ 1 (observe that c is invertible).
We can now proceed to carry out the steps listed in §2.11, replacing the imposed axial stretch λ by the dimensionless mean strain :

1. A and B are symmetric and ker A, which is the line spanned by Ξ = (1, 1, 0, 0) is indeed independent of (note that ker A corresponds to rigid-body translations, y I 1 = y II 1 , as announced earlier); 2. B * ( ) is the restriction of B to ker A: it is a 1 × 1 matrix, i.e. a scalar whose value is calculated as

B * ( ) = (1, 1, 0, 0) • B • (1, 1, 0, 0) = 2 .
The condition det B * = 0 yields the load (0) as (0) = 0 and we recover equation (3.6); 3. at = (0) , the kernel of B * ( (0) ) = 0 is of dimension 1, so n * = 1. A general element Ξ (0) of ker B * (0) = ker A is the rigid-body translation Ξ (0) = α (0) (1, 1, 0, 0) where the magnitude α (0) ∈ R of the mode will remain undetermined;

4. a particular solution of the problem (2.9) at order q 2 is found in terms of α (0) as Ξ

(1)

p = α (0) δ 4 , -δ 4 , 1, 1 ; 5. with B (0) the block matrix B (0) = b (0) 0 0 0 and b (0) the upper-left block b (0) = db d = ( 1 0 0 1 ), we
have from equation (2.14), (1) = δ 2 8 -1. With (0) = 0, the osculating parabola = (0) + (1) q 2 to the curve of marginal stability at the point of intersection with the -axis is therefore

= q 2 δ 2 8
-1 (osculating parabola).

(3.9)

These results are fully consistent with the direct method of §3.4. First, the dotted parabola in figure 3.2 fits the curve of marginal stability accurately for δ = 5 near its intersection with the -axis. Second, equation (3.9) predicts the correct critical dimensionless mismatch strain δ = √ 8 corresponding to the microscopic-to-macroscopic transition, as obtained in equation (3.7) by a direct method.

Discussion

The double-beam model does not aim at reproducing the experiments accurately. Nevertheless it successfully captures the transition from a macroscopic buckling mode for low mismatch strain δ, to a microscopic buckling mode for higher mismatch strain δ. The nature of the instability is governed by a dimensionless parameter δ = G δ √ C R , showing that there is a competition between the stiffness R of the springs, favoring long wavelength, and the mismatch strain δ, favoring short wavelength. The spring constant R is itself meant to represent the deformability of the cross-sections in the 3-d bistrip.

For non-deformable cross-sections, R → ∞ or δ → 0, the Euler-Bernoulli beam model is recovered. The latter always predict a macroscopic instability, as is well known: (1) = -1 < 0 for δ = 0. This confirms that microscopic buckling cannot be explained by the classical Euler-Bernoulli beam model: the deformability of the cross-sections is a key ingredient.

The critical mismatch strain δ * has been calculated by two methods. As a closed-form equation for the curves of marginal stability are available for this particular model, a direct method was possible. The general small-wavenumber expansion yields the same results.

A plate model with inhomogeneous pre-stress

In this section, the bistrip experiment sketched earlier in figure 1.1 is revisited based on a plate model. An elastic plate with width v, thickness h, infinite length L, Young's modulus E and Poisson's ratio ν is subject to an inhomogeneous pre-stress, see figure 4.1. The plate model is accurate when the bistrip is thin, in the limit h v. Its predictions will be quantitatively correct in this limit, which is an improvement over the qualitative double-beam model.

The uniaxial membrane pre-stress is inhomogeneous and aligned with the direction z:

N 0 zz ( ; x) = Eh + χ(x) δ 2 , (4.1a)
where χ(x) = ±1 is now the piecewise continuous function

χ(x) = -1 in half-plate I, 0 < x < v/2, +1 in half-plate II, v/2 < x < v. (4.1b)
The difference membrane pre-stress N 0,II zz -N 0,I zz = Eh δ comes from the mismatch strain δ, i.e. from the strain δ imposed to II prior to gluing with I. By contrast, the average membrane pre-stress (N 0,I zz +N 0,II zz )/2 = Eh arises from the mean strain , as imposed the prescribed distance between the remote ends of the plate. During the experiment, the mismatch strain δ remains fixed; the imposed mean strain is initially positive, and then decreases as the ends are brought closer to one another, until a bifurcation is observed.

Plate model

As a thin plate buckles with infinitesimal strain, the mean strain , the mismatch strain δ, and the in-plane displacements will all remain small. The buckling bifurcation is described using the Föppl-von Kármán plate model, linearized near the planar pre-stressed configuration. Near any planar configuration, the buckling mode is purely transverse and we denote by w(x, y) the deflection and by ŵ(x, y) a transverse virtual displacement. In weak form, the linearized equilibrium equation of the plate reads, see for instance [START_REF] Audoly | Elasticity and geometry: from hair curls to the nonlinear response of shells[END_REF],

∀ ŵ(x, y), 1 L L 0 v 0 ŵ,αβ (x, y) m αβ (x, y) + N 0 zz ( ; x) ŵ,z (x, y) w ,z (x, y) dx dz = 0, (4.2a)
where L → ∞ is the length of the plate, a comma in subscript denotes a partial derivative, Greek indices such as α and β are restricted to the in-plane directions (they can take on the values x or y) and an implicit summation is implied over repeated indices according to Einstein's summation convention.

In equation (4.2a), m αβ denotes the bending stress. For a homogeneous isotropic material, it is given by the constitutive relations

m αβ = D (1 -ν) w ,αβ + ν δ αβ w ,γγ . (4.2b)
Here, D = E h 3 12 (1-ν 2 ) denotes the plate's bending modulus and δ αβ is Kronecker's symbol. The next step is to consider perturbations that are harmonic in the axial direction z. At the same time, we move to dimensionless quantities and use the width v as the natural scale for in-plane coordinates. In terms of the rescaled axial wavenumber q = q v and of the rescaled axial coordinate z = z/v, we denote by ξ 1 (x) and ξ(x) the complex amplitudes corresponding to the real and virtual deflections, respectively, w(x, y) = ξ 1

x v e i q z (4.3a) ŵ(x, y) = ξ x v e i q z . (4.3b)

By the same argument as earlier, the complex amplitudes can be assumed to be real without any loss of generality, as this amounts to set the phase of the bifurcating mode.

When the harmonic deflections are inserted into (4.2a) and the integration over z is carried out, one obtains the equation for the critical mode as a quartic eigenvalue problem, ∀ ξ(x), a( ξ, ξ 1 ) + q 2 b ( ξ, ξ 1 ) + q 4 c( ξ, ξ 1 ) = 0.

(4.4a)

The plate model has infinitely many degrees of freedom: in the eigenvalue problem, the eigenvector is now the function ξ 1 (x). The bilinear operators are found as

a( ξ, ξ 1 ) = 1 0 ξ ξ 1 dx (4.4b) b ( ξ, ξ 1 ) = 1 0 2 ξ ξ 1 + + χ(x) 2 δ ξ ξ 1 -ν d 2 ( ξ ξ 1 ) dx 2 dx (4.4c) c( ξ, ξ 1 ) = 1 0 ξ ξ 1 dx, (4.4d)
where χ(x) = χ(v x) has, like the original χ(x), the value -1 in the half-plate labelled I (0 < x < 1/2) and +1 in the half-plate labelled II (1/2 < x < 1). The dimensionless mismatch strain δ and mean imposed strain have been identified as As with the double-beam, we work with infinitesimal strain, and the mean imposed axial stretch λ has been replaced by the mean imposed dimensionless strain . In addition, we are now working in an infinite-dimensional setting, with ξ 1 being a function and not a vector: the dot notation used earlier for the evaluation of operators, such as ξ •a•ξ 1 , has been replaced by the parentheses notation a( ξ, ξ 1 ), which involves an integral over the cross-section, see equations (4.4b-4.4d).

δ = δ δ † , = δ † , where δ † = 1 12 (1 -ν 2 ) h 2 v 2 . ( 4

Direct approach: numerical shooting method

To analyze the solutions of the eigenvalue problem (4.4a), we start by integrating by parts so as to eliminate the virtual deflection ξ(x). This yields an equilibrium equation in the interior,

-ξ 1 (x) + 2 q 2 ξ 1 (x) -q 4 + q 2 + χ(x) δ 2 ξ 1 (x) = 0 (0 ≤ x ≤ 1), (4.5a)
stress-free boundary conditions,

ξ 1 -ν q 2 ξ 1 x i 0 = 0 (4.5b) ξ 1 -q 2 2 -ν ξ 1 x i 0 = 0, (4.5c)
as well as kinematic and dynamical conditions at the interface,

ξ 1 , ξ 1 , ξ 1 , ξ 1 are continuous at x = 1/2. (4.5d)
In (4.5b-4.5c), x I 0 = 0 and x II 0 = 1 denote the dimensionless coordinates of the stress-free edges. An alternative derivation of the linearized equilibrium (4.5a) is to combine the harmonic Ansatz (4.3a) for the deflection with the linearized von Kármán equation,

-D ∆ 2 w(x, z) + ∂ N 0 zz ( ; x) w ,z (x, z) ∂z = 0,
and by using dimensionless variables: this yields the same result. Equations (4.5b) and (4.5c) can also be obtained similarly from the boundary conditions at the stress-free edge of an elastic plate, w ,xx + ν w ,zz = 0 and w ,xxx + 2 -ν w ,xzz = 0, as derived for instance by [START_REF] Landau | Theory of Elasticity (Course of Theoretical Physics)[END_REF].

To solve the eigenvalue problem (4.5a-4.5d), we used a numerical shooting method. For any value of the parameters (Poisson's ratio ν, mismatch strain δ, mean imposed strain , wavenumber q), a 4 × 4 matrix S(ν, δ, , q), called the shooting matrix, is filled in, based on the numerical integration of the differential equation (4.5a) with particular initial conditions. For fixed values of ν and δ, the curves of marginal stability in the plane ( , q) are then found by solving numerically the implicit equation det S(ν, δ, , q) = 0.

(4.6)

Further details on the construction of the shooting matrix S are given in Appendix A.

Results

In the following, we consider an incompressible material, ν = .5. Curves of marginal stability corresponding to various values of δ are plotted in figure 4.2(a) based on the shooting method in (4.6). The phase diagrams are similar to those obtained for the double-beam model, and must be read in the same way. The plate is initially stretched ( > 0) and the imposed mean dimensionless strain is progressively decreased until a bifurcation to a non-planar solution is observed. The first bifurcation takes place at ( c , q c ), as denoted by the black disks in the figure. Again, we find that the first bifurcation mode is macroscopic when the mismatch strain is less than a critical value, δ < δ * (ν), and microscopic when it is greater, δ > δ * (ν).

The numerical results suggest that the critical value is δ * (ν = .5) ≈ 90. A closed-form expression for this δ * will be obtained later based on the small-wavenumber expansion, see §4.4. In figure 4.2(b), the critical wavenumber is plotted as a function of the mismatch strain. In figure 4.2(c), a measure of the curvature of the cross-section, ξ 1 (0)/ξ 1 (0) is plotted as a function of δ: the cross-section undergoes a rigid-body displacement when the mode is macroscopic (ξ 1 is a linear function of x, corresponding to an infinitesimal translation and an infinitesimal rotation, as in configuration labelled 'A' in the figure), and it bends in its own plane when the mode is macroscopic (ξ 1 (x) = 0, configurations 'B' and 'C').

The limits of small and large dimensionless mismatch strain δ are shown in figure 4.3. For small mismatch strain, δ δ * , the pre-stress is homogeneous and comes exclusively from the mean imposed strain : the limit curve has been found by setting2 directly δ = 0 in (4.6) (dashed curve in figure 4.3a). In the opposite limit δ δ * , the plate II is under large tensile stress N 0,II zz and the buckling mode is limited to region I: this is equivalent to a homogeneous half-plate I subject to pre-stress N 0,I zz = -δ/2 having clamped boundary conditions at the interface x = 1/2 (dashed curve in figure 4.3b); the corresponding asymptotic value of the wavenumber is q c = 3.9 for ν = .5 (δ δ * ).

Small-wavenumber expansion

With the aim to calculate analytically the critical mismatch strain at the macroscopic to microscopic transition, we apply the small-wavenumber expansion of §2.11.

The quartic eigenvalue problem (4.4a) has exactly the same form as that of the double-beam model, compare with (3.4a): it can be cast into a quadratic eigenvalue problem exactly in the same way. Therefore, we define two new operators as earlier,

A = a 0 0 -c , B = b c c 0 . (4.7a)
Unlike the original operators a, b , c which act on two functions (a virtual one ξ and an increment ξ 1 ), the operators A and B act on two pairs of functions, the virtual pair Ξ and the pair of increments Ξ. The pair of increments, for instance, reads as earlier in (3.8a):

Ξ = ξ 1 q 2 ξ 1 . (4.7b)
Here is an example of how the block notation in (4.7a) must be understood: for two virtual functions ξ and ξ † ,

ξ ξ † • B • ξ 1 q 2 ξ 1 = b ξ, ξ 1 + c ξ, q 2 ξ 1 + c ξ † , ξ 1 .
With these notations, the eigenvalue problem (4.4a) governing the bifurcation of the plate is rewritten as a quadratic eigenvalue problem that matches up with the general form (2.3a):

∀ Ξ, Ξ • A • Ξ + q 2 Ξ • B • Ξ = 0.
This is a two-way equivalence: by the same argument as earlier, any eigenvector Ξ of the equation above can be shown to be of the special form (4.7b) and to correspond to a solution of the original problem (4.4a).

We are now ready to calculate the critical value of the mismatch strain by the recipe given in §2.11:

1. A and B are symmetric. The kernel ker A has dimension 2 and it is spanned by pairs of real function ϕ 1 , 0 such that the first function is an affine function of x and the second function is identically 0. Like A, ker A is independent of ; 2. the restriction B * of B to ker A is calculated as

B * = 1 8 δ + 4 1 8 δ + 4 1 24 3 δ + 8(6 + -6ν)
, where we use the basis vectors 1 0 and x 0 in ker A (remember that each vector in the basis is a pair of functions). The first entry of B * , for instance, is calculated as

(B * ) 11 = B 1 0 , 1 0 = b (1), (1) = 1 0 + χ(x) 2 δ dx = .
Next, (0) (ν, δ) is found by solving det B * = 0 which yields the quadratic equation,

(0) 2 + 24 (1 -ν) (0) - 3 16 δ 2 = 0. (4.8)
To avoid cumbersome expressions, we will not attempt to solve this equation for (0) in closed-form; instead, we view it as an implicit equation for (0) . The first critical load is given by the positive root (0) > 0. The other root is negative and therefore corresponds to a secondary critical load; 3. the null space of B * at (0) is of dimension n * = 1, and is made up of pairs of functions of the form

Ξ (0) = α (0) ξ (0) 1 0 , where ξ (0) 1 (x) = 1 + β x and β = - 8 (0) δ + 4 (0) .
(4.9)

At this point, we have found the bifurcating mode, up to a factor α (0) which will remain undetermined;

4. a particular solution Ξ

(1) p of the eigenproblem (2.9) at order q 2 is found as Ξ

(1)

p = α (0) f (x) 0
where f satisfies the linear initial-value problem,

-f (x) = (0) + χ(x)) δ 2 1 + β x (0 ≤ x ≤ 1) f (0) = 0 f (0) = 0 f (0) = ν f (0) = (2 -ν) β.
Note that the initial conditions f (0) = f (0) = 0 have been chosen arbitrarily to specify a particular solution. The solution is:

f (x) = 1 1920 1 -χ(x) 2 (1 -2x) 4 β x + 2 β + 5 δ -8β δ + 2 (0) x 5 -40 δ + 2 (0) x 4 • • • + 40 β(δ -8ν + 16) + 4δ x 3 -40 (β + 3)δ -24ν x 2 + 5 δ (3β + 8) x + (-2β -5)δ ;
5. with n * = 1, the coefficient (1) yielding the curvature of the bifurcated branch in the phase diagram is given directly by (2.14) as

(1) = g ν ( (0) , δ), (4.10a)

where the function g ν is available in closed form: for ν = .5,

g ν=.5 ( , δ) = -49 δ 4 + 8 δ 2 71 2 + 700 -2520 -512 2 3 2 + 56 + 490 8960 3 δ 2 + 16 2 . (4.10b)
To evaluate the denominator in (2.14), we have used

B (0) Ξ (0) , Ξ (0) = b (0) α (0) ξ (0) 1 , α (0) ξ (0) 1 = α (0) 2 1 0 1 + β x
2 dx, and have used the value of β found in (4.9).

In figure 4.4, we have plotted both (0) (ν = 0.5, δ) as given by the implicit equation (4.8) (green curve), and the sign of (1) = g ν ( (0) , δ) from equation (4.10b) ( (1) = g ν ( (0) , δ) > 0 in the light red domain). As explained in the introductory paragraph of §2, the small-wavenumber expansion predicts that the first unstable mode is macroscopic whenever the green curve lies outside the light red domain ( (1) = g ν ( (0) , δ) < 0) and microscopic whenever it lies inside the light red domain ( (1) = g ν ( (0) , δ) > 0). For a graphical interpretation of (1) , see the dotted curve in figure 4.2(a). The corresponding ranges of the mismatch 4.2(a) (stars, disks and thin black curve) and predictions of the small-wavenumber expansion (thick green curve and red domain). The bifurcation mode is macroscopic when the green and black curves coincide, which is also when the green curve is outside the red region; it is microscopic otherwise. The prediction (0) EB of the Euler-Bernoulli rod model (dashed blue curve) is asymptotically consistent with that of the plate model in the limit of weak pre-stress, δ → 0. (b,b') Sketch of the curves of marginal stability in the q, plane corresponding to the macroscopic and microscopic cases, respectively. The expected domain of stability of the planar solution, as indicated by the 'stable' and 'unstable' labels, could be confirmed by a proper stability analysis. The critical dimensionless mismatch strain δ * (ν) at which the macroscopic to microscopic transition occurs can be found analytically by solving the two equations (4.8) and g ν ( (0) , δ) = 0 for the two unknowns δ and (0) . This yields a polynomial equation (not included here, for the sake of brevity) whose relevant root can be calculated to any desired accuracy. For ν = .5, we find: δ * (ν = .5) = 88.247 (4.12)

All these results are in good agreement with those obtained earlier by the direct method. The bifurcation is indeed macroscopic for low values of the mismatch strain δ, and microscopic for high values of δ. Moreover, the estimate δ * ≈ 90 found by the shooting method is close to the exact value (4.12) obtained by the smallwavenumber expansion. Finally, in figure 2.1(a) the curve passing through the black disks ( c ) indeed splits off from the curve passing through the green stars ( (0) ) at the critical value δ * predicted by (4.12).

4.5. Equivalence with a naturally curved rod model in the limit of a small pre-strain As discussed in the introduction, the naturally curved Euler-Bernoulli rod model has been shown to be applicable to the bistrip in the limit of a weak pre-stress using rigorous arguments [START_REF] Cicalese | On local and global minimizers of prestrained thin elastic rods[END_REF]. In line with this result, we show here that the predictions of the plate model are indeed asymptotically consistent with those of an Euler-Bernoulli rod model in the limit of small pre-stress.

A linear bifurcation analysis of the naturally curved Euler-Bernoulli rod yields the critical value of the mean strain as

EB c = (M 0 y ) 2 µJ E h v , (4.13)
as derived by [START_REF] Goriely | Spontaneous helix hand reversal and tendril perversion in climbing plants[END_REF]; [START_REF] Liu | Structural transition from helices to hemihelices[END_REF]. At this critical value EB c , long-wavelength helical solutions appear by a bifurcation from the straight solution. In equation above, M 0 y is the residual bending moment arising from the pre-stress,

M 0 y = v 0 N 0 zz ( ; x) x - L 2 dx = E h v 2 δ 8 , (4.14)
and (µJ) is the effective twisting modulus of the plate. The latter be found by identifying the bending energy per unit length associated with a twisting deformation, w(x, z

) = τ (x -L/2) z, with 1 2 µJ τ 2 : the result is µJ = 2 D (1 -ν) v = E h 3 v 6(1+ν)
. Combining (4.13), (4.14) and the above expression of µJ, we obtain the critical mean imposed strain predicted by the Euler-Bernoulli model as

EB c = (0) EB = (1+ν) 2 v h 2 3 δ 2
16 . In terms of the dimensionless parameters introduced in (4.4e), this can be rewritten as

24 (1 -ν) (0) EB = 3 16 δ 2 . (4.15)
This prediction is asymptotically equivalent with that derived earlier in (4.8) from the plate model: the only difference is that the term (0) 2 is absent from the left-hand side of (4.15); this term is indeed negligible in the limit of a weak pre-strain, as both δ and (0) are then small. The Euler-Bernoulli prediction (4.15) for macroscopic buckling is shown in figure 4.4 using a blue dashed curve: the agreement with the plate model is good for asymptotically small δ, but becomes poor as soon as δ reaches ∼ 10, i.e. for a pre-strain that is still much below the critical pre-strain δ * ≈ 90 at which microscopic modes appear. This shows that the Euler-Bernoulli model is of limited interest for pre-strained structures, even for the analysis of macroscopic modes.

Discussion

The plate model displays the same type of behavior as the double-beam model, with the advantage that it is expected to be quantitatively correct for slender bistrips, h v. We have identified a dimensionless mismatch strain δ that governs the nature of the first bifurcation mode: it is macroscopic for small mismatch strain δ and microscopic for large mismatch strain. The parameter δ measures the antagonistic effects of the mismatch strain (shorter wavelengths are preferred for larger mismatch strain) and of the bending rigidity of the cross-sections (longer wavelengths are preferred for more rigid cross-sections). The buckling behavior of the plate has been analyzed by a numerical shooting method, and the critical mismatch strain δ * corresponding to the macroscopic-to-microscopic transition has been obtained analytically using the small-wavenumber expansion, see (4.12). As a final note concerning the plate model, observe that the strain has been rescaled by the quantity δ † which goes to zero in the limit of a small aspect-ratio h v by equation (4.4e). As a result, the finite values of and δ calculated above correspond to infinitesimal values of the physical mean strain and strain mismatch δ: this confirms that it is sufficient to consider a small-strain plate model with linear constitutive assumptions.

A 3-d hyper-elastic model with inhomogeneous pre-stress

In this section, the experimental bistrip of [START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF]; [START_REF] Liu | Structural transition from helices to hemihelices[END_REF] is approached based an hyper-elastic, infinitely long, rectangular cuboid made of a Gent material. By using Fourier analysis in the axial direction, the bifurcation problem is formulated as an eigenvalue problem in the cross-section.

(1)

(2)

(3) (4) The bifurcations from the cylindrically invariant solutions are first analyzed by the direct method, which combines the finite-element method with numerical eigenvalue analysis: qualitatively, this yields similar results as with the simpler double-beam and plate models. This direct method is however computationally intensive. The small-wavenumber expansion is presented next, and provides an efficient and accurate method for characterizing the macroscopic-to-microscopic transition.

Model

Different configurations of the bistrips are shown in figure 5.1:

• In their natural configuration (1), the strips, labelled I and II, are rectangle cuboids with different cross-section dimensions h × (v/2) and h × (v /2), respectively.

• In configuration (2), the strip II is pre-stretched with axial stretch p and the strips are glued together.

The original dimensions h and v /2 have been chosen such that both strips have the same dimensions h × (v/2) in this configuration (2). The coordinates (x, ỹ, z) in this pre-stretched configuration are used as Lagrangian coordinates, and we denote by D the cross-section D = (0, v) × (0, h): (x, ỹ) ∈ D.

• The remote ends are moved towards one another and as a result an additional stretch λ is imposed (with λ < 1 as the bistrip is released). This yields the cylindrically invariant configuration (3), which may be stable or unstable: with respect to their natural configuration (1), the stretch is λ in the first strip I and p λ in the second strip II.

• A buckled, non-cylindrically symmetric, configuration (4) can appear through a bifurcation.

Different sets of coordinates are associated with each of these configurations: X = (X, Y, Z) for the natural configuration (1), x = (x, ỹ, z) for the pre-stretched configuration (2), x = (x, y, z) for the equilibrium (unbuckled or buckled) configurations (3-4). Consistent with standard notations for 3-d elasticity, we denote vectors by a simple underline and matrices by a double underline in this section. We further denote by e X , e Y and e Z the vectors of the Cartesian basis which remain unchanged in natural and pre-stretched configurations (1) and (2).

By contrast with the previous double-beam and plate models, we now work with finite strain. The pre-stretch p plays the role formerly assigned to the mismatch strain δ, and the stretch λ imposed by the displacements of the remote ends plays the role formerly assigned to the average strain . In the following, we refer to λ as the 'mean imposed stretch', using configuration (2) as a reference.

Let G denote the transformation gradient from configuration (1) to configuration (2):

G = 1 in region I, 0 < x < v/2, r(p) e X ⊗ e X + e Y ⊗ e Y + p e Z ⊗ e Z in region II, v/2 < x < v.
(5.1)

Here, r(p) denotes the transverse stretch associated with the axial stretch p under simple traction for the particular material law considered: denoting by W 3D (r, p) the strain energy per unit reference volume for an equi-biaxial stretch, r(p) is defined by the implicit equation ∂W 3D ∂r r(p), p = 0. As we use a nearly incompressible Gent model (see below), det G ≈ 1 and r(p) is close to p -1/2 . We denote by ϕ(x) the displacement from the pre-stretched configuration (2) to the final one (4): the final position reads3 x = x + ϕ(x) and the transformation gradient is F = 1 + ∇ϕ(x). With reference to the natural configuration (1), the total transformation gradient reads:

H(x) = ∂x ∂X = F • G. (5.2a)
This multiplicative decomposition is classical in the elastic theory of growth.

A Gent hyper-elastic material model is used and the strain energy per unit reference volume has the form

W 3D (H) = µ 2 -J m ln 1 - I c -3 J m -µ ln J + K 2 - µ J m J -1 2 (5.2b)
where K, µ, J m are material constants and I c and J are the invariants of the total transformation gradient,

I c = H : H, J = det H. (5.2c)
We use the standard semi-column notation for the doubly contracted product, A : B = tr(A T • B . In all the forthcoming numerical simulations, we used the following set of material constants: shear modulus µ = 1, bulk modulus K = 10 and J m = 100. The relatively large ratio K/µ = 10 implies that the material is weakly compressible: the initial Poisson's ratio can be calculated as ν = 0.45. In addition, the large value J m = 100 makes the constitutive law practically equivalent to a neo-Hookean model, up to values of the principal stretches as large as ∼ √ J m = 10, or as small as ∼ 1/ √ J m = 0.1. These values of the material parameters are not meant to match exactly those relevant to the experiments. All the methods, models and conclusions presented in this paper are to some extent independent of the choice of material parameters. In particular, all our results remain virtually unchanged if we use an incompressible neo-Hookean instead.

General equations for equilibrium

The average strain energy per unit length is given by the integral of the strain energy W 3D dX dY dZ =

W 3D det G dx dỹ dz as 1 L L 0 D W 3D (H(x)) det G dx dỹ dz, (5.3)
where we consider the limit of an infinitely long solid, L → ∞. The energy has to be minimized under the constraint that the displacement is prescribed at the remote ends, i.e. that the mean stretch matches the imposed one λ.

For any particular solution ϕ(x), consider the Green-Lagrange strain, e = 1 2 (H T • H -1), and the second Piola-Kirchhoff stress Σ = 1 det G ∂W 3D ∂e , which are both symmetric. For any virtual motion ϕ(x), define the virtual increment e of the Green-Lagrange strain as the symmetric part of H T • H, where H = F • G = ∇ ϕ(x) • G. The non-linear equilibrium solutions ϕ(x) are found by cancelling the first variation of the total energy for any admissible virtual displacement:

∀ ϕ(x, ỹ, z), 1 L L 0 D Σ(ϕ) : H T (ϕ) • H dx dỹ dz = 0.
(5.4)

As the endpoint positions are fixed, we consider admissible virtual displacements that have a zero average axial strain.

Fundamental branch: cylindrically invariant solutions

We will first compute the fundamental branch of solutions: for any value of the mismatch stretch p and of the stretch λ imposed by the loading, we seek a cylindrically invariant solution of the nonlinear equilibrium equation (5.4) in the form ϕ λ 0 (x) = ϕ 0 (λ; x, ỹ) + (λ -1) z e Z .

(5.5)

The associated transformation gradients read F 0 = 1 + ∇ϕ λ 0 (x) and H 0 = F 0 • G, see figure 5.1. The second Piola-Kirchhoff stress tensor has a block-diagonal form imposed by the symmetry of the homogeneous solution Σ 0 = Σ(ϕ λ 0 ) with:

Σ(ϕ λ 0 ) =   Σ 0 (λ; x, ỹ) 0 0 Σ ⊥ 0 (λ; x, ỹ)   (5.6)
where the upper-left 2 × 2 block Σ 0 denotes the cross-sectional stress in the (x, ỹ) plane, and the lower-right scalar Σ ⊥ 0 denotes the stress along the z-axis. Both transformation gradients F 0 and H 0 have a similar block decomposition.

When this particular form of solution is inserted into the general equilibrium equilibrium (5.4), one obtains a 2-d problem of non-linear elasticity whose main unknown is the cross-sectional displacement ϕ 0 (λ; x, ỹ). The axial direction z is absent from this problem, except for the imposed axial stretch λ which enters as a parameter. A numerical approximation of the solution ϕ 0 (λ; x, ỹ) can be obtained by the nonlinear finite-element method, for any set of values of the load parameters p and λ, see §5.5 below for the details of the implementation. Note that the dependence of the solution ϕ 0 (λ; x, ỹ) on the mismatch stretch p is implicit in our notation. We work in units such that the width is v = 1. The domain D = (0, 1) × (0, h) is rectangular, see configuration (2) in figure 5.1. Typical pre-stress distributions Σ 0 (cross-sectional) and Σ ⊥ 0 (axial) are plotted in figure 5.2. The axial stress is nearly constant in each domain I and II and the solid is in a stress state close to simple traction. In the present context of finite elasticity, Poisson's ratio depends on the pre-stretch and is therefore slightly different in the domains I and II. As a result, the simple-traction solution is (slightly) geometrically incompatible and a small amount of cross-sectional stress appears near the interface, see figure 5.2(b).

Bifurcation analysis: direct approach using the finite-element method

The invariant solution ϕ λ 0 is available numerically, and we proceed to the bifurcation analysis. We start with the direct method: Fourier analysis is used in the axial direction z, and the existence of an adjacent equilibrium is formulated as a quartic eigenvalue problem and then solved numerically using the finite element method.

Let ϕ 1 (x) denote a small increment of displacement adding up to an invariant solution: the displacement is expanded as ϕ 5.6). Aspect ratios h/v = 0.1, 0.2, 0.4 (from top to bottom); material parameters: µ = 1, K = 10, Jm = 100; mismatch stretch p = 1.11; mean imposed stretch λ = .9. The axial stress (left) is nearly uniform in each half of the cross-section: it is compressive in region I and tensile in region II (because of the compressive axial stress, this cylindrically invariant configuration is actually unstable, see figure 5.3a). The cross-sectional stress (right) arises from geometric incompatibility between the two domains, and is limited to a neighborhood of their interface; it is much smaller than the axial stress (note the widely different scales in the color bars).

(x) = ϕ λ 0 (x)+ϕ 1 (x)+• • • When linearized,
Here, H 1 = ∇ϕ 1 (x) • G is the incremental transformation gradient and L 0 is the tensor of tangent moduli evaluated in the invariant solution,

(L 0 ) ijkl (λ) = H 0 is ∂Σ sj ∂e rl H 0 kr ϕ λ 0 , (5.8)
where H 0 is denote the components of H 0 in the Cartesian frame. The complete expressions of L and Σ are given in Appendix B.1 for the Gent model.

In view of the symmetry of the base solution ϕ λ 0 , we consider bifurcation modes ϕ 1 that are harmonic in the axial variable z with a wavenumber q; such modes are coupled to virtual displacements that are harmonic and have the same wavenumber q:

ϕ 1 (x, ỹ, z) = ξ (x, ỹ) + i ξ ⊥ (x, ỹ) e Z e i q z , ϕ(x, ỹ, z) = ξ (x, ỹ) + i ξ⊥ (x, ỹ) e Z e i q z .
(5.9)

Here, ξ (x, ỹ) is a vector-valued function defined on the reference cross-section D which yields the complex amplitude of the cross-sectional displacement; ξ ⊥ (x, ỹ) is a complex function defined on D which yields the complex amplitude of the axial displacement. Their virtual counterparts are ξ (x, ỹ) and ξ⊥ (x, ỹ). In equation above, the imaginary unit number i that multiplies the axial displacements introduces a convenient phase shift, such that the real parts are entirely decoupled from the imaginary parts. As a result, we can assume all functions to be real, without loss of generality. The incremental and virtual displacements are discretized by the finite-element method on the 2-d domain D = (0, 1) × (0, h). By convention, degrees of freedom corresponding to displacement in the cross-section (x, ỹ) are listed first, while degrees of freedom corresponding to axial displacement appear last: in the discrete setting, the complex amplitudes are represented by vectors (using block notation)

ξ 1 =       ξ ξ ⊥       , ξ =       ξ ξ⊥       . (5.10a)
Here, ξ and ξ ⊥ represent the degrees of freedom of the functions ξ (x, ỹ) and ξ ⊥ (x, ỹ), respectively, and a similar convention is used for the virtual displacements. The additional vertical spacing around the columnvectors representing cross-sectional displacement ( ) emphasizes that they contain twice as many degrees of freedom as those representing axial displacement ( ⊥ ).

In equation (5.7) for the adjacent equilibria, the Fourier analysis effectively removes the integration over z, and transforms any gradient with respect to z into a multiplication by i q. In addition, the discretization replaces the cross-sectional integration by dot products. This transforms the equation into a polynomial eigenvalue problem of the form

∀ ξ, ξ • K λ + q C λ + q 2 M λ • ξ 1 = 0, (5.10b)
where the discretized operators K λ , C λ and M λ can be found by identification with (5.7). A complete definition of these operators is given in Appendix B.2. As shown there, the symmetry of the base solution warrants that the operators K λ and M λ associated with even powers of q are block-diagonal while operator C λ associated with the odd power q is block 'anti-diagonal',

K λ =       k λ 0 0 k ⊥ λ       , C λ =       0 c T λ c λ 0       , M λ =       m λ 0 0 m ⊥ λ       .
(5.11)

In addition, the square blocks k λ , k ⊥ λ , m λ and m ⊥ λ are all symmetric. The larger square blocks k λ and m λ act on cross-sectional displacements, in accord with the ordering conventions used in (5.10a); the smaller square blocks k ⊥ λ and m ⊥ λ act on axial displacement; the rectangular block c λ couples cross-sectional and axial displacements. All these matrices can readily be calculated in terms of the base solution ϕ 0 (λ, x, ỹ), as explained in the appendix; they encode information about both the pre-stress Σ 0 (λ, x, ỹ) that drives the bifurcation, and the tangent moduli L 0 (λ).

Note that the eigenvalue problem (5.10b) matches the generic form (2.2) announced earlier.

Numerical implementation

For given loading parameters (p, λ), the invariant solution ϕ 0 (λ; x, ỹ) is first computed by solving the 2-d non-linear elasticity problem described in §5.3. To this end, the finite-element method is used and the displacement ϕ 0 is interpolated using linear Lagrange elements. Convergence is obtained by a standard Newton-Raphson method. Our implementation makes use of the finite-element library FEniCS (Ølgaard et al., 2009).

In a second step, the incremental displacement (ξ , ξ ⊥ ) is discretized: we use again a finite-element discretization of the 2-d domain D using linear Lagrange elements and the FEniCS library. The corresponding basis of functions are used, together with the known invariant solution ϕ 0 , to fill in the matrices K λ , C λ and M λ , as explained in Appendix B.2.

The discrete form of the quadratic eigenvalue problem (5.10b) is then solved using the SLEPc library [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. A two-level orthogonal Arnoldi method is used, see [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF] for a review on numerical solution of quadratic eigenvalue problems; it is combined with shift-andinvert preconditioning, allowing one to compute eigenpairs corresponding to small real eigenvalues q more the critical wavenumber predicted by the 3D model using the direct method (dots connected by thin curves) converges to that predicted by the plate model (thick curve). For all aspect-ratios, the small-wavenumber expansion (star symbols on the horizontal axis) correctly predicts the critical mismatch strain δ * at which the bifurcation mode switches from macroscopic to microscopic. (c) First critical mode corresponding to p = 1.11 with h/v = 0.1, shown with an arbitrary amplitude.

efficiently [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF]. To remove eigenpairs that correspond to rigid-body modes, the symmetric part of H T 1 • H 0 , which is a measure of the linearized strain, is calculated for each solution (q, ξ 1 ) of (5.10b): whenever the norm of the linearized strain is below a small threshold, the solution is identified as a rigid-body mode and eliminated.

Results

The numerical results are shown in figure 5.3(a). The bifurcation diagram must be read like those in figures 2.1 (generic model), 3.2(a) (double-beam model) and 4.2 (plate model). In the experiment, the value of the mismatch stretch p is first fixed upon gluing the two strips. Next, the tension is released: the mean stretch λ imposed by the position of the remote ends decreases progressively from its initial value λ = 1, with a given step-size4 δλ. A bifurcation takes place at the point of the curve of marginal stability corresponding to the largest abscissa λ: this yields a macroscopic instability if the corresponding wavenumber q c is zero, and a microscopic instability otherwise.

As seen in the figure, the 3-d hyper-elastic model features both macroscopic and microscopic buckling, depending on the amount of pre-stretch p. As with the double-beam and plate models, a smaller pre-stretch favors macroscopic buckling (consistent with the Euler-Bernoulli beam model) while a larger pre-stretch favors microscopic buckling (consistent with a strut-on-elastic-foundation model). With the particular aspect-ratio h/v = 0.1 and material parameters chosen in this example (see legend), the transition occurs at the critical value of the mismatch stretch p * ≈ 1.1.

(5.12)

The results of the 3-d hyper-elastic model are therefore qualitatively similar to those of the simpler doublebeam and plate models ( §3-4). The critical modes predicted by the 3-d model are visually similar to a plate buckling mode when the aspect-ratio is small (see figure 5.3c for h/v = .1): the bistrips deforms essentially by bending and the displacement is mainly out-of-plane.

A quantitative agreement between the plate and the 3-d models can in fact be obtained: in figure 5.3(b), we check that the results of the 3-d hyper-elastic model converge towards the predictions of the plate model in the limit of a small aspect-ratio, h/v → 0. The reference curve shows the critical wavenumber of the plate as a function of the mismatch strain (thick curve): this plot uses the natural variables of the plate model, namely the dimensionless mismatch strain δ and the rescaled wavenumber q c defined in §4. The thin curves correspond to the results of the 3-d Gent model, for different aspect-ratios. The latter set of curves are plotted as follows. First, the eigenvalue analysis of §5.4 yields data-points (p, λ, q c ). Second, the stretches (p, λ) characterizing the loading in a finite-elasticity context are converted into strains (δ, ) relevant to the small-strain limit by identifying the stretches in figures 4.1 and 5.1(3) as 1 + + δ 2 = λ p and 1 + -δ 2 = λ. Third, the dimensionless mismatch strain δ and mean imposed strain are calculated using (4.4e). As we work in units such that v = 1, the dimensionless wavenumber q c = q v is identical to the wavenumber q c predicted by the 3-d model. We use an incompressible plate model (ν = .5) for comparison, consistent with the fact that we use an almost incompressible Gent model.

Small-wavenumber expansion

We now apply the small-wavenumber expansion presented in §2.11 to the 3-d hyper-elastic model: this will allow us to predict whether the buckling mode is macroscopic or microscopic without carrying out a computationally intensive eigenvalue analysis, and will provide a more accurate value for the critical mismatch strain at the macroscopic to microscopic transition. This section is independent from the previous ones, and can be skipped in a first reading.

The first step is to cast the bifurcation equation (5.10b) into the canonical form (2.3a) that has no linear term in q. This is a tricky step. The idea is to duplicate and arrange the degrees of freedom into a new unknown vector Ξ: in block notation,

Ξ =                 λ -1 ξ q -1 ξ ⊥ ξ q λ -1 ξ ⊥                 .
(5.13a) Each half of the vector Ξ obeys the ordering convention in (5.10a), i.e. the first two thirds of each half encodes cross-sectional displacements and the last third of each half encodes axial displacement. In this block-vector notation, the double-struck line splits the generalized degrees of freedom into two halves, while the single-struck lines separate the the cross-sectional and the axial degrees of freedom. We carry on with this convention in the following. With the block matrices A λ and B λ defined as (5.13b) and the sub-blocks defined themselves using block notation as (5.13c) one can cast the eigenvalue problem (5.10b) governing bifurcation into an equivalent form that matches up the canonical form proposed earlier in (2.3a),

A λ =    0 A λ A T λ 0    , B λ =    0 B λ B T λ 0    ,
A λ =    λ k λ λ 2 c T λ 0 λ k ⊥ λ    , B λ =   λ m λ 0 c λ λ m ⊥ λ   ,
∀ Ξ, Ξ • A λ • Ξ + q 2 Ξ • B λ • Ξ = 0. (5.13d)
Here, Ξ is an arbitrary vector having the same dimension as Ξ in (5.13a). Note that it is straightforward to assemble the block matrices A λ and B λ once the base solution ϕ λ 0 is known, as this involves rearranging the blocks k λ , . . . , m ⊥ λ calculated earlier in (5.11). The small-wavenumber expansion can now be carried out by applying the recipe of §2.11 to (5.13d):

1. the matrices A λ and B λ are obviously symmetric. The null space ker A λ is worked out in Appendix C: it is spanned by the 8 following vectors,

      0       . (5.14)
Each of these vectors has the same length and ordering convention as the vector Ξ in (5.13a). The eight vectors represent special types of adjacent equilibria: six of them are rigid-body modes (3 translations: T * 1 , T * 2 , and T 4 ; and 3 rotations: T * 3 , T 1 , T 2 ), and two are homogeneous infinitesimal transformations: a stretching mode (T * 4 ) and a twisting mode (T 3 ), see the appendix for details. As a result, ker A λ is of dimension 8, and is indeed independent of λ; 2. in the basis (5.14), the restriction B * (λ) of B λ to ker A is a symmetric 8 × 8 matrix which is antidiagonal: using block-matrix notation,

B * (λ) = 0 b λ b T λ 0 .
(5.15)

The entries of b λ can be calculated as, for instance, b λ 14 = B * 15 (λ) = T * 1 • B T λ • T 4 . Note that the 4 × 4 matrix b λ is non-symmetric. The critical load λ (0) is found by solving det B * (λ (0) ) = 0: in view of the block form of B * , this is equivalent to finding the root λ = λ (0) of the (simpler) equation det b λ = 0; 3. at the critical value λ = λ (0) , the non-symmetric matrix b

λ = b (0) is singular. Let β = β 1 , • • • , β 4
denote a null vector of b (0) , and

β † = β † 1 , • • • , β † 4 a null vector of its transpose b T (0)
. In view of (5.15), the null space of B * (λ (0) ) is of dimension n * = 2 and is spanned by . Each of these vectors has eight entries, which are the coordinates of a null vector in the basis (5.14): the null vectors can be spelled out in the original space by multiplying the coordinates by the basis vectors. This shows that the following two vectors form a basis of ker B * (λ (0) ):

      V 0       ,       0 V *       , (5.16) where V = 4 i=1 β i T i and V * = 4 i=1 β * i T * i .
By equation (2.12a), the limit Ξ (0) of the bifurcation mode for q → 0 lives in the null space of B * (λ (0) ): for some coefficients α (0) 1 and α

(0) 2 , Ξ (0) = α (0) 1       V 0       + α (0) 2       0 V *       .
(5.17) 4. a particular solution solution Ξ (1) p of the eigenvalue problem (2.9*) at order q 2 is found as 0) , where A -1 (0) is a pseudo-inverse of the singular operator A (0) ; 5. the solvability condition (2.13) at order q 4 can then be written as

Ξ (1) p = -A -1 (0) • B (0) • Ξ (
∀ α1 , α2 ,       α1 V α2 V *       • -B (0) • A -1 (0) • B (0) + λ (1) B (0) •       α (0) 1 V α (0) 2 V *       = 0.
Expanding the matrices inside the square brackets using their block form, this can be rewritten as

∀ α1 , α2 , α1 α2 • -a 1 + a 2 λ (1) 0 1 1 0 • α (0) 1 α (0) 2 = 0,
where

a 1 = V • B (0) • A -1 (0) • B (0) • V * , a 2 = V • B (0) • V * and A -1 (
0) denotes a pseudo-inverse of the block A (0) in (5.13b); 6. the value of λ (1) is found by the condition that the matrix in square brackets in the equation above is singular, viz. λ (1) = a 1 /a 2 :

λ (1) = V • B (0) • A -1 (0) • B (0) • V * V • B (0) • V * .
(5.18)

Note the similarity with the equation (2.14) valid for n * = 1 (which is not directly applicable as n * = 2 here).

Practically, the quantities λ (0) and λ (1) are calculated as follows. For any value of the mean applied stretch λ, the homogeneous solution ϕ 0 (λ; x, ỹ) is calculated by solving a 2-d non-linear elasticity problem using the finite-element method, as explained in §5.3. The matrices k λ , . . . , m ⊥ λ appearing in (5.11) are assembled as described in Appendix B.2. These blocks are arranged into the matrices A λ and B λ , see equation (5.13c). The 4 × 4 matrix b λ is calculated as explained below equation (5.15). Iterating over λ, numerical root-finding is used to find the root λ = λ (0) of det b λ = 0. The corresponding right and left null vectors β and β * are calculated numerically, as well as V and V * appearing immediately after equation (5.16). A pseudo-inverse A -1 (0) of A (0) is calculated. The operator B (0) = dB λ dλ λ (0) appearing in the denominator of (5.18) is estimated numerically using finite differences in λ. Finally, λ (1) is calculated from equation (5.18).

This method is first applied to the example problem from figure 5.3(a), using the same aspect-ratio h/v = 0.1 and the same set of material parameters as for the direct method: in the figure, the dotted parabolas λ = λ (0) + q 2 λ (1) are the prediction of the small-wavenumber expansion, for different values of the mismatch strain p. In the neighborhood of their intersection with the q = 0 axis, the parabolas provide a good approximation of the complete curves of marginal stability obtained earlier by the direct method.

For this particular example, λ (1) , as calculated by (5.18), goes from positive for p = 1.11 to negative for p = 1.056: this is in accord with the macroscopic to microscopic transition observed at p * ≈ 1.1 using the full eigenvalue analysis (direct method), see equation (5.12). By doing a linear interpolation of the values of λ (1) just calculated from (5.18) at p = 1.056 and p = 1.11, one can obtain an accurate estimate of p * as p * = 1.090.

(5.19)

When converted into the dimensionless mismatch strain of the plate model, as explained as the end of §5.6, this yields δ = 89.6. This is close to the value δ = 88.25 calculated by the plate model in (4.12): this confirms that the plate model agrees with the 3-d model in the limit h/v → 0; the small difference between the two values is due to the fact that the aspect-ratio is small but finite in the 3-d model (h/v = 0.1). This procedure for calculating the critical value of the mismatch stretch p * was repeated for different aspect-ratios: these values of p * were converted into values of δ, and then displayed as stars in figure 5.3(b). Each of these stars falls exactly at the point where the curve of marginal stability found by the direct method merges with the axis q c = 0. This confirms that the predictions of the small-wavelength expansion are valid.

Summary

In this section, we modelled the bistrip experiment using a 3-d hyper-elastic model. By solving a nonlinear elasticity problem in 2-d using the finite-element method, we calculated a family of base solutions that are cylindrically invariant. Using numerical eigenvalue analysis, we obtained the curves of marginal stability for different values of the mismatch stretch p (direct method); we reached the same conclusion as with the double-beam and plate models: the first bifurcation mode goes from macroscopic to microscopic when the mismatch stretch goes beyond a critical value, p = p * . The predictions of the 3-d hyper-elastic model agree with those of the plate model when the aspect ratio is small, h v. We also applied the general small-wavelength expansion to calculate the critical mismatch strain p * corresponding to the macroscopic-to-microscopic transition, for different values of the aspect-ratio h/v. The values of p * obtained in this way agree with those calculated by the direct method. The small-wavelength expansion has two important benefits over the direct method. It provides more accurate values of p * , and it is much less intensive computationally: it requires only standard linear algebra operations on the matrices K λ , C λ , M λ that are known in terms of the invariant solution ϕ 0 . By contrast, the direct method requires solving a polynomial eigenvalue problem for each value of λ, and then iterating over λ to bracket the critical stretch λ c and the wavenumber of the first bifurcating mode q c . This determines whether the system bifurcates towards a microscopic or a macroscopic mode. Then a second level of iteration over the pre-stretch p is required to bracket the critical value p for a particular value of the aspect-ratio h/v.

In the limit of a weak pre-strain, p → 1, the predictions of the 3-d model studied in this section are equivalent with those of an Euler-Bernoulli rod model having a natural curvature consistent with the prestrain distribution. This can be shown by an argument similar to that given earlier in §4.5 for a plate; it is more technical, however, and is not given here.

Discussion, conclusion and perspectives

In previous work on thin elastic structures, a weakly inhomogeneous pre-strain has been taken into account by modifying the reference metric or the natural curvature of classical rod or plate models [START_REF] Goriely | Spontaneous helix hand reversal and tendril perversion in climbing plants[END_REF][START_REF] Liu | Structural transition from helices to hemihelices[END_REF]. The present work shows that this is not possible in the general case of a finite, non-uniform pre-strain: as evident from the analysis of bistrips, large pre-strain can lead to buckling modes that violate the kinematic assumptions underlying classical rod models. Accordingly, we have approached the stability of a pre-strained rod-like elastic solids based on structural models with deformable cross-sections: we used different structural models featuring different degrees of sophistication, all of which yielded qualitatively similar results. We found that the first bifurcation mode switches from long-wavelength for weak pre-stress (q c = 0), to short-wavelength for larger pre-stress (q c = 0). The transition is governed by a dimensionless parameter comparing the antagonistic effects of pre-stress (which favors short wavelengths) and of the stiffness of the cross-section (which favors long wavelengths). For large or small pre-stress, our models are asymptotically equivalent with the buckling of an Euler-Bernoulli beam with natural curvature, and with a strut on an elastic foundation, respectively. Drawing an analogy with Koiter's post-bifurcation analysis, we have proposed a small-wavenumber expansion that captures the long-wavelength to shortwavelength transition. This method is especially powerful when applied to the full 3-d model ( §5): it involves elementary linear algebra in terms of the unbuckled solution, and removes the need for a computationally extensive eigenvalue analysis.

Our initial motivation for this work was to explain the selection of the number of perversions in a pre-stressed bistrip based on a static, linear bifurcation analysis; this goal has been achieved. The 3-d model studied in §5 makes use of the same constitutive law and modelling assumptions as the finite-element simulations of [START_REF] Liu | Structural transition from helices to hemihelices[END_REF], which they already compared to the bistrip experiments: there is no need for us to carry out a detailed comparison with the experiments again. We simply note here that we found the critical wavenumber q c to be an increasing function of the dimensionless pre-strain, and that this accounts for the increase of the number of perversions with the amount of pre-strain and with the inverse crosssection aspect-ratio v/h, as seen in the experiments. As our approach makes uses of Fourier analysis in the axial direction, the macroscopic-to-microscopic transition manifests itself as a bifurcation on the critical wavenumber q c . To extend our work to a finite length L, the dispersion relation derived in this paper needs be combined with the relevant boundary conditions which, typically, provide a quantification condition on the wavenumber q. For very slender rods, L h and L v, the wavenumber is quantified in vanishingly small increments and this warrants consistency with the case L = ∞ addressed in the present work.

We emphasize that the present work is concerned with a linear bifurcation analysis, i.e. our equations are linear with respect to the amplitude of buckling but non-linear with respect to the wavenumber q. In future work, it would be interesting to consider the case of a finite buckling amplitude as well: this would allow one to confirm that the bifurcation is super-critical, and to compute the amplitude along the post-bifurcated branch; other typically non-linear effects such as the coalescence and subdivision of perversions may also be addressed [START_REF] Liu | Emergent perversions in the buckling of heterogeneous elastic strips[END_REF]. Another interesting extension would be to push the small-wavelength extension of §2 to the next order: this would allow one to derive an 'amplitude equation' for the critical wavenumber of the form q 4 c -α (p -p * ) q 2 c = 0 near the critical pre-strain p * , corresponding to the asymptotic law q c ∝ (p -p * ) 1/2 seen e.g. in the numerical results of figure 5.3(b) at the transition.

In our work, the macroscopic and the microscopic buckling modes originate from a single mode whose wavenumber bifurcates: this is entirely different from the problem of mode jumping where different buckling modes interact non-linearly [START_REF] Stein | The phenomenon of change in buckle pattern in elastic structures[END_REF][START_REF] Schaeffer | Boundary conditions and mode jumping in the buckling of a rectangular plate[END_REF][START_REF] Holder | Boundary conditions and mode jumping in the von Kármán equations[END_REF][START_REF] Everall | Mode jumping in the buckling of struts and plates: a comparative study[END_REF][START_REF] Chien | Mode jumping in the von Kármán equations[END_REF], and from other problems involving an interaction between different macroscopic and microscopic modes [START_REF] Bai | Mode interaction in thin-walled I-section struts with semi-rigid flange-web joints[END_REF].

In the small-wavenumber expansion we assumed that the curves of marginal stability are always of one of the two types shown in figure 2.1. This assumption turned out to be correct for all the examples that we studied and, as a result, the criterion λ (1) = 0 (or (1) = 0 in §3-4) successfully captured the macroscopicto-microscopic transition. If, however, the marginal curves have at least three critical points (λ c , q c ) where dλ/dq = 0, different types of bifurcations may take place and the criterion λ (1) = 0 may not be applicable (this is similar to the well-known limitation of a linear bifurcation analysis, which overestimates the load at which buckled solutions can appear by a sub-critical pitchfork bifurcation).

Motivated by the bistrip experiments, we have considered a piecewise constant pre-strain, see (4.1a) and (5.1). It is straightforward to extend the present work to an arbitrary distribution of pre-strain in the cross-section. As an illustration, let us consider the plate model, now with a pre-stress depending linearly on the transverse coordinate x: we replace (4.1a) with N 0 zz ( ; x) = Eh + 3 2 x L -1 2 δ . In this formula, the numerical factor 3/2 has been chosen such that the residual bending moment is still given by equation (4.14): this warrants that the mismatch parameter δ of the new (linear) pre-stress distribution is directly comparable with that of the older (piecewise constant) pre-stress distribution. Instead of carrying out the buckling analysis by the direct method (shooting method, §4.2), we use the small-wavenumber expansion: it gives the same results more easily. Adapting the analysis of §4.4, we find that equation (4.8) is where H 1 is the incremental transformation gradient defined below (5.7). For the Gent model, the full expression of the moduli B 1 , B 2 and B 3 is

B 1 = 2 µ J m det G (J m -I 0 c + 3) 2 B 2 = 2 µ -(K -2µ Jm )J 0 (J 0 -1) det G B 3 = K -2µ Jm J 0 2 J 0 -1 det G ,
where I 0 c = tr C 0 and J 0 = det C 0 are the deformation invariants evaluated in the invariant solution ϕ λ 0 .

Appendix B.2. Quadratic forms entering in the eigenvalue problem Let us return to the weak form of the linearized equilibrium (5.7). We do not assume that the real and virtual displacements ϕ 1 (x, ỹ, z) and ϕ(x, ỹ, z) are harmonic in the axial coordinate z, as we did earlier in (5.9); instead, we consider an arbitrary dependence on z. This will allow us to derive the representation of rigid-body motions in Appendix C below. The linearized equilibrium (5.9) depends on the first gradients of ϕ 1 and ϕ along both the cross-sectional directions x and ỹ and along the axial direction z. By grouping the terms depending on the total order of derivation with respect to z, one can rewrite the linearized equilibrium as

∀ ϕ, 1 L L 0 K λ ( ϕ, ϕ 1 ) + C λ ( ϕ, ϕ 1,z ) + C λ (ϕ 1 , ϕ ,z ) + M λ ( ϕ ,z , ϕ 1,z ) dz = 0, (B.2a)
where the symbol z appearing in a subscript after a comma denotes an axial derivative ∂/∂ z. The bilinear forms K λ , C λ and M λ , which are spelled out below, depend on the pre-stretch λ through the base solution ϕ λ 0 , and contain gradients of the displacement with respect to the cross-sectional directions x and ỹ, noted as grad, but not with respect to the axial direction z. Stated differently, cross-sectional gradients ∂/∂ x and ∂/∂ ỹ (denoted by grad) are implicitly contained in the operators in equation above, but the axial gradients ∂/∂ z appear explicitly in the arguments of the operators.

Next, the real and virtual displacements are decomposed into a cross-sectional projection ϕ 1 and ϕ belonging to the (xy) plane, and an axial component as

ϕ 1 (x, ỹ, z) = ϕ 1 (x, ỹ, z) + ϕ ⊥ 1 (x, ỹ, z) e Z , ϕ(x, ỹ, z) = ϕ (x, ỹ, z) + ϕ ⊥ (x, ỹ, z) e Z .
The operators in (B.2a) are found by identification with (5.7): the way they couple the cross-sectional and axial components is dictated by the symmetry of the invariant solution as follows,

K λ ( ϕ, ϕ 1 ) = k λ ( ϕ , ϕ 1 ) + k ⊥ λ ( ϕ ⊥ , ϕ ⊥ 1 ), (B.2b) C λ ( ϕ, ϕ 1,z ) = c 1 λ ( ϕ ⊥ , ϕ 1,z ) + c 2 λ ( ϕ , ϕ ⊥ 1,z ), (B.2c) M λ ( ϕ ,z , ϕ 1,z ) = m λ ( ϕ ,z , ϕ 1,z ) + m ⊥ λ ( ϕ ⊥ ,z , ϕ ⊥ 1,z ). (B.2d)
From Appendix B.1, the detailed expressions of the different 'blocks' are, for the Gent model,

k λ ( ϕ , ϕ 1 ) = D r2 (gradϕ 1 • Σ 0 ) : grad ϕ + r4 B 1 F 0 : gradϕ 1 F 0 : grad ϕ • • • + r4 B 2 (F 0 ) T • gradϕ 1 • (C 0 ) -1 : (F 0 ) T • grad ϕ • (C 0 ) -1 • • • + r4 B 3 (C 0 ) -1 : ((F 0 ) T • gradϕ 1 ) (C 0 ) -1 : ((F 0 ) T • grad ϕ ) dx dỹ, (B.2e) k ⊥ λ ( ϕ ⊥ , ϕ ⊥ 1 ) = D r2 Σ 0 • gradϕ ⊥ 1 .grad ϕ ⊥ + r2 B 2 2 grad ϕ ⊥ . (C 0 ) -1 • gradϕ ⊥ 1 dx dỹ, (B.2f) c 1 λ ( ϕ ⊥ , ϕ 1,z ) = D r2 B 2 2λ grad ϕ ⊥ .((F 0 ) -1 • ϕ 1,z ) dx dỹ, (B.2g) c 2 λ ( ϕ , ϕ ⊥ 1,z = D r2 p2 B 1 λ F 0 : ϕ ⊥ ,1z grad ϕ + r2 B 3 λ (F 0 ) -T : ϕ ⊥ 1,z grad ϕ 1 dx dỹ, (B.2h) m λ ( ϕ ,z , ϕ 1,z ) = D p2 Σ ⊥ 0 ϕ 1,z . ϕ ,z + B 2 λ 2 1 2 ϕ 1,z • ϕ ,z dx dỹ, (B.2i) m ⊥ λ ( ϕ ⊥ ,z , ϕ ⊥ 1,z ) = D p2 Σ ⊥ 0 ϕ ⊥ 1,z ϕ ⊥ ,z + p4 B 1 λ 2 ϕ ⊥ 1,z ϕ ⊥ ,z + B 2 λ 2 ϕ ⊥ 1,z ϕ ⊥ ,z + B 3 λ 2 ϕ ⊥ 1,z ϕ ⊥ ,z dx dỹ. (B.2j)
Note that the operators K λ and M λ , which are associated with an even total order of axial differentiation ∂/∂ z, couple cross-sectional real and virtual displacements on one hand, and axial real and virtual displacement on the other hand; by contrast, the operator C λ which is associated with a total order of axial differentiation ∂/∂ z equal to one (odd order), couples a cross-sectional virtual displacement with an axial real displacement and vice versa.

As we work with a set of material parameters that make the solid nearly incompressible, we approximate det G ≈ 1 in the linear bifurcation analysis.

Appendix B.3. Connection with the Fourier form of the operators

With the aim to derive the Fourier form of the eigenvalue problem announced in (5.11), we now consider displacements that are pure Fourier modes, i.e. that depend harmonically on z, as in (5.9): ϕ 1 (x, ỹ, z) = ξ (x, ỹ) e i q z ϕ ⊥ 1 (x, ỹ, z) = i ξ ⊥ (x, ỹ) e i q z ϕ (x, ỹ, z) = ξ (x, ỹ) e i q z ϕ ⊥ (x, ỹ, z) = i ξ⊥ (x, ỹ) e i q z

For large L, the averages over the axial variable z in equation (B.2a) can be simplified with the help of the identity 1 L L 0 S e i q z : T e i q z dz = 1 2 S : T † (for q = 0) for any tensors S and T of complex amplitudes. Here, denotes the real part, and T † is the complex conjugate of T . When this identity is applied to (B.2a), the variable z disappears, the wavenumber q appears, and one arrives at the eigenproblem announced in (5.10b): ∀ ξ, ξ • K λ + q C λ + q 2 M λ • ξ 1 = 0, using the block notations of §5.4. In this equation, the operators K λ and M λ are the discrete form of the operators K λ and M λ that appeared in (B.2b) and (B.2d), while C λ is a discrete symmetric matrix coupling cross-sectional and axial displacements: For the 3-d prismatic solid, the linearized equilibrium equations (5.7) have obvious solutions corresponding to rigid-body displacements, as well as to uniform twist and uniform (increment of) axial stretching.

Here, we seek a representation of these modes in our notation. As we shall show, all these modes turn out to be null vectors of the stiffness matrix A λ , and as such they play a key role in the bifurcation analysis. By analyzing these special modes, our goal is to derive the expressions of the eight vectors forming a basis of A λ , see equation (5.14).

We consider a base solution x 0 = x + ϕ 0 (x, ỹ) that is invariant in the axial direction. We start from the linearized equilibrium in (B.2a), which is equivalent to (5.7): using the block structure of the operators given in (B.2b-B.2d), and carrying out an integration by part, we write the linearized equilibrium as ∀ ϕ (x, ỹ, z), ϕ ⊥ (x, ỹ, z) ,

1 L L 0 k λ ( ϕ , ϕ 1 ) + k ⊥ λ ( ϕ ⊥ , ϕ ⊥ 1 ) -m λ ( ϕ , ϕ 1,z z ) -m ⊥ λ ( ϕ ⊥ , ϕ ⊥ 1,z z ) -((c 1 λ ) T -c 2 λ )( ϕ , ϕ ⊥ 1,z ) + 1 λ -(c 2 λ ) T )( ϕ ⊥ , ϕ 1,z ) dz = 0.
The boundary terms coming from the integration by parts have been as the prismatic solid is considered L → ∞.

Thanks to the integration by part, the derivatives of the virtual displacement with respect to the axial variable z have been eliminated in equation above. This allows one to remove both the integral and the virtual displacement: we arrive at the linearized equilibrium in strong form,

k λ 0 0 k ⊥ λ • ϕ 1 + 0 -c T λ c λ 0 • ϕ 1,z - m λ 0 0 m ⊥ λ • ϕ 1,z z = 0 (C.1)
where we have defined the operator c λ = c 1 λ -(c 2 λ ) T defined below (B.3). In the following, we review 8 particular solutions ϕ 1 to this equation, namely 3 rigid-body translations, 3 rigid-body rotations, and 2 non-rigid modes (corresponding to uniform stretching, and to uniform twisting, respectively).

Appendix C.1. Modes having a cross-sectional displacement independent of z The modes whose cross-sectional displacement is independent of z are analyzed first.

• A rigid-body translation along the axis x corresponds to an incremental displacement ϕ 1 (x, ỹ, z) = e X .

This ϕ 1 has to be a solution of (C.1): inserting into this equation, we find

k λ • e X = 0,
where we have identified the unit vector e X with the constant function that takes on the value e X . Using the block-vector notations introduced in §5.7, see equation (5.13c) we rewrite this as:

A λ • T * 1 = 0, where T * 1 =       e X 0       =⇒       0 T * 1       ∈ ker A λ .
The block matrix A λ has been defined in (5.13b). We have just shown that the first vector listed in equation (5.14) is indeed a null vector of A λ . • Let us consider an infinitesimal rotation about the axis Z. The displacement corresponding to this infinitesimal rotation is contained in the plane of the cross-section. an arbitrary normalization, it writes ϕ 1 (x, ỹ, z) = e Z × x 0 (x, ỹ), where × denotes the cross product of two vectors and x 0 = (x, ỹ) + ϕ 0 (λ; x, ỹ) is the projection onto the cross-section of a material point in the invariant solution labelled (3) in figure 5.1: note that the rotation operates on this configuration and not on (1) or (2).

Introducing the notation ρ z (x, ỹ) = e Z × x 0 (x, ỹ), and repeating the same argument as earlier, we find that the rigid-body rotation is represented by a vector that lives in the null space of A λ ,

T * 3 =       ρ z (x, ỹ) 0       =⇒       0 T * 3       ∈ ker A λ .
• Finally, we consider a uniform increment of stretch along the axis Z. Differentiating the invariant solution with respect to the stretch λ and multiplying by a conventional factor λ, we find the corresponding mode as ϕ 1 (x, ỹ, z) = u str (x, ỹ) + λ z e Z , where the increment of cross-sectional displacement reads u str (x, ỹ) = λ dϕ 0 (λ;x,ỹ) dλ

. Inserting into (C.1), this yields k λ • u str (x, ỹ) -λ c T λ • e Z = 0. This equation is well-known in the analysis of the stretching of thin rods: this is the partial differential equation in the 2-d cross-section, which must be solved for the cross-sectional displacement u str to capture the incremental deformation of the cross-sections by Poisson's effect. Using notations introduced in equation (5.13c), §5.7, and observing that k ⊥ λ • e Z = 0 as a consequence of the symmetries-as confirmed by the detailed expression of k ⊥ λ in (B.2f)-, we rewrite this as

A λ • T * 4 = 0, where T * 4 =       -u str (x, ỹ) e Z       =⇒       0 T * 4       ∈ ker A λ .
Appendix C.2. Modes having an axial displacement independent of z We review the four remaining modes. Their axial displacement happens to be independent of z.

• The rigid-body translation along Z is described by a constant perturbation along the axis, ϕ 1 (x, ỹ, z) = e Z . Inserting into (C.1), we rewrite this as k ⊥ λ • e Z = 0, 39 which we can rewrite using the notation of §5.7 as

A T λ • T 1 = 0, where T 1 =       0 e Z       =⇒       T 1 0       ∈ ker A λ .
• The infinitesimal rigid-body rotation about the transverse axis X writes ϕ 1 (x, ỹ, z) = e X × x 0 (x, ỹ, z) =

-λ z e Y + ρ ⊥ x (x, ỹ) e Z where ρ ⊥ x (x, ỹ) = e X × x 0 (x, ỹ) • e Z = e Y • x 0 (x, ỹ). Inserting into (C.1) and observing that k λ • e Y = 0, we have

-λ c λ • e Y + k ⊥ λ • ρ ⊥ x = 0.
Using the notation of §5.7, we rewrite this as

A T λ • T 2 = 0, where T 2 =       e Y -ρ ⊥ x       =⇒       T 2 0       ∈ ker A λ .
• An infinitesimal rigid-body rotation about the second transverse axis Y is treated similarly: with • Finally, the mode corresponding to a uniform increment of twist writes ϕ 1 (x, ỹ, z) = λ z ρ z (x, ỹ) + u ⊥ warp (x, ỹ) e Z , where u ⊥ warp is known as the warping function and ρ z is the infinitesimal rotation about the axis introduced earlier in the analysis of rigid-body rotations. Note that the rotation about the axis z is now multiplied by an angle λ z proportional to the axial coordinate (uniform twist). Inserting the expression of ϕ 1 into (C.1) , we obtain k ⊥ λ • u ⊥ warp + λ c λ • ρ z = 0, which is the discrete version of the classical partial differential equation that yields the warping function u ⊥ warp in the theory of thin elastic rods. Adapting the argument used for the analysis of rigid-body rotations about the axes X and ỹ, this can rewritten be as

      T 4 0       ∈ ker A λ , where T 4 =       ρ z u ⊥ warp       .
In total, we have found eight infinitesimal modes in the null space of A λ . We can offer no proof for the fact that these eight vectors span the entire null space of A λ , but we have checked numerically that ker A λ is indeed of dimension eight in our numerical examples.

Figure 1 . 1 :

 11 Figure1.1: Principle of the experiment, and sketch of the main experimental and numerical findings of[START_REF] Huang | Spontaneous and deterministic three-dimensional curling of pre-strained elastomeric bi-strips[END_REF];[START_REF] Liu | Structural transition from helices to hemihelices[END_REF]. (a-c) Preparation of the bistrip with inhomogeneous pre-stress in the cross-section; (d) the terminal tensile load is decreased and the bistrip buckles. (e) Incipient buckling mode (top) and post-buckled solution (bottom), for macroscopic (left) and microscopic (right) buckling. The position of the perversions in the post-buckled solutions is determined by the crests and valleys of the incipient buckling pattern (star symbols): a macroscopic buckling mode yields at most one perversion in the entire strip (left) while a microscopic buckling mode yields two perversions per wavelength (right).

Figure 2 . 1 :

 21 Figure2.1: Two buckling scenarios for a prismatic solid with non-uniform pre-stress in the cross-section: typical curves of marginal stability predicted by the generic equation (2.3a). The mean imposed stretch λ is progressively decreased (meaning that the diagram has to be read from right to left) and the homogeneous solution bifurcates along the solid curve. (a) Case of a microscopic buckling instability: the system first bifurcates at λ = λc and a wavenumber qc = 0 is selected. We denote by λ (0) the stretch such that a linear mode with an infinite wavelength is available near the fundamental branch: λ (0) < λc for microscopic buckling. (b) Case of a macroscopic buckling instability: a zero wavenumber is selected, qc = 0, and λ (0) = λc.

  Figure3.1: The double-beam model: δ is a mismatch strain, and an imposed mean strain (imposed by the displacement of the the remote endpoints). For the sake of legibility, an extensional pre-stretch is shown, + δ/2 > -δ/2 > 0; in reality the structure does not becomes unstable unless of the pre-stretch is compressive, -δ/2 < 0.

Figure 3

 3 Figure3.2: Bifurcation analysis of the double-beam. (a) Curves of marginal stability as a function of the imposed average strain for different values of the dimensionless mismatch strain δ = 1, 2, √ 8, 3, 4, 5 (from dark to bright). The orange curve for δ = δ * = √ 8 separates the regime of macroscopic buckling from the regime microscopic buckling. The dotted parabola is the result of the asymptotic analysis in equation (3.9) for δ = 5 and is checked to be consistent with the direct method (grey curve). The 'stable' and 'unstable' keywords indicate the expected region of stability and instability of the unbuckled solution. The critical values for the imposed average strain c and the associated wavenumber q c are denoted by black disks. (b) Wavenumber q c of the first unstable mode as functions of dimensionless mismatch strain δ, featuring the macroscopic to microscopic transition at the critical value of δ * . The dashed grey line represents the limit value of q c when δ → ∞; this

  the mismatch strain δ > δ * : the double-beam model successfully reproduces the macroscopic-to-microscopic transition reported in the experiments of Huang et al. (2012) and Liu et al. (2014) and in their stability analysis based on a 3-d finite-element model.

  Figure3.3: Convergence of the bifurcation modes of the bi-rod model. (a) Equivalence with the Euler beam in the limit where the dimensionless mismatch strain is small, δ δ * : the curves of marginal stability in the plane ( , q) collapse onto

Figure 4

 4 Figure 4.1: A plate model with piecewise constant pre-strain: the homogeneous solution shown here can become unstable when the lowest pre-stress N 0,I zz is compressive enough.

  .4e) Equation (4.4a) is an eigenvalue problem of the general type (2.2) when the operators are identified as d 0 λ = a, d 2 λ = b and d 4 λ = c, and d 1 λ = d 3 λ = 0 by symmetry.

Figure 4

 4 Figure 4.2: Bifurcation from the planar solution of the plate model, with ν = .5. (a) Critical wavenumber as a function of mean imposed dimensionless strain for different values of the dimensionless mismatch strain δ.The solid curves are the predictions of the direct method (4.6); the dotted curve is the prediction of the small-wavenumber expansion (4.10b) for δ = 800. The labels 'stable' and 'unstable' indicate the expected region of stability and instability of the planar solution: this would need to be confirmed by a proper stability analysis. The critical values for the imposed average strain c and the associated wavenumber q c are denoted by black disks. (b) Wavenumber, (c) measure of the curvature of the cross-section, and (d) eigenvector corresponding to the first bifurcating mode = c, for different values of the dimensionless mismatch strain δ. The points labelled A, B and C in different plots correspond to different levels of mismatch strain, δ = 50, 200, 400, respectively.

Figure 4

 4 Figure 4.3: Convergence of the bifurcation modes of the plate model for small and large dimensionless mismatch strain δ, with ν = .5. (a) The case of small mismatch strain (δ δ * ) is equivalent to a rectangular plate with uniform pre-stress. (b) The case of large mismatch strain (δ δ * ) is equivalent to a half-plate (I) with one free and one clamped edge.

  Figure 4.4: Bifurcation diagram for the planar solution of the plate model, for ν = .5. (a) Phase diagram: results of the shooting method from figure4.2(a) (stars, disks and thin black curve) and predictions of the small-wavenumber expansion (thick green curve and red domain). The bifurcation mode is macroscopic when the green and black curves coincide, which is also when the green curve is outside the red region; it is microscopic otherwise. The prediction

  parameter δ are denoted by the keywords 'macro' and 'micro' below the phase diagram in figure4.4(a): the small-wavenumber expansion predicts that the first bifurcation mode is macroscopic if |δ| < δ *

Figure 5

 5 Figure 5.1: The 3D hyper-elastic model: (1) natural configuration, (2) pre-stretched configuration, (3) cylindrically invariant solution with mismatch stretch p and mean imposed stretch λ, (4) buckled configuration.

Figure 5

 5 Figure 5.2: Pre-stress in the homogeneous solution: axial pre-stress Σ ⊥ 0 (λ; x, ỹ) (left) and norm of cross-sectional stress |Σ 0 (λ; x, ỹ))| (right), see equation (5.6). Aspect ratios h/v = 0.1, 0.2, 0.4 (from top to bottom); material parameters: µ = 1, K = 10, Jm = 100; mismatch stretch p = 1.11; mean imposed stretch λ = .9. The axial stress (left) is nearly uniform in each half of the cross-section: it is compressive in region I and tensile in region II (because of the compressive axial stress, this cylindrically invariant configuration is actually unstable, see figure 5.3a). The cross-sectional stress (right) arises from geometric incompatibility between the two domains, and is limited to a neighborhood of their interface; it is much smaller than the axial stress (note the widely different scales in the color bars).

Figure 5

 5 Figure 5.3: Bifurcation analysis of the 3D hyper-elastic model, for material parameters µ = 1, K = 10, Jm = 100 and typical mesh size e = 0.02. (a) Curves of marginal stability for an aspect ratio h/v = 0.1 and for different values of the mismatch stretch p: results of the eigenvalue analysis (direct method, §5.4-5.5, symbols) and comparison with the small-wavenumber expansion (dotted parabolas, §5.7). The critical values for the mean stretch λc and the associated wavenumber qc are denoted by black disks. (b) Numerical test of convergence towards the the plate model in the limit of a small aspect-ratio h/v → 0: the critical wavenumber predicted by the 3D model using the direct method (dots connected by thin curves) converges to that predicted by the plate model (thick curve). For all aspect-ratios, the small-wavenumber expansion (star symbols on the horizontal axis) correctly predicts the critical mismatch strain δ * at which the bifurcation mode switches from macroscopic to

  C λ is obtained by discretizing the operatorC λ ( ξ, ξ) = c 1 λ ( ξ⊥ , ξ ) -(c 2 ) T λ ( ξ⊥ , ξ ) + (c 1 ) T λ ( ξ , ξ ⊥ ) -c 2 λ ( ξ , ξ ⊥ ). (B.3)Therefore, the lower-left block c λ appearing in (5.11) is associated with the continuous operator c λ = c 1 λ -(c 2 λ ) T . The symmetry and the diagonal or anti-diagonal block-structure of the operators stated in equation (5.11) are a consequence of their definitions (B.2b-B.2j).Appendix C. Rigid and neutral modes of the 3-d operators

•

  A similar argument shows that the rigid-body translation along the axis Y , ϕ 1 (x, ỹ, z) = e Y , is represented by the second null vector of A λ listed in equation (5.14), the vector e Y appearing in the first block of T * 2 is the function taking on the constant value e Y , by a slight abuse of notation.

  ϕ 1 (x, ỹ, z) = e Y ×x (x, ỹ, z) = λ z e Y +ρ ⊥ y (x, ỹ) Z and ρ ⊥ y (x, ỹ) = e Y × x 0 (x, ỹ) •e Z = -e X •x 0 (x, ỹ

For ν = 0, the corresponding solution ξ 1 is not constant in this case, even though the differential equation has constant coefficients, see §335D in the book by[START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF].

Recall that ϕ denotes the displacement, consistent with the fact the transformation is written as x = x + ϕ(x) and the transformation gradient as F = 1 + ∇ϕ(x). Even though this decomposition x = x + ϕ(x) is most often used in linearized or incremental elasticity, we stress that we are in the framework of finite elasticity for the moment, i.e. we retain all nonlinearities with respect to ϕ.

As seen in the figure, the numerical error on qc is larger (for fixed δλ) when p is close to p * , as the curvature of the curve of marginal stability approaches zero; accordingly, we refined the step size δλ to keep the numerical error on qc below ±.05. By running a convergence test, we have also checked that our typical mesh size, e = .02 (in units such that the width is v = 1), warrants that the error on qc due to discretization remains below this threshold.

This is because the pre-stress enter into this equation through the residual moment M 0 y only (see §4.5), which has been preserved by design.

Abdelkhalek, S., Zahrouni, H., Legrand, N., Potier-Ferry, M., 2015. Post-buckling modeling for strips under tension and residual stresses using asymptotic numerical method. Internayional Journal of Mechanical Sciences 104, 126-137.

unmodified 5 , while the quantity (1) = g ν ( (0) , δ) in (4.10b) becomes g ν ( (0) , δ) = (0) 2 1260 + 5ν+2 210 (0) -(1 -ν 2 ). Since equation (4.8) still holds, the thick green curve in the phase diagram in figure 4.4 remains unchanged. With a linear pre-strain distribution, we find that a macroscopic to microscopic transition still takes place but the critical pre-stress δ * is different: solving (1) = g ν = 0 with the expression of g ν just given yields (0) = 20.1 in the incompressible case (ν = .5), which by equation (4.8) corresponds to δ * = 63.8. This is significantly less than the value found earlier for a piecewise constant pre-strain, compare with (4.12). All the results presented in this paper can be similarly extended to an arbitrary pre-stress distribution. Thin elastic structures comprising incompatible strain are a relatively new topic which is still imperfectly understood. The thin rod model with natural curvature has been used to model these structures but the recent work of [START_REF] Cicalese | On local and global minimizers of prestrained thin elastic rods[END_REF] and our own results show that it is applicable only when the pre-strain is infinitesimal, i.e. for |p -1|

1 or δ 1 in our notation. Still, our results show that for a pre-strain as large as p * or δ * (which are both of order 1), and even slightly above, the wavelength of the buckling modes ∼ 1/q c is much larger than the cross-section dimension ∼ v ∼ h. This suggests that it must be possible to derive 1-d models governing thin structures with finite pre-strain, even though such models have not been established to date (they are certainly not of the Euler-Bernoulli type). In future work, we will address this question, using the analysis of a prismatic solid presented in §5 as a starting point.

Appendix A. Numerical shooting method for the plate model

Here, we describe a numerical solution of the eigenvalue problem (4.5a-4.5d) for the elastic plate with non-uniform pre-stress, that governs the initial buckling bifurcation. The eigenvalue problem depends on Poisson's ratio ν, on the mismatch strain δ and on the mean imposed strain ; its eigenvector is the complex deflection amplitude ξ 1 (x), and its eigenvector is the wavenumber q.

The eigenproblem can be solved by a shooting method as follows. In each half-plate j = I, II, and for a given value of q, the solutions ξ j 1 (x) that satisfy both the differential equation (4.5) and the boundary conditions (4.5b-4.5c) is of dimension 4 -2 = 2. Two particular solutions ξ j,1 1 and ξ j,2 1 spanning this space are found by integrating the differential equation with initial conditions (ξ j,k 1 ), (ξ j,k 1 ) , (ξ j,k 1 ) , (ξ j,k 1 ) x j 0 = 1, 0, νq 2 , 0 for first function, k = 1 0, 1, 0, q 2 (2 -ν) for second function, k = 2.

(A.1)

Here x I 0 and x II 0 denote the starting point of the progressive integration, which takes place in either half of the plate: x I 0 = 0 and x II 0 = 1. For each half-plate j = I, II, one can define a 4 × 2 shooting matrix S j (ν, δ, , q) by filling with the values at the interface of the successive derivatives of the numerical solution ξ j,k 1 , S j (ν, δ, , q) ik = d i ξ j,k 1 dx i (1/2). By construction, this matrix yields ξ 1 and its successive derivatives at the interface x = 1/2 in terms of the initial values ξ 1 (x j 0 ) and ξ 1 (x j 0 ) on the stress-free edge:

Finally, the continuity conditions (4.5d) can be expressed as:

det S(ν, δ, , q) = 0, (A.2a)

where the 4 × 4 shooting matrix S is obtained by assembling the 4 × 2 matrices S I and S II from either half-plate: in block-matrix notation, S(ν, δ, , q) = S I (ν, δ, , q) -S II (ν, δ, , q) . (A.2b)

Appendix B. Linearized equilibrium of the Gent model

In this section, we give the detailed expressions of the tensors appearing in the non-linear equilibrium (5.4) and in the linearized equilibrium (5.7) of the 3-d hyper-elastic model analyzed in §5.

Appendix B.1. Pre-stress and tangent moduli

We start from the strain energy per unit reference volume W 3D of a Gent material, as defined in (5.2b). The second Piola-Kirchhoff stress tensor Σ =

where C = H T • H is the Cauchy-Green deformation tensor. We rewrite the definition (5.1) of the transformation gradient associated with pre-stretch G in compact form as follows:

where r and p denote the piecewise constant principal stretches

For the special case of an invariant solution ϕ λ 0 (x), as defined in (5.5), the transformation gradients take the block-diagonal form

where the first two rows and columns correspond to the cross-sectional direction x and ỹ, and the last row and column correspond the axial direction z. The cross-sectional blocks F 0 = 1 + grad ϕ 0 and C 0 = r2 (F 0 ) T • F 0 depend on the gradients of the cross-sectional displacement ϕ 0 (λ; x, ỹ). We denote by grad the gradient with respect to the coordinates within the cross-section (x, ỹ), for instance for a scalar field v(x, ỹ): grad v = ∂v ∂ x , ∂v ∂ ỹ . As mentioned in equation (5.6), the pre-stress Σ 0 = Σ(ϕ λ 0 ) in the invariant configuration ϕ λ 0 splits up into similar blocks,

The tensor of tangent moduli appearing in (5.8) reads, when evaluated in the invariant solution,