Sébastien Neukirch

Joël Frelat

Alain Goriely

Corrado Maurini

Vibrations of post-buckled rods: the singular inextensible limit

Keywords: Vibrations, Kirchhoff elastic rods, Buckling, Bifurcation

The small-amplitude in-plane vibrations of an elastic rod clamped at both extremities are studied. The rod is modeled as an extensible, shearable, planar Kirchhoff elastic rod under large displacements and rotations, and the vibration frequencies are computed both analytically and numerically as a function of the loading. Of particular interest is the variation of mode frequencies as the load is increased through the buckling threshold. While for some modes there are no qualitative changes in the mode frequencies, other frequencies experience rapid variations after the buckling threshold, the thinner the rod, the more abrupt the variations. Eventually, a mismatch for half of the frequencies at buckling arises between the zero thickness limit of the extensible model and the inextensible model.

Introduction

The first step in the study of vibrating elastic structures [START_REF] Virgin | Vibration of axially loaded structures[END_REF] focuses on the dynamical response of the system around its unstressed configuration. In vibration analysis, the dynamics of infinitely small-amplitude disturbances around the fundamental equilibrium state are generally first considered, leading to a linear problem. However, nonlinear effects are known to play a key role in many elastic systems, and in the context of vibrations, nonlinearities can be included by studying large-amplitude oscillations around the fundamental state [START_REF] Touzé | Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes[END_REF]. The second step in the analysis of vibrations is to study the effect of external loads. They have a direct influence on the dynamical response of the system, as easily demonstrated by tuning the natural frequencies of a string by putting it under tension. Similarly, in compression, the natural frequencies of a rod decrease. Again nonlinear effects become important when external loads not only change the vibration response of the rod but also alter its overall stability through buckling. Several studies have investigated dynamical responses of post-buckled elastic rods around their post-buckled state, see e.g. Chap. 12 of [START_REF] Virgin | Vibration of axially loaded structures[END_REF] and [START_REF] Perkins | Planar vibration of an elastica arch: Theory and experiment[END_REF][START_REF] Santillan | Post-buckling and vibration of heavy beam on horizontal or inclined rigid foundation[END_REF][START_REF] Plaut | Vibration and snap-through of bent elastica strips subjected to end rotations[END_REF]. The present work also focuses on the problem of small-amplitude vibrations around a pre-strained deformed nonlinear elastic rod, and shows that vibration frequencies behave singularly at buckling. This problem is relevant for a number of applications including the manufacturing of piano (or violin) soundboards where the wooden board is bent before being clamped in the rigid metal frame [START_REF] Bilhuber | The influence of the soundboard on piano tone quality[END_REF][START_REF] Mamou-Mani | Numerical simulation of a piano soundboard under downbearing[END_REF][START_REF] Mamou-Mani | Prestressed soundboards: Analytical approach using simple systems including geometric nonlinearity[END_REF]. Other systems where prestress and/or pre-strain play an important role for the vibration response are gongs, cymbals, and steel drums where plastic deformations of the metal plates are used to introduce separation of response modes [START_REF] Chaigne | Nonlinear vibrations and chaos in gongs and cymbals[END_REF].

Here, we consider the problem of in-plane vibrations of a pre-and post-buckled Kirchhoff extensible shearable elastic rod. First, we study the equilibrium configurations of a clampedclamped rod as the ends are gradually moved together. The rod has a straight natural shape and for small axial displacements the rod remains straight until a critical displacement is reached where the rod buckles in the plane. For each value of the axial displacement we study the smallamplitude vibrations around the equilibrium state and we follow how frequencies of the natural modes change as the load is increased. For both equilibrium and vibrations we enforce displacement control boundary conditions, that is, the axial displacement rather than the axial load is imposed.

The paper is organized as follows: In Section 2 we present Kirchhoff model for elastic rods, in Section 3 we derive the equations for the small-amplitude vibrations of a rod around its postbuckled equilibrium. We compute equilibriums in Section 4, vibrations in the extensible case in Section 5 and in the extensible case in Section 6. We then compare the two cases in Section 7 and a show a discrepancy between them, which is further analyzed analytically in Section 8. Discussion (Section 9) and conclusion (Section 10) follow. We consider an elastic rod with a rectangular cross-section of width b and thickness h, total length L and arc length S in its unstressed reference state. In this state the rod lies along the e x axis, from the origin O = (0, 0, 0) to the point at (L, 0, 0). The position vector of the center of the rod cross-section is noted R(S ) and we have R(0) = (0, 0, 0) and R(L) = (L, 0, 0) in the reference state.

Model

Kinematics

We use the special Cosserat theory of rods [START_REF] Antman | Nonlinear problems of elasticity[END_REF] where the rod can suffer bending, extension, and shear deformations. We work under the assumption that the rod cross-section remains planar (and rectangular) as the rod deforms and use a set of three Cosserat directors (d 1 (S ), d 2 (S ), d 3 (S )) embedded in each cross-section: d 1 is perpendicular to the section plane, d 2 is along the small span (of length h) of the section, and d 3 is along the wide span (of length b) of the section. In the undeformed state, d 1 (S ) ≡ e x , d 2 (S ) ≡ e y , and d 3 (S ) ≡ e z . We only consider deformed states that are (i) planar (where the rod center line R(S ) lies in the (x, y) plane, the rod being bent along its small span h), and (ii) twist-less (where the director d 3 (S ) ≡ e z ). Note that in the presence of extension and shear, S may no longer be the arc-length of the curve R(S ) in the deformed state.

We introduce extension and shear strains, e 1 and e 2 , such that:

R (S ) def = dR/dS = (1 + e 1 ) d 1 + e 2 d 2 .
(1)

In the absence of extension (e 1 = 0) and shear (e 2 = 0), the director d 1 is the unit tangent to the We use the Kirchhoff dynamical equations for elastic rods [START_REF] Antman | Nonlinear problems of elasticity[END_REF], where the stresses in the section are averaged to yield an internal force N(S ) and an internal moment M(S ). These internal forces and moments are the loads exerted on the section at S by the part of the rod at S > S . In the absence of body force and couple, the linear and angular momentum balance then read

N (S , T ) = ρ h b R(S , T ), (3) M (S , T ) + R (S , T ) × N(S , T ) = ρ I θ(S , T ) , (4) where (

Constitutive law

We use the standard linear constitutive relationship relating the bending strain κ(S

) def = θ (S )
to the bending moment

M 3 def = M • d 3 : M 3 = E I κ, (5)
where E is Young's modulus. Note that κ is not the curvature in general. In a similar way, the

N 1 = E h b e 1 , (6)
N 2 = G h b e 2 . ( 7
)
where G is the shear modulus.

Equations in component form

In the planar case considered here, we have Z(S , T ) ≡ 0, N z (S , T ) ≡ 0, M x (S , T ) ≡ 0, and M y (S , T ) ≡ 0, ∀(S , T ) so that the equations for the six remaining unknowns are

X = (1 + e 1 ) cos θ -e 2 sin θ, (8a)
Y = (1 + e 1 ) sin θ + e 2 cos θ, (8b) θ = M/(EI), (8c)
M = e 2 N 1 -(1 + e 1 ) N 2 + ρI θ, (8d)
N x = ρhb Ẍ, (8e)
N y = ρhb Ÿ, (8f)
where M = M z = M 3 , N 1 = N x cos θ + N y sin θ, and N 2 = -N x sin θ + N y cos θ. The strains (e 1 , e 2 ) are given by Eqs ( 6) and ( 7) as functions of N x and N y .

Dimensionless variables

We scale all lengths with L, time with τ def = L 2 ρhb/(EI), forces with EI/L 2 , and moments with EI/L. This naturally introduces a parameter

η def = I hbL 2 = 1 12 h L 2 , (9)
which takes small values in the present case of slender rods. Dimensionless variables will be written lowercase, e.g. x def = X/L, or m def = ML/(EI). The constitutive relations ( 6) and [START_REF] Mamou-Mani | Numerical simulation of a piano soundboard under downbearing[END_REF] read:

e 1 = η n 1 and e 2 = 2(1 + ν) η n 2 ,
where the Poisson ratio ν arises from the relation

E = 2(1 + ν) G.
Rods with η > 0 are extensible, shearable rods for which the rotational inertia is accounted for. They will be simply called extensible rods. Rods with η = 0 will be simply called inextensible rods, although they really are inextensible, unshearable rods for which the rotational inertia is ignored. Such rods are frequently called elastica [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimivi propreitate gaudentes[END_REF][START_REF] Goss | The history of the planar elastica: Insights into mechanics and scientific method[END_REF].

Small-amplitude vibrations around pre-and post-buckled equilibrium

The systems of equations [START_REF] Mamou-Mani | Prestressed soundboards: Analytical approach using simple systems including geometric nonlinearity[END_REF] in dimensionless form reads

x (s, t) = cos θ + η (n 1 cos θ -2(1 + ν) n 2 sin θ) , (10a)
y (s, t) = sin θ + η (n 1 sin θ + 2(1 + ν) n 2 cos θ) , (10b) θ (s, t) = m, (10c)
m (s, t) = -n 2 + η (1 + 2ν) n 1 n 2 + θ , (10d)
n x (s, t) = ẍ, (10e)
n y (s, t) = ÿ, (10f)
with n 1 = n x cos θ + n y sin θ and n 2 = -n x sin θ + n y cos θ. In our problem, we consider a clampedclamped rod and control the end-shortening

d def = 1 -(x(1, t) -x(0, t))
. This setup implies the following boundary conditions:

x(0, t) = 0 x(1, t) = 1 -d , (11a)
y(0, t) = 0 y(1, t) = 0 , (11b) θ(0, t) = 0 θ(1, t) = 0 . (11c)
For each given value of the end-shortening d, we find the equilibrium configuration (x e , y e , θ e , m e , n xe , n ye ) by solving system [START_REF] Antman | Nonlinear problems of elasticity[END_REF] with ẍe = 0 and ÿe = 0. Then we look for small-amplitude vibrations around this equilibrium configuration, that is, we set

x(s, t) = x e (s) + δ x(s) e iωt , (12a)

y(s, t) = y e (s) + δ ȳ(s) e iωt , (12b) θ(s, t) = θ e (s) + δ θ(s) e iωt , (12c)
m(s, t) = m e (s) + δ m(s) e iωt , (12d)
n x (s, t) = n xe (s) + δ nx (s) e iωt , (12e)
n y (s, t) = n ye (s) + δ ny (s) e iωt , (12f)
where δ 1 is a small parameter, and ω is the frequency of the vibration. Inserting ( 12) into [START_REF] Antman | Nonlinear problems of elasticity[END_REF] and keeping only linear terms in δ, we obtain equations for the spatial modes ( x, ȳ, θ, m, nx , ny ):

x (s) = -θ sin θ e + η (n 1 cos θ e -2(1 + ν) n2 sin θ e ) + η θ (-n 1e sin θ e -2(1 + ν) n 2e cos θ e ) , (13a) ȳ
(s) = θ cos θ e + η (n 1 sin θ e + 2(1 + ν) n2 cos θ e ) + η θ (n 1e cos θ e -2(1 + ν) n 2e sin θ e ) , (13b) θ
(s) = m, (13c) m
(s) = -n 2 + η (1 + 2ν)(n 1 n 2e + n 1e n2 ) -ω 2 θ , (13d)
n x (s) = -ω 2 x, (13e)
n y (s) = -ω 2 ȳ, (13f)
with n1 = nx cos θ e + ny sin θ e + θ -n xe sin θ e + n ye cos θ e and n2 = -n x sin θ e + ny cos θ e + θ -n xe cos θ en ye sin θ e . The boundary conditions on the spatial modes are

x(0) = 0 x(1) = 0 , (14a) ȳ
(0) = 0 ȳ(1) = 0 , (14b) θ(0) = 0 θ(1) = 0 . (14c)
For given parameters η and ν and given end-shortening d, the equilibrium (x e , y e , θ e , m e , n xe , n ye ) is first computed from [START_REF] Antman | Nonlinear problems of elasticity[END_REF] with ẍe = 0 and ÿe = 0. Then the 6D system [START_REF] Yamaki | Non-linear vibrations of a clamped beam with initial deflection and initial axial displacement, part I: Theory[END_REF] with the six boundary conditions ( 14) is a well-defined boundary value problem, but with the additional unknown ω. For computational purpose, we normalize the linear solution of this problem by imposing the condition

m2 (0) + n2 x (0) + n2 y (0) = 1. ( 15
)

Equilibrium

We use a 'home-made' predictor-corrector path-following code to address the problem numerically. For each value of p In Fig. 2 equilibrium paths are given for both inextensible (i.e. η = 0) and extensible (η > 0)

rods. In the latter case, we see that the displacement d starts to increase as soon the curve leaves the origin: due to extensibility the rod shorten before buckling. Buckling happens at p = 4π 2 for the inextensible case and at lower values for extensible cases. These equilibrium paths show that extension and shear play a minor role in the buckling load and that the inextensible rod solution is obtained in the limit η → 0 of the extensible rod solution, as expected. Closed-form formula for the dotted (red) curve is given in Appendix A.

Vibrations in the extensible case

Vibrations around the straight state

We first consider the equilibrium solution where the axially loaded rod is straight, but compressed:

y e (s) = 0, x e (s) = (1 -ηp) s, n ye (s) = 0, n xe (s) = -p, θ e (s) = 0, m e (s) = 0. ( 16
)
The vibrations around this straight equilibrium are either extensional of flexural. Extensional vibrations modes are solutions of:

x (s) = η nx (17a) n x (s) = -ω 2 x (17b)
with boundary conditions x(0) = 0 = x(1). This yields vibrations frequencies of the form ω = j π/ √ η (with j = 1, 2, . . .), which do not vary with the load p. The upper curve, before buckling of Fig. 3-(a) is such a solution. Flexural vibrations modes are solutions of:

ȳ + µ p ȳ -µ 2 ω ȳ = 0, (18)
with boundary conditions ȳ(0) = 0 = ȳ(1) and

η 2 ω 2 + [1 -p(η -η )] 2 ȳ (0, 1) + η ȳ (0, 1) = 0. We introduced the notations η = 2η(1 + ν), µ p = (1 -p(η -η ))p + ω 2 (η + η ), and
µ 2 ω = ω 2 [1 -p(η -η )] 2 -[1 -p(η -η )]pη -η η ω 2 .
The general solution is of the form:

ȳ(s) = A cos k + s + B cosh k -s + C sin k + s + D sinh k -s (19)
where we used the two wave numbers

k ± = (1/ √ 2) µ 2 p + 4µ 2 ω ± µ p . The boundary condi- tions require that 2 a + a -k + k -cos k + cosh k --1 + (a + k + ) 2 -(a -k -) 2 sin k + sinh k -= 0 , (20)
with a± = [1 -p(η -η )] 2 + η ∓(k ± ) 2 + η ω 2 ,
which is an equation for ω. The analytical solutions ω(p) given by (20) match numerical solutions given in Fig. b,c). Buckling occurs when ω = 0 = µ ω , i.e. for k + = √ µ p = 2π and k -= 0. This yields:

p = 1 + 16π 2 (η -η) -1 2(η -η) = 4π 2 -16π 4 (η -η) + O(η 2 ) (21)

Vibrations around the buckled state

Once the equilibrium solution is known, we solve the boundary value problem ( 13)-( 14)

numerically with a shooting method: (i) we first use a guess for the unknown parameters χ = ( m(0), nx (0), ny (0), ω) and we integrate the system (13) up to s = 1; (ii) we then check if the boundary conditions ( 14)-( 15) are satisfied. If not, we change the guess χ accordingly (using a Newton-Raphson scheme) until the boundary conditions at s = 1 are satisfied. In Fig. 3, frequencies for the first four modes are given as a function of the parameter p, for different values of η. It should be noted that the computations performed here are for a displacement control experiment, that is d (and not p) is controlled. Nevertheless, as there is a one-to-one correspondence between d and p, for each value of the applied longitudinal displacement d, the equilibrium axial load p is read from Fig. 2 and then the frequency is computed and plotted in Fig. 3. We see in Fig. 3 that as η → 0, frequencies globally increase and tend toward limiting curves. Finally we note that every curve ω(p) is continuous, but that the curves for the odd modes experience a rapid increase just after buckling. As η → 0 this rapid increase becomes more abrupt to eventually turn into a discontinuity in the case η = 0, see Section 6. In each of We now focus on the first mode, which emerges from ω = 0 at buckling. In Fig. 7-(a), we plot ω as function of the rise of the buckled rod at its mid-point: y e (1/2) = Y e (L/2)/L. Each curve corresponds to a different value of the parameter η, from η = 1/1200 (i.e. L = 10 h) to η = 1/480 000 (i.e. L = 200 h). We see that all curves emerge from ω = 0 at buckling (y e (1/2) = 0) and asymptotically tend to the curve computed in the inextensible case (see Section 6) when y e (1/2) becomes large. For very small η values, curves rise sharply from ω = 0 and quickly approaches the inextensible asymptote. As a matter of fact these curves can be made to almost collapse on a master curve if the horizontal axis is plotted in unit of the rod thickness h: in Fig. 7-(b) we plot the frequency ω as function of Y e (L/2)/h, for the same set of η values. All curves nearly collapse on a master curve which has a (numerically determined) slope 28 at the origin. Using h = 12η L we obtain:

ω 8.1 √ η Y e (1/2) L , (22)
which means that the limit η → 0 is singular for this first mode. The same phenomenon happens for all odd modes.

In dimensional form (i.e. Ω = ω EI/(ρhb) /L 2 ) we have the following expression for the frequency Ω (in rad/s):

Ω 8.1 Y e (L/2) L 2 E ρ , for Y e (L/2) 2h . ( 23
)
The presence of the celerity c = E/ρ of compression elastic waves in this expression shows that directly after buckling and for a short loading interval (i.e. Y e (L/2) = 0 to Y e (L/2) 2h), the lowest mode of vibration of a buckled rod is of extension-compression type. We also see in The curves nearly collapse on a master curve whose slope at the origin is 28. The dashed curve (red online), whose slope at the origin is 2 √ 2 π 2 27.9 and which is hardly distinguishable from the previous ones, is the first mode solution of (25). (c) Same data as in (a) but with d = D/L on the horizontal axis. (d) Same data but with the horizontal axis rescaled with (4π 2 η); As in (b) the curves nearly collapse on a master curve. The dashed curve (red online), which is hardly distinguishable from the previous ones, is the first mode solution of (25).

of the rod could be a way to define the notion of a shallow (resp. deep) buckled equilibrium shape: a shallow (resp. deep) equilibrium shape has a vibrational response that is primarily extensional (resp. flexural).

We now focus on the second mode and show that it does not suffer the same singularity as the first one. The mode frequency emerges from a finite value ω b = ω b (η) at buckling, and we see in Fig. 3 that its variation (with p) after buckling is much slower than for the first mode. In Figure 8 the variation of the second mode frequency is plotted as a function of y e (1/2) = Y e (L/2)/L. We numerically extract the approximation:

ω ω b (η) -(9 -5400 η) Y e (L/2) L 2 (24)
In conclusion, the limit η → 0 is smooth for the second mode (and in fact for all even modes).

An alternative equation is used in (for example) [START_REF] Woinowsky-Krieger | The effect of an axial force on the vibration of hinged bars[END_REF][START_REF] Yamaki | Non-linear vibrations of a clamped beam with initial deflection and initial axial displacement, part I: Theory[END_REF][START_REF] Yamaki | Non-linear vibrations of a clamped beam with initial deflection and initial axial displacement, part II: Experiment[END_REF][START_REF] Nayfeh | Exact solution and stability of postbuckling configurations of beams[END_REF] for the vibrations of postbuckled extensible rods. The transverse displacement Y(S , T ) is there solution of the equation: with

EI Y + ρ h b Ÿ + PY = 0 (25a) 10
P = Ehb L D - 1 2 L 0 Y 2 dS . (25b)
In this model a certain number of assumptions on the nonlinear terms are made, see e.g. Eq. (B.3).

Calculating the frequency of the first mode just after buckling, we find:

Ω 2/3 π 2 Y e (L/2) L 2 E ρ (26)
where √ 2/3 π 2 8.06, in agreement with our result (23). We plot in Fig. 7-(b) and (d) the first mode solution of Eq. ( 25), and we note that this solution is independent of η, which is not the case for Kirchhoff equations. Moreover the second mode of (25) is found to be independent of Y e (L/2) (see Fig. 8), in contradiction with (24). A more comprehensive comparison of the two models is the subject of a forthcoming paper.

Vibrations in the inextensible case

Vibrations around the straight state

We consider an inextensible unshearable rod of length L, strongly held at both sides by clamps separated by a distance L. The rod is held straight and hence no flexural dynamics at all can take place, i.e. flexural vibrations (that would be given by the η = 0 version of Eq. ( 18))

are impossible here as they would require shortening of the ends.

Axial vibrations are given by:

x (s) = 0 (27a) n x (s) = -ω 2 x (27b)
with boundary conditions x(0) = 0 = x(1). The solution is x(s) ≡ 0, and nx (s) ≡ const., where we see that any ω is admissible. Even if no extensional deformation is present, the load in this statically indeterminate problem can fluctuate with any frequency. This is illustrated in Fig. 9 by the (red) shaded region in the interval 0 ≤ p ≤ 4π 2 .

Vibrations around the buckled state

As in the extensible case we solve the boundary value problem ( 13)-( 14) numerically with a shooting method. The results are shown in Fig. 9 and in particular we see that at the buckling threshold, the first four frequencies are: ω 44.36, 74.

Comparison of the extensible and inextensible results

We here compare the post-buckled vibrations frequencies obtained numerically for the extensible (η > 0) and inextensible (η = 0) cases. On the one hand, we see in Fig. 9 that in the inextensible case there is no curve emerging from ω = 0 at buckling, whereas there always is one for all η values in the extensible case (Fig. 3). On the other hand, in the inextensible case there is a curve emerging from ω 74.4, nowhere near any of the extensible curves. We plot in Fig. 12-(a) a comparison of the inextensible η = 0 and extensible η = 1/19 200 cases, for the first two modes. We see that the second extensible mode (emerging from ω 44 at buckling) is always very close to its inextensible counterpart, but that the first extensible mode (emerging from ω = 0 at buckling) is first very far from its inextensible counterpart and only approaches it later in the post-buckling regime. The consequence is that for a short interval after buckling the inextensible model wrongly predicts a fundamental mode at ω 44 whereas for any 'real life' elastic rod there is a vibration mode with lower frequency.

A full comparison between extensible and inextensible cases shows that all extensible odd modes experiences a similar mismatch with their inextensible counterpart, see Fig. 12-(b) for the first four modes. For the solution at buckling, with p given by (21), we further investigate numerically the variation with η of the frequencies of the first four modes in the extensible case and compare them with the first four frequency of the inextensible case η = 0. In Fig. 13, we see that as η → 0 the even extensible modes are converging to solutions of the inextensible case, but that the limit η → 0 of the odd extensible modes does not correspond to the frequency values found in the inextensible case.

Analytical study for the inextensible case

In order to understand the mismatch of the odd modes' frequencies at buckling between the extensible and inextensible cases, we look at the problem analytically. In the inextensible case, buckling happens for p = 4π 2 . We look for the frequencies that emerge from the continuum present for p ≤ 4π 2 . If the η → 0 limit was not singular we would just look for the solutions of (20) for η = 0 and p = 4π 2 , that is:

2 k + k -cos k + cosh k --1 + 4π 2 sin k + sinh k - (28)
with

k ± = √ 4π 4 + ω 2 ± 2π 2 .
The solutions are listed in Table 1 and we see that, for the odd modes, they do not match what is found numerically.

To investigate matter further we set η = 0 in Eqs. [START_REF] Antman | Nonlinear problems of elasticity[END_REF], and look for the frequency values that exist just after buckling (i.e. around post-buckled configurations). In the post-buckled configuration, the equilibrium equations and the boundary conditions are given by n ye = const. (29a) 20) for p = 4π 2 and η = 0. This also corresponds to the first eight solutions of P 1 (ω 0 ) = 0 (see Eq. ( 38)).

n xe = -p (
and without loss of generality, we choose x e (0) = y e (0) = 0. For equilibrium modes whose shapes are invariant when reflected along the line parallel to the e y axis and containing the point (x(1/2), 0), we have n ye = const. = 0 in Eq. (29a). The first bifurcated equilibrium mode, represented in Fig. 1 and on which we focus, is such a mode. We address the behavior of the solutions after but close to buckling. Therefore, we expand the variables θ e (s) and y e (s) in powers of , a small parameter measuring the distance from buckling:

θ e (s) = θ 1 (s) + 2 θ 2 (s) + 3 θ 3 (s) + O 4 (30a) x e (s) = x 1 (s) + 2 x 2 (s) + 3 x 3 (s) + O 4 (30b) y e (s) = y 1 (s) + 2 y 2 (s) + 3 y 3 (s) + O 4 (30c) p = p 0 + p 1 + 2 p 2 + 3 p 3 + O 4 (30d)
We substitute these expansions in the equilibrium equations (29), which have to be satisfied to all orders in . The solution up to order 3 reads: θ e (s) = sin 2πs +

p = 4π 2 + 2 π 2 /2 + O 4 (31d)
In order to relate to the control parameter d and the amplitude after bifurcation, we compute the end-shortening

d = 1 -(x e (1) -x e (0)) = 2 /4 + O 4 = 2 p 4π 2 -1 + O 4 , (32)
and the rod maximum deflection

y e (1/2) = π 1 - 5 48 2 + O 4 . ( 33
)
Vibration around the post-buckled equilibrium

We expand all modal variables ( x, ȳ, θ, m, nx , ny ) and the frequency ω in powers of . For instance, we have ω = ω 0 + ω 1 + 2 ω 2 + 3 ω 3 + O 4 , and so on. We can now solve equations [START_REF] Yamaki | Non-linear vibrations of a clamped beam with initial deflection and initial axial displacement, part I: Theory[END_REF] (with η = 0) with boundary conditions [START_REF] Yamaki | Non-linear vibrations of a clamped beam with initial deflection and initial axial displacement, part II: Experiment[END_REF], using the equilibrium solution (31). To order 0 we have:

x 0 = 0 with x0 (0) = 0 = x0 (1) (34a)
n x0 = -ω 2 0 x0 (34b)
ȳ 0 + 4π 2 ȳ 0 -ω 2 0 ȳ0 = 0 with ȳ0 (0) = ȳ0 (1) = ȳ 0 (0) = ȳ 0 (1) = 0 (34c)

The first two equations describe the longitudinal mode and are decoupled from the third one which is associated with the transverse mode. More precisely, the longitudinal mode is given by x0 (s) = 0 and nx0 (s) constant .

whereas for the transverse mode, the solution ȳ0 (s) is

ȳ0 (s) = A 0 cos k + 0 s -cosh k - 0 s cos k + 0 -cosh k - 0 - k - 0 sin k + 0 s -k + 0 sinh k - 0 s k - 0 sin k + 0 -k + 0 sinh k - 0 , (36) with k
± 0 def = 4π 4 + ω 2 0 ± 2π 2 .
The boundary conditions impose that

A 0 P 1 (ω 0 ) = 0 with (37) P 1 (ω 0 ) def = 2 k + 0 k - 0 cos k + 0 cosh k - 0 -1 + 4π 2 sin k + 0 sinh k - 0 (38)
which is an equation for ω 0 . The first eight solutions of P 1 (ω 0 ) = 0 are listed in Table 1. Note that to order 0 there is no frequency jump. To order 1 we have

x 1 = θ0 sin 2πs with x1 (0) = 0 = x1 (1) (39a)
n x1 = -ω 2 0 x1 -2 ω 0 ω 1 x0 (39b) ȳ 1 + 4π 2 ȳ 1 -ω 2 0 ȳ1 = 2 ω 0 ω 1 ȳ0 + 2π nx0 cos 2πs (39c) with ȳ1 (0) = ȳ1 (1) = ȳ 1 (0) = ȳ 1 (1) = 0 . (39d)
We start by solving equation (39a). The boundary condition x1 (1) = 0 implies

A 0 P 2 (ω 0 ) = 0 with (40) P 2 (ω 0 ) def = k + 0 k - 0 (k + 0 2 -2π 2 )(cos k + 0 -cosh k - 0 ) + 2π 2 (cos k + 0 cosh k - 0 -1) +(k + 0 2 k - 0 2 + 8π 4 ) sin k + 0 sinh k - 0 . (41)
The solutions have to satisfy Eqs. (37) and (40), which are transcendental equations for ω 0 . A numerical root finding analysis reveals that P 1 (ω 0 ) = 0 and P 2 (ω 0 ) = 0 share half of their roots, see Table 2 where columns with an even index correspond to common roots and match numerical values at p = 4π 2 for the continuous curves plotted in Fig. 3-(d). In the case of a common root, Eqs. (37) and (40) are fulfilled for non-vanishing A 0 , and the corresponding modes have frequencies that are continuous in the control parameter close to buckling.

In the case of distinct roots, we are compelled to set A 0 = 0. In this case where P 1 (ω 0 ) 0, we solve Eqs. (39a) and (39b) to obtain x1 (s) = 0 and nx1 (s) constant. The general solution of Eqs. (39c), (39d) is then where

ȳ1 (s) = 2π nx0 k + 0 2 k - 0 2 P 1 (ω 0 ) c 1 k + 0 cos k + 0 s -c 1 k + 0 cosh k - 0 s -c 2 k - 0 sin k + 0 s + c 2 k + 0 sinh k - 0 s +P 1 (ω 0 )(cos k + 0 s -cos 2πs) (
c 1 = k - 0 (cos k + 0 -1)(cosh k - 0 + 1) + k + 0 sin k + 0 sinh k - 0 ( 43
)
c 2 = k + 0 sin k + 0 (cosh k - 0 -1) + k - 0 (cos k + 0 -1) sinh k - 0 (44)
In order to select a mode, we need to proceed to order 2 , which reads

x 2 (s) = -ȳ 1 (s) sin 2πs with x2 (0) = 0 = x2 (1) (45)
The boundary conditions at s = 1 impose P 3 (ω 0 )/P 1 (ω 0 ) = 0 where

P 3 (ω 0 ) = 2k - 0 2 k + 0 2 (cosh k - 0 cos k + 0 -1) -4k + 0 (k - 0 2 + k + 0 2 ) sin k + 0 (cosh k - 0 -1) +(4k - 0 (k - 0 2 + k + 0 2 )(1 -cos k + 0 ) + 4π 2 k - 0 k + 0 sin k + 0 ) sinh k - 0 (46)
Oddly enough this function P 3 (ω 0 ) has the same set of common roots as P 1 (ω 0 ) and P 2 (ω 0 ), see Table 3. The other roots correspond to the frequency values at which the (discontinuous) odd mode curves emerge from p = 4π 2 in Fig. 3-(d which implies that the common roots must verify:

k + 0 2 tanh k - 0 2 = k - 0 2 tan k + 0 2 (48)
These roots correspond to frequencies that do not vary abruptly after buckling has occurred. An approximate formula is k + 0 π/2 + 2 jπ (with positive integers j), which yields ω 0 (π/2) 2 (4 j + 1) (4 j -3)(4 j + 5) (corresponding to columns with i > 1 even in Table 3).

Formulas for the roots of the three functions P 1 , P 2 , and P 3 in the limit of large k + 0 are given in Appendix C. In particular it is shown that the set of roots of P 3 which is not in common with P 1 and P 2 is such that k + 0 3π/2+2 jπ (with positive integers j). This implies that the frequencies emerging from buckling are such that ω 0 (π/2) 2 (4 j + 3) (4 j + 7)(4 j -1) (corresponding to columns with i > 1 odd in Table 3).

Discussion

In [START_REF] Marigo | Hierarchy of one-dimensional models in nonlinear elasticity[END_REF] it is shown that, under certain hypotheses, in the zero thickness limit a slender cylinder behaves either as a flexible rod or as an extensible string and that one does not need to consider single bump wave to a double bump wave. Anecdotically, we note that piano soundboards are precisely tuned in this parameter range (e.g. h = 1 cm, L = 2 m, and Y e (L/2) = 1 cm). To which extent this toy model is relevant for the real problem of the piano soundboard and whether piano manufacturers are using the distinction between shallow shapes (where extension prevails) and deep buckled shapes (where bending prevails) to enrich the sound is a tantalizing idea that deserves further attention. We also note that in the case of the piano, the soundboard is a plate and not a rod, and that furthermore the soundboard is linked to the piano strings by a bridge, so really the whole system should be studied.

Conclusion

We have studied the in-plane vibrations of a slender extensible and shearable elastic rod around its post-buckled equilibrium configuration, in the displacement control case. We have shown that after buckling there is a narrow window in the loading parameter values in which half of the vibration frequencies vary abruptly. As the thickness of the rod is decreased, the vibration frequencies tend toward limiting values which do not correspond to what is found with a fully inextensible, unshearable model (i.e. an elastica). This mismatch present for half of the vibration frequencies has been found numerically and proved analytically. We conclude that we have identified a loading setup where the elastica model for slender rods fails to give the correct answer.

Figure 1 :

 1 Figure 1: Clamped-clamped rod buckled in the (x, y) plane. The end-shortening D is controlled. The point A in the reference configuration moves to point a in the deformed configuration, introducing horizontal U ≤ 0 and vertical V displacements. The origin is taken at the fixed point O at the left end of the rod.

  ∂T , T is time, ρ the mass per unit volume of the material, and I the second moment of area of the cross-section (in the present case I = h 3 b/12).

tension N 1 def= N • d 1 and the shear force N 2 def= N • d 2

 122 are linked to the extension e 1 and shear strains e 2 through

def=Figure 2 :

 2 Figure 2: Fundamental and post-buckled equilibrium path of a clamped-clamped rod with η = 1/1200, η = 1/4800, and η = 0 (top to bottom). As the controlled end-shortening d = D/L is gradually raised, an increasing axial load p = -N xe L 2 /(EI) is recorded.

Figure 3 :

 3 Figure 3: Frequencies for the lowest four vibration modes of a clamped-clamped rod around its fundamental and postbuckled equilibrium configurations in the extensible case: (a) for η = 1/1200, (b) for η = 1/4800, (c) for η = 1/19200 (plain). The labels A i , B i , C i with i = 1, 2, 3, 4 correspond to the shapes given in Figs. 4, 5, 6.Once a solution χ i is found for a given p = p i , we set p = p i+1 and use the value χ i as starting guess (predictor step) for the shooting method at p = p i+1 (corrector step). In this setup, each curve ω = ω(p) represents a path in the numerical bifurcation diagram, and we have numerically computed the four first paths (i.e. lowest four curves ω = ω(p)) for several values of the parameter η = I/(hwL 2 ).

Figs. 4 , 5 ,Figure 4 :

 454 Figs.[START_REF] Santillan | Post-buckling and vibration of heavy beam on horizontal or inclined rigid foundation[END_REF], 5, and 6, dynamical shapes (x(s, t) = x e (s) + x(s) cos ωt, y(s, t) = y e (s) + ȳ(s) cos ωt) of the vibrating rod are plotted.

Figure 5 : 4 Figure 6 :

 546 Figure 5: First four modes at d = 0.05, for η = 1/4800: configurations D 1 (a), D 2 (b), D 3 (c), and D 4 (d) in Figure 3.

Fig. 7 -Figure 7 :

 77 Fig. 7-(b) that for Y e (L/2) 8h the behavior is of flexural type (i.e. curves have reached the inextensible asymptote). This separation between two different behaviors in the elastic response

Figure 8 :

 8 Figure 8: Post-buckled frequency of the second mode (for η > 0). Continuous (black)) curves are, from bottom to top, for 1/η = 1200, 4800, 19200. The (red) dashed curve is the second mode solution of (25).

3 'Figure 9 :

 39 Figure9: Frequencies for the lowest four vibration modes of a clamped-clamped rod around its fundamental and postbuckled equilibrium configurations in the inextensible case η = 0. The labels B i , C i with i = 1, 2, 3, 4 correspond to the shapes given in Figs.10, 11.

4 'Figure 10 : 3 . 4 'Figure 11 :

 4103411 Figure 10: First four modes at d = 0.05, in the inextensional case: configurations B 1 (a), B 2 (b), B 3 (c), and B 4 (d) in Figure 3.

Figure 12 :

 12 Figure 12: Frequencies of a clamped-clamped rod around its post-buckled equilibrium configurations. Comparison of the extensible case (η = 1/19 200, plain curve) with the inextensible case (dotted, red) (a) for lowest two vibration modes, and (b) for lowest four vibration modes and plotted as a function of Y e (L/2), the rise of the buckled equilibrium shape at its middle point.

Figure 13 :

 13 Figure 13: The curves show the variation with η of the frequencies of the first four modes in the extensible case, at the buckling threshold. The (red) dots at η = 0 show the frequency values of the first four modes of the inextensible case, at the buckling threshold. A mismatch for half of the modes exists.

  29b) θ e = -p sin θ en ye cos θ e with θ e (0) = 0 = θ e (1) (29c) x e = cos θ e with x e (1)x e (0) = 1d (29d) y e = sin θ e with y e (0) = y e (1

  + 23 cos(2πs) -3 cos(6πs)) + O 4 (31c)

  ). Moreover one can verify thatP 3 (ω 0 ) = k + 0 k - 0 P 1 (ω 0 ) -

Figure D. 14 :

 14 Figure D.14: Frequencies for the lowest four vibration modes of a clamped-clamped rod around its fundamental and post-buckled equilibrium configurations, in load-control boundary conditions. Plain curves are for the extensible case with η = 1/4800, and dotted (red) curves are for the inextensible case (η = 0).

Table 1 :

 1 Frequencies and wave numbers for the lowest eight modes of vibration as given by Eq. (

	),	(29e)
	14

Table 2 :

 2 First eight solutions of P 2 (ω 0 ) = 0 (see Eq. (41)).

	42)

Acknowledgements

It is a pleasure to thank Olivier Thomas for discussions on Eq. ( 25). This publication is based in part upon work supported by Award No. KUK-C1-013-04, made by King Abdullah University of Science and Technology (KAUST) (AG). AG is a Wolfson/Royal Society Merit Award holder.

both bending and extension in a slender rod problem. Marigo and Meunier [START_REF] Marigo | Hierarchy of one-dimensional models in nonlinear elasticity[END_REF] result was obtain for the statics of a clamped-free rod subjected to body and concentrated forces. We have showed here that for displacement control clamped-clamped boundary conditions a model including both bending and extension is necessary to correctly describe the vibrations of the rod. One could then ask which constraint should be relaxed in order to fall back into the results of [START_REF] Marigo | Hierarchy of one-dimensional models in nonlinear elasticity[END_REF]. We simply remark that in the load control case, see Appendix D, the inextensible limit is smooth and hence a model including bending alone correctly captures the vibrational behavior of the rod.

From a theoretical point of view, the singular inextensible limit in the displacement control case is rather surprising and unexpected. Indeed, Equations ( 13) can be easily recast in the form of a classical eigenvalue problem

where X is the six-dimensional vector build from the six normal mode variables, L is a firstorder linear operator in L 2 ([0, 1]), (the set of square integrable functions on the unit interval)

and M is an inertia matrix. Therefore, a naive application of the classical theory of perturbation of eigenvalues for linear operators would suggest that once the eigenvalues have been found

for a value of the parameter η, they can be locally continued in this parameter. That is we would expect the curves ω(p) of Fig. 3 to vary continuously as the parameter η is decreased to zero. While it is true for some curves, other curves behaves singularly in the η → 0 limit. The fundamental mathematical reason for this phenomenon is that both the null spaces of the inertia matrix M and the linear operator L have a non-empty intersection in the inextensible (η = 0) case. Mathematically, the classical theory does not apply and new conditions for the analytic continuation of frequencies with respect to the parameters emerge. Some frequencies satisfy these relations (and hence can be analytically continued) while others do not. Finally, also due to the joint degeneracy of M and L, there exists for η = 0 a continuous family of solutions at the critical load, given by an arbitrary increase along nx (i.e. solutions of Eq. ( 27)). For larger values of the load but still at η = 0, some of these solutions are selected and emerge, apparently out of the blue. The mathematical structure of these linear problems and how they are related to various limits (nearly inextensible rods) deserves further attention.

A tentative classification for the modes is given by the number of nodes present in 0 < s < 1, a mode i having i -1 nodes and i antinodes. We see that this classification works for the modes in the pre-buckling regime (Fig. 4), but fails for the third mode in the inextensible (η = 0) case (see shape B 3 in Fig. 10 and C 3 in Fig. 11), as well as in the extensible rod (η > 0) case provided we are far enough in the post-buckling regime (see shape C 3 in Fig. 6).

An interesting situation arises when the first mode intersects with the second mode. This can easily be obtained by computing an approximation of the buckled rod rise Y e (L/2) for which the first mode frequency ω (1) 28 Y e (L/2)/h meets the second mode frequency ω (2) 44. This happens for Y e (L/2) 1.6h. At this height the shape of the fundamental mode changes from a Appendix A. Closed form solution for the planar elastica

In the inextensible, unshearable case, the solution of system [START_REF] Antman | Nonlinear problems of elasticity[END_REF] with ẍe = 0 and ÿe = 0 corresponds to the equilibrium of the planar elastica, first studied by Euler [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimivi propreitate gaudentes[END_REF] (see [START_REF] Goss | The history of the planar elastica: Insights into mechanics and scientific method[END_REF] for an historical account).

Closed-form solutions of these equations can be written in terms of elliptic functions, see e.g. [START_REF] Nizette | Toward a classification of Euler-Kirchhoff filaments[END_REF]. In particular the dotted (red) curve in Fig. 2 has the parametric expression:

with λ ∈ [0, 1). The elliptic integrals are defined as

dθ and E(λ) = notations one uses the displacements (see Fig. 1):

also functions of the arc-length S of the reference configuration. In the current configuration the derivative of the current position with regard to S will not yield a unit vector if extension occurs.

If we restrict to the case where shear is not present (i.e. e 2 = 0), we have X (S ) 2 + Y (S ) 2 =

(1 + e 1 ) 2 which yields:

Usually the right-hand side is neglected and the following von Karman approximation is used:

Appendix C. Analytical formulas for the roots of the functions P 1 , P 2 , and P 3

In the limit of large k + 0 , that is in the limit of large k - 0 = k + 0 2 -4π 2 and high frequencies

we have:

Function P 1 (ω 0 ) (see Eq. ( 38)) then reads:

and the solutions to P 1 (ω 0 ) = 0, for large k + 0 , are:

where j is a large integer, and i is the mode number (i.e. column number in Table 1). In the same limit of high frequencies, function P 2 (ω 0 ) (see Eq. ( 41)) reads:

and the solutions to P 2 (ω 0 ) = 0, for large k + 0 , are:

where j is a large integer, and i is the mode number (i.e. column number in Table 2). In the same limit of high frequencies, function P 3 (ω 0 ) (see Eq. ( 46)) reads:

and the solutions to P 3 (ω 0 ) = 0, for large k + 0 , are:

where j is a large integer, and i is the mode number (i.e. column number in Table 3). We see that in this limit the three functions P 1 , P 2 , and P 3 share half of their roots, namely those given by (C.3), (C.7), or (C.9).

Appendix D. Smooth inextensible limit in the load control case

We here show that in the case where the axial load p is controlled, the inextensible limit (η → 0) is no longer singular. The boundary conditions for this case read:

x(0, t) = 0 n x (1, t) = -p , (D.1a) y(0, t) = 0 y(1, t) = 0 , (D.1b) θ(0, t) = 0 θ(1, t) = 0 . (D.1c)

Vibrations of an inextensible rod are then possible before buckling as the right end is allowed to move axially. We see in Fig. D.14 that the inextensible case η = 0 is obtain as the smooth limit η → 0 of extensible case, and that no mismatch is present (in particular there is a mode starting from ω = 0 after buckling in the inextensible case).