Marwan Hariz

Loïc Le Marrec

Jean Lerbet

Explicit analysis of large transformation of a Timoshenko beam : Post-buckling solution, bifurcation and catastrophes

Keywords: Timoshenko beam, Large transformation, Boundary/Cauchy value problem, Dead and follower loads, Quasi static perturbation, Bifurcation and catastrophe

Explicit analysis of large transformation of a Timoshenko beam : Post-buckling solution, bifurcation and catastrophes

Introduction

The study of deformation of beams has started in the 18th century with Euler and Bernoulli's classical beam theory [START_REF] Han | Dynamics of transversely vibrating beams using four engineering theories[END_REF]. This standard theory provides reasonable approximations for many problems especially for slender structure. This model uses kinematical hypotheses: cross-section remaining normal to the center-line and sometimes inextensibility of the center-line. Large transformation of Euler-Bernoulli beam was already addressed in Euler's work at 1744. However such kinematical hypotheses are not easily justified under such transformations. Additionally Lagrange multipliers may increase the difficulty to formulate the non-linear problem [START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with lagrange multipliers and a perturbation solution[END_REF]. Later on, Timoshenko relaxes the kinematical hypotheses introducing stress-strain relations for all degrees of freedom [START_REF] Timoshenko | on the transverse vibrations of bars of uniform cross-section[END_REF][START_REF] Timoshenko | Theory of elastic stability[END_REF]. In the linear case the Timoshenko model allows the consideration of shorter beam and less-standard material [START_REF] Elishakoff | Handbook on Timoshenko-Ehrenfest beam and Uflyand-Mindlin plate theories[END_REF]. Moreover it provides an adapted formulation for large transformation where non-linear coupling between several strains intervenes. Among others, Mohyeddin and Fereidoon [START_REF] Mohyeddin | An analytical solution for the large deflection problem of timoshenko beams under three-point bending[END_REF] offered analytical solutions for the large deflection of a simply supported Timoshenko beam. Li investigated a closed problem both analytically and numerically [START_REF] Li | Large deflection and rotation of timoshenko beams with frictional end supports under three-point bending[END_REF]. Even if these works offer analytical expressions, such formulation cannot be extended easily to more general problems as their formulations are not covariant. Indeed a geometrical-material language is more adapted.

In such point of view the Timoshenko beam can be seen as a one-dimensional Cosserat body [START_REF] Cosserat | Théorie des corps déformables[END_REF]. Following this approach it is possible to extend the Timoshenko initially linearized theory to large transformations [START_REF] Eugster | A variational formulation of classical nonlinear beam theories[END_REF]. Accordingly Antman developed non-linear theory of elastic bodies such as strings, rods, beams [START_REF] Antman | Nonlinear problems of elasticity[END_REF]. Reissner adopted the same formulation to give the principle of virtual work leading to constitutive and equilibrium equations [START_REF] Reissner | On one-dimensional finite-strain beam theory: the plane problem[END_REF]. Following the same approach, Simo examined a numerical formulation of a finite strain beam theory and gave a stress-strain relation [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. i[END_REF]. Later on, these results where expanded by Rakotomanana who regarded waves and vibrations of strings, beams and shells [START_REF] Rakotomanana | Eléments de dynamique des solides et structures déformables[END_REF]. Henceforth, Le Marrec et al. gave an exact theory of Timoshenko beam undergoing three-dimensional finite transformation and subjected to dynamical perturbations [START_REF] Marrec | Vibration of a timoshenko beam supporting arbitrary large pre-deformation[END_REF] . Stability of the equilibrium solutions is of crucial importance under such large transformations (e.g. [START_REF] Bigoni | Nonlinear solid mechanics: bifurcation theory and material instability[END_REF]).

Buckling of straight beam was extended to ring by Reissner [START_REF] Reissner | Some remarks on the problem of column buckling[END_REF]. Bažant and Cedolic investigated stability of elastic structures using energy methods [START_REF] Cedolin | Stability of structures: elastic, inelastic, fracture and damage theories[END_REF]. The link with incremental equilibrium equations was investigated in [START_REF] Bažant | A correlation study of formulations of incremental deformation and stability of continuous bodies[END_REF] where an example of flexure and shear of a Timoshenko beam is presented. However energetic approach is not the unique approach. For example the numerical solutions of the extensible Timoshenko beam model under distributed load performed in [START_REF] Dell'isola | Extensible beam models in large deformation under distributed loading: A numerical study on multiplicity of solutions[END_REF] motivates Corte et al. to identify sequences of equilibria among which two at most are stable [START_REF] Della Corte | Large deformations of timoshenko and euler beams under distributed load[END_REF]. In the previous mentioned works analytical expression of a solution is not the main objective. In practice such analytical solutions mainly invoque Jacobian elliptic functions. These functions appear as solutions of many important problems in classical mechanics. Mathematical background can be found in [START_REF] Meyer | Jacobi elliptic functions from a dynamical systems point of view[END_REF][START_REF] Baker | Jacobi elliptic functions and the complete solution to the bead on the hoop problem[END_REF] (among others) and fundamental relations are reported in [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF]. Ohtsuki [START_REF] Ohtsuki | An analysis of large deflection in a symmetrical three-point bending of beam[END_REF] gave analytical and numerical solutions for large deflections of a symmetric three point bending of a simply supported beam subjected to a central concentrated load. Chucheespakul et al. [START_REF] Chucheepsakul | Double Curvature Bending of Variable-Arc-Length Elasticas[END_REF] set up Euler-Bernoulli elastica for pinned-pinned beam. Whereas Magnusson et al. provided the behaviour of the extensible elastica solution for an Euler-Bernoulli beam [START_REF] Magnusson | Behaviour of the extensible elastica solution[END_REF]. The Timoshenko beam was treated by Humer who adapted Reissner's beam approach and gave buckling and postbuckling solutions for cantilever beam subjected to follower force [START_REF] Humer | Exact solutions for the buckling and postbuckling of shear-deformable beams[END_REF][START_REF] Humer | Exact solutions for the buckling and postbuckling of a shear-deformable cantilever subjected to a follower force[END_REF]. In pursuit Batista specified analytical solutions of cantilever beam [START_REF] Batista | Large deflections of shear-deformable cantilever beam subject to a tip follower force[END_REF][START_REF] Batista | Analytical treatment of equilibrium configurations of cantilever under terminal loads using jacobi elliptical functions[END_REF]. Other boundary conditions were presented in [START_REF] Batista | A closed-form solution for reissner planar finite-strain beam using jacobi elliptic functions[END_REF] but no general formulation was presented for general boundary conditions. In this paper, we are interested in straight Timoshenko beam supporting a large and quasi-static plane transformation. The hypotheses are the following: the beam is elastic, isotropic, homogenous and has linear constitutive relation, lastly loading is imposed only at the boundaries. General problem under these hypotheses is presented in a dimensionless form in the first section. In order to embrace wide applications the situations for quasi-static follower or dead load is examined and domain of variation of each dimensionless parameters is examined in the second section. In the third section, our approach leads to a Cauchy initial problem on contrary to most previously conducted studies (based on boundary value problem). This imposes a meticulous analysis (values, variation domain) of each component of the problem. Two invariants of the problem are exhibited. Existence and unicity of the solution for a prescribed load is addressed in the same section. The next section focuses on explicit and analytical solutions of the problem for any given load (force and moment) at one end. The problem of regularity of these solutions in regards to a smooth (and quasi-static) evolution of the load at one end is tackled through a deep analysis of the analytical expressions. After an illustrating example (sec-6) the problem of a pure-shear follower load is presented and shows how asymptotic solutions can be recovered through Taylor expansion of the attainable expression. In the following section the problem of quasi-static stability problem is addressed as a driven parametric oscillator in a general situation. At last the section 7 shows how the proposed approach is able to face problem of quasi-static instability more general than bifurcation (buckling) as catastrophe. The conclusion underlines the main points of the work and proposes some further considerations.

Problem statements

Cosserat formulation

A Cosserat beam model is used [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. i[END_REF]. A material curve C lying in the Euclidian space corresponds to the positions of the center of mass G of each section S. The reference configuration corresponds a stress-free state for which C is a straight segment of length L. A fixed origin O and a Cartesian frame {e i } := (e x , e y , e z ) is chosen such that OG = Se z in the reference configuration. Here S ∈ [0, L] is a material curvilinear coordinate of C. After transformation, the placement of G is defined by the map S → ϕ(S) := OG(S) from [0, L] to the Euclidian space. Note that the material coordinate S always belongs to [0, L] even if the true length of C is changed. The sections S of the beam are supposed to be rigid and normal to C in the reference configuration. Sections are supposed to be uniform (same size and shape). For such a Cosserat-like structure the orientation of the section is prescribed by a moving orthonormal frame {d i } := (d 1 , d 2 , d 3 ) for which d 3 is always normal to the section. In the reference configuration {d i } coincides with the Cartesian frame {e i } and after transformation the directors d i (S) do not depend on the placement of C : in contrary to Euler-Bernoulli model d 3 is not necessarily tangent to the material curve.

75

All along the paper, plane motion in (e x , e z )-plane is considered. Hence for any transformation d 1 and d 3 lie in (e x , e z )-plane and d 2 = e y .

Curvature and strains

A rotation R(S) relates each director to the reference one: d i (S) = R(S)e i . For this plane motion

R =   cos θ 0 sin θ 0 1 0 -sin θ 0 cos θ   (1)
where the angle θ(S) = � e z , d 3 = � e x , d 1 measures in a trigonometric way the rotation of the section around d 2 . Spatial derivation of directors is obtained through (e.g [START_REF] Marrec | Vibration of a timoshenko beam supporting arbitrary large pre-deformation[END_REF]):

∂d i ∂S = κ × d i (2)
where the twist vector (or generalized curvature) κ is:

κ(S) = κ 2 d 2 , with κ 2 (S) = ∂θ ∂S . (3)
According to [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. i[END_REF], the strain measures of the beam are the curvature κ and the generalized strain:

�(S) := ∂ϕ ∂S -d 3 (4)
The directors components of this strain have distinct physical meanings: � • d 1 is the shear strain whereas � • d 3 is the normal strain (by construction � • d 2 = 0). For later convenience, let us introduce

ε(S) := ∂ϕ ∂S = � + d 3 (5)
Introducing directors components ϕ i = ϕ • d i (of course ϕ 2 = 0) straight forward computations show that the components of ε in the moving directors frame are:

ε 1 (S) = ∂ϕ 1 ∂S + ϕ 3 κ 2 , ε 2 (S) = 0, ε 3 (S) = ∂ϕ 3 ∂S -ϕ 1 κ 2 (6)

Internal energy, forces and moments, equilibrium relations

A Kirchhoff-Saint Venant model of isotropic material is used for which the Helmholtz free energy per unit length is quadratic with regard to strain measures:

Ψ = 1 2 AGε 2 1 + 1 2 AE(ε 3 -1) 2 + 1 2 EIκ 2 2 (7)
Hence, the stress resultants depend linearly on the conjugate strains:

N 1 = GA ε 1 shear force N 3 = EA (ε 3 -1) normal force M 2 = EI κ 2 bending moment (8)
Where A and I are the area and the quadratic moment of the section, E and G are the Young modulus and shear modulus G (including eventually a shear correction factor [START_REF] Timoshenko | on the correction for shear of the differential equation for transverse vibrations of prismatic bars[END_REF]) of the beam material. In the directors frame, force and moment vectors are:

N(S) = N 1 d 1 + N 3 d 3 , M(S) = M 2 d 2 .
Equilibrium relations are (e.g. [START_REF] Marrec | Vibration of a timoshenko beam supporting arbitrary large pre-deformation[END_REF]):

∂N ∂S = 0 ∂M ∂S + ε × N = 0 (9)
Using ( 2) and projecting the first equation onto d 1 and d 3 and the seconds onto d 2 one obtains respectively:

GA ∂ε 1 ∂S + EA(ε 3 -1)κ 2 = 0 EA ∂ε 3 ∂S -GAε 1 κ 2 = 0 EI ∂κ 2 ∂S -EAε 1 (ε 3 -1) + GAε 1 ε 3 = 0 (10)
where ( 8) has been used.

Dimensionless form

Dimensionless formulation of the problem is performed by defining first a dimensionless bulk-shear ratio g and the gyration radius of the beam � [START_REF] Marrec | Three-dimensional vibrations of a helically wound cable modeled as a timoshenko rod[END_REF][START_REF] Forgit | A timoshenko-like model for the study of three-dimensional vibrations of an elastic ring of general cross-section[END_REF]:

g = E G , � = � I A (11)
As g � 2(1 + ν) where ν is the Poisson's ratio, 2 � g � 3 for standard material. The gyration radius � has a dimension of a length, in term of magnitude � = O(R), where R is a typical size of the section. The dimensionless curvilinear abscissa and length (slenderness ratio) are:

s = S � , � = L � (12)
By convention let us denote any physical variables previously mentioned v(S) for which a dimensionless companion v(s) can be associated as follow:

ε i (s) = ε i (S), κ i (s) = � κ i (S), ϕ i (s) = 1 � ϕ i (S), θ(s) = θ(S) (13)
Using the fact that ∂f ∂S = 1 � ∂f ∂s and the convention f � := ∂f ∂s for any function f (S) and injecting (13) into (10) one gets:

ε � 1 + g(ε 3 -1)κ 2 = 0 gε � 3 -ε 1 κ 2 = 0 gκ � 2 + ε 1 ε 3 -gε 1 (ε 3 -1) = 0 (14)
Using [START_REF] Rakotomanana | Eléments de dynamique des solides et structures déformables[END_REF], dimensionless energy per unit length [START_REF] Li | Large deflection and rotation of timoshenko beams with frictional end supports under three-point bending[END_REF] becomes:

Ψ = 1 2 ε 2 1 + 1 2 g(ε 3 -1) 2 + 1 2 gκ 2 2 (15)
It is interesting to observe that physical force and moment N(S) and M(S) have also dimensionless companions N(s) and M(s) related by

N = 1 GA N, M = 1 � 1 GA M
In particular, in term of dimensionless components

N 1 = ε 1 N 3 = g(ε 3 -1) M 2 = gκ 2 (16)
the equilibrium equation ( 9) becomes

N � = 0 M � + ε × N = 0 (17)
Consitutive relation [START_REF] Reissner | Some remarks on the problem of column buckling[END_REF] and equilibrium relation [START_REF] Cedolin | Stability of structures: elastic, inelastic, fracture and damage theories[END_REF] (or its projection ( 14)) are our starting point for the analysis of solutions of the static equilibrium problem of beam.

Remark on the boundary conditions

The static problem is posed up to a translation and a rigid rotation. It is then natural to observed that the problem (17) is perfectly described thanks to dynamical (stresses or strains) quantities.

Parametrization of the boundary conditions

The boundary conditions are prescribed at a given end, say s = �. Then consider the known quantities

N � = N(�), M � = M(�) .
Of course, M � = M � d 2 is oriented along e y and N � belongs to the (e x , e z )-plane (equivalently (d 1 , d 3 )plane). Hence orientation and magnitude of N � has to be described properly. Let us define φ � = � d 3 (�), N � the relative angle of the load respectively to the normal d 3 (�) of the last section. By convention φ � is measured in a trigonometric way, such that:

N 1 (�) = N � sin (φ � ) , N 3 (�) = N � cos (φ � ) , (18)
where N � = �N � � and φ � ∈]π, π]. This convention emphasize the crucial role of φ � that prescribes the shear or normal character of the external force at this specific end. For φ � = ±π/2 the external force is a shear, for φ � = π it is a compression, last for φ � = 0 it is a traction.

Remark 3.1 (Cauchy problem). It must be stressed that under this point of view the boundary condition is the set (N � , φ � , M � ): an assignation of the variables (or their derivatives) at a specific end. In that context the problem has the structure of a Cauchy initial value problem from which the solution is known as unique [START_REF] Coddington | Theory of ordinary differential equations[END_REF]. It is clear that for most of the physical problem, distinct constraints are imposed at each end and then the problem has mainly the structure of a Cauchy boundary value problem for which the solution is not necessarily unique. The next sections are mainly related to Cauchy initial value problem associated to the set (N � , φ � , M � ). The last section would be devoted to the Cauchy boundary value problem.

Follower and dead load

At the end s = �, the rotation of the section θ(�) and the angle φ � between the section and the load N � are related by

φ � = φ -θ(�) (19)
According to Fig. 1, φ = � e z , N � is the angle between the external force and the normal of the last section at rest (in such a sense, it can be seen as a Lagrangian quantities).

N � N � θ(s) φ(s) s = � G C S e z e z e x e x d 3 (s) d 1 (�) φ d 1 (s) d 3 (�) φ φ � θ �
Let us suppose that the magnitude N � and M � of the external efforts varies in a quasi-static way. A priori this may imply a possible reorientation of the last section or of the load. Indeed, prescribing the angle φ � during this variation has a large physical impact. More precisely, during this slow evolution, standard situations appear for such special cases:

• If φ � is held constant, the external force is a follower load. In particular if φ � is maintained equal to π the external force is a purely compressive follower load and if φ � = ±π/2 this force acts as a pure shear follower load. Note that during this evolution the value of φ changes in order to respect [START_REF] Dell'isola | Extensible beam models in large deformation under distributed loading: A numerical study on multiplicity of solutions[END_REF].

• If φ is held constant, the external force is a dead load. In particular if φ is maintained null or equal to π the external force is a pure vertical (along e z ) dead load and if φ = ±π/2, this force acts as a pure horizontal dead load. Here φ � has to change during the loading process in order to respect [START_REF] Dell'isola | Extensible beam models in large deformation under distributed loading: A numerical study on multiplicity of solutions[END_REF].

Domain of variation

A priori N � ≥ 0, φ � ∈ [0, 2π[ and M � ∈ R but physical considerations may help to limit these bounds. First of all it can be observed that non-overlapping of the section of the beam during a large bending imposes |κ 2 | � 1/2. Second, it is justified to limit the analysis to a total force lower than a certain multiple of the Euler critical load. For pinned boundary conditions this later is physically P e = EIπ 2 /L 2 , then its dimensionless form is P e = gπ 2 /� 2 . As beam model is justified for slenderness ratio � � 20 and 2 � g � 3, one gets P e � 0.1. Hence it is exhaustive to focus on the following bounds for boundary conditions variables

0 ≤ N � ≤ 0.1, -0.1 ≤ M � ≤ 0.1 . (20)
These bounds are particularly large, especially for elongated structures. It must not be considered as an indicator of the order of magnitude of these variables, but rather as a maximum bounds (physically unreachable in general). However the bounds N � ≤ 0.01 and |M � | ≤ 0.01 are attainable in most of the cases.

These bounds may be used to prescribe some bounds for strains at s = �, and more generally at any s according to [START_REF] Reissner | Some remarks on the problem of column buckling[END_REF]. One get:

-0.1 ≤ ε 1 (s) ≤ 0.1, 0.95 ≤ ε 3 (s) ≤ 1.05, -0.05 ≤ κ 2 (s) ≤ 0.05

Problem analysis

First integration

The system (17) may be easily integrated. Considering the first line, one gets directly

N(s) = N � , ∀ s ∈ [0, �]
In other words, the internal load is constant along the beam. However, as the orientation of the directors d 1 (s) and d 3 (s) are not uniform, the shear and longitudinal component are not uniform along s, but controlled by

N 1 (s) = N � • d 1 (s), N 3 (s) = N � • d 3 (s) (21)
Considering now the second equation of ( 17), its integration is trivial too. First as ε = ϕ � , one gets

M � + ϕ � × N � = 0 then, after integration M(s) -M(0) + (ϕ(s) -ϕ(0)) × N � = 0, ∀ s ∈ [0, �] (22)
In particular, one gets a momentum relation between the bending moment at both ends:

M(0) = M � + � ϕ(�) -ϕ(0) � × N �
Inspired by ( 18) and ( 21), a new variable φ(s) is introduced such that

N 1 (s) = N � sin (φ(s)) N 3 (s) = N � cos (φ(s)) (23)
Then, φ(s) is the relative angle between the normal of the section at s and the external load. Of course φ(�) = φ � . The relation [START_REF] Dell'isola | Extensible beam models in large deformation under distributed loading: A numerical study on multiplicity of solutions[END_REF] can be extended too (see Fig. 1):

φ(s) = φ -θ(s) (24)
As θ � = -φ � , the internal couple M 2 = gθ � along the beam becomes

M 2 (s) = -gφ � (s) (25)
Last all strains may be written in terms of φ according to ( 16):

ε 1 (s) = N � sin (φ(s)) ε 3 (s) = 1 + N � g cos (φ(s)) κ 2 (s) = -φ � (s) (26)
To simplify further notation φ � � will be used instead of φ � (�) in the following.

Non-homogeneous equation

The two first equations of ( 14) are directly satisfied, the last one becomes

g 2 φ �� -gN � sin (φ) + (g -1)N 2 � sin (φ) cos (φ) = 0 (27)
Let us consider that φ �� � = 0. By multiplying ( 27) by 2φ � one obtains after integration

(gφ � ) 2 + 2gN � cos (φ) -(g -1)N 2 � cos 2 (φ) = µ ( 28
)
where µ is a constant related to set (N � , φ � , M � ) of boundary conditions

µ = M 2 � + 2gN � cos (φ � ) -(g -1)N 2 � cos 2 (φ � ) (29)
Remark 4.1 (Invariants). The magnitude N � is the first invariant of the beam as, for any s:

N 2 � = N 1 (s) 2 + N 3 (s) 2
In the first integral (28) of ( 27) the parameter µ appears as a second invariant of the beam configuration. Indeed, in terms of internal load, as M 2 (s) = -gφ � and N 3 (s) = N � cos (φ(s)), one gets all along the beam:

µ = M 2 (s) 2 + 2gN 3 (s) -(g -1)N 3 (s) 2
These two invariants are presented graphically in Fig.

3 The equation ( 28) is a scalar, first order and non-linear ordinary differential equation. Its coefficient is written in terms of the invariants µ and N � of the problem. Its resolution can be performed thanks to some special change of variable that will be developed in the next section. � is prescribed (e.g. [START_REF] Coddington | Theory of ordinary differential equations[END_REF]). Existence and unicity is then obtained if the invariants (µ, N � ) are prescribed. As the load properties

(N � , φ � , M � ) at s = � defines (N � , µ, φ � ) or (N � , µ, φ �
� ) in a unique manner, existence and unicity of the solution are ensured by the prescription of the load at one end.

Remark 4.3 (Discussion on the regularity of the solution). As the set (N � , φ � , M � ) contributes to define both the initial condition and the coefficients µ, N � of the differential equation, the regularity of the solution may be (a prioiri) strongly affected by even a smooth change of N � , φ � or M � . This very special character of the problem explains (at least in a part) the attention of scientists on behavior of beam under large transformation. This problem will be adressed more deeply in section 5.6.

Homogeneous equation

Let us consider the special case φ �� = 0 for all s (then φ � = cste). From ( 27) and according to the domain of variation of g, this situation appears only if N � sin (φ) = 0, what corresponds to two distinct cases:

• If N � = 0 the last end supports only a non-null couple M � . In that case φ(s) = as + b where a and b depends on boundary conditions. Since M � = -gφ � (�), a = -M � /g and imposing arbitrarily φ(�) = 0 implies b = M � �/g and by [START_REF] Dell'isola | Extensible beam models in large deformation under distributed loading: A numerical study on multiplicity of solutions[END_REF] one gets φ = θ(�), therefore by ( 24) one obtains

θ(s) = M � g (s -�) + θ(�) (30)
Dynamical variables are given by N 1 (s) = 0, N 3 (s) = 0 and M 2 (s) = M � . Using ( 6),( 16) one gets:

ϕ 1 (s) = ϕ 1 (0) + g M � (cos (θ(s)) -1) ϕ 3 (s) = ϕ 3 (0) + g M � sin (θ(s)) ( 31
)
This problem is a standard solution in Elastica theory.

• If sin(φ(s)) = 0 then φ(s) = 0 or π: the load is a pure longitudinal force. Therefore by using ( 24) one gets N 1 (s) = 0, N 3 (s) = ±N � and M 2 (s) = 0. So kinematical variables are given for pure traction or compression by:

ϕ 1 (s) = ϕ 1 (0) θ(s) = θ(0) ϕ 3 (s) = ϕ 3 (0) + (1 ± N � g )s
For other situations, φ � is necessarily not uniform and the problem consists in the resolution of (28).

Analysis of µ

As it has been seen µ plays an important role as first variable integration depending on (N � , φ � , M � ). More precisely, the shear load does not affect this parameter that is controlled by M � and the longitudinal part of the force N � cos (φ � ). The behavior of µ with respect to these parameters is given in Fig. 3. In practice µ is positive only for longitudinal traction or eventually moderate longitudinal compressive load associated to non null couple M � . The following bounds for µ are obtained

M 2 � -2gN � -(g -1)N 2 � ≤ µ ≤ M 2 � + 2gN � -(g -1)N 2 � ( 32
)
and -0.5 � µ � 0.5 according to the numerical values proposed in section 3.3. Let us introduce

µ a = -2gN � -(g -1)N 2 � , µ c = 2gN � -(g -1)N 2 � . (33)
Remarks that in contrary to µ, the parameters µ a and µ c depends only on N � . The equality µ = µ a is obtained only for pure compressive load. For a pure traction µ = µ c however the equality µ = µ c may be observed for other configuration for which φ � � = 0 and M � � = 0. There is always µ a ≤ µ but in the special cases for which M � = 0 the bounds are more restrictive: µ a ≤ µ ≤ µ c .

Jacobian elliptic functions

In this section the problem ( 28) is solved using a series of transformations that leads to Jacobian elliptic functions. The study concerns non-homogeneous solutions of (28).

Problem statement

Let us use first the tangent half-angle substitution for φ:

t(s) = tan � φ(s) 2 � ( 34
)
where t(s) is a real valued function. Therefore

cos (φ) = 1 -t 2 1 + t 2 , sin (φ) = 2t 1 + t 2 , φ � = 2t � 1 + t 2 (35)
Injecting ( 35) into [START_REF] Humer | Exact solutions for the buckling and postbuckling of a shear-deformable cantilever subjected to a follower force[END_REF], the differential equation is written as:

t �2 = a t 4 + b t 2 + c (36)
• If a � = 0, the roots of the t-polynomial at the right-hand side can be computed in order to obtain:

t �2 = a (t 2 -α -)(t 2 + α + ) (37)
• If a = 0, then

t �2 = b t 2 + c (38)
All these formulations need a deep analysis of each parameter a, b, c and α ± that are intrinsically related to the set of boundary conditions (N � , φ � , M � ), and more precisely to two independent coefficients µ and N � :

a = µ -µ a 4g 2 b = 2µ + µ a + µ c 4g 2 c = µ -µ c 4g 2 and α + = g+ √ g 2 -(g-1)µ g-1 -N � g+ √ g 2 -(g-1)µ g-1 + N � α -= N � - g- √ g 2 -(g-1)µ g-1 N � + g- √ g 2 -(g-1)µ g-1 (39)

Parameter analysis 145

According to (4.4) µ a ≤ µ then a ≥ 0.

Remark 5.1. Equality a = 0 corresponds to a pure longitudinal compression for which an homogeneous solution has been found already. Notice that if a = 0, then µ = µ a and, as µ a ≤ 0 ≤ µ c , both b and c are strictly negative. In other words, the right hand side of (38) is negative and then no real non-homogeneous solutions t(s) may exists. For pure compressive external load, the homogeneous solution is the only real-150 valued solution.

It is now justified to focus hereafter on [START_REF] Poston | Catastrophe theory and its applications[END_REF] with strictly µ > µ a then a > 0 . However sign of b and c may change according to M � , φ � and N � .

The roots α ± are real as g 2 -(g -1)µ is always positive. The values of α + according to N � and µ is presented in Fig. 4-left (the values of α + are presented in a domain such that µ > µ a ). This roots is strictly positive and close to 1. First order Taylor expansion gives:

α + = 1 - g -1 g N � + O(N 2 � , µ 2 ) (40)
α -= 1 - 2µ µ + 2gN � + O(N � , µ)

Resolution of the elliptic differential equation

As a > 0 and α + > 0 in all the domain of variation, the following formulation of ( 37) is proposed:

� 1 √ aα + t � √ α + � 2 = � ( t √ α + ) 2 + 1 � � ( t √ α + ) 2 - α - α + � (41)
This motivates the following change of variable associated to a rescaling of the curvilinear abscissa

h(ζ) = t(s) √ α + where ζ = √ aα + (s + s 0 ) ( 42
)
where the constant s 0 will be related to the boundary conditions. As

d ds = √ aα + d dζ , (41) becomes: � dh dζ � 2 = (h 2 + 1)(h 2 - α - α + ), (43)
for which solutions are a Jacobian elliptic function [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF]:

h(ζ) = ±cs(ζ | m), with m = 1 + α - α + (
m = 4gN � µ + 2gN � + O(µ, N � ) (45)
This is illustrated in Fig. 5-right, where it can be observed that this approximation is well justified. In a classical formulation of Jacobian elliptic functions, the second argument m belongs in [0, 1]. Here, it is observed that m belongs to ]0, +∞[. Here some analysis is given in order to link the two formulations:

• If µ a < µ < µ c , we have α -> 0 and α + > 0 then m > 1. For that situation, one may use [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF] cs

(ζ | 1 m ) = 1 √ m ds( ζ √ m | m)
to obtain a more standard formulation:

h(ζ) = ± � α + + α - α + ds( � α + + α - α + ζ | α + α + + α - ) then t(s) = ± � α -+ α + ds( � a(α -+ α + )(s + s 0 ) | α + α + + α - ) (46)
• If µ = µ c , α -= 0 and m = 1, then

h(ζ) = ±cs(ζ | 1) = ± 1 sinh(ζ)
Alternatively, this can be obtained by observing that c = 0 then the differential equation [START_REF] Thom | Structural stability and morphogenesis[END_REF] becomes

t �2 = t 2 (at 2 + b). • If µ > µ c we have 0 < m < 1 then t(s) = ± √ α + cs( √ aα + (s + s 0 ) | α + + α - α + ) (47)
This last formulation is always justified if one consider for convention m ∈ R.

In Fig. 6

√ α -)[ if µ a ≤ µ ≤ µ c . ζ -2K 0 2K 4K -5 0 5 cs(ζ| 1
2 ) ds(ζ| 1 2 ) 1/ sinh(ζ) The shifting parameter s 0 has to be adjusted such that

t(�) = tan ( φ � 2 ), 2t � (�) 1 + t(�) 2 = φ � � (48)
where φ � � is just determined by the bending moment according to [START_REF] Chucheepsakul | Double Curvature Bending of Variable-Arc-Length Elasticas[END_REF]:

M � = -gφ � � .
In practice s 0 is determined numerically by an optimisation algorithm associated to the following statement

find s 0 ∈ D such that �t(�) -tan( φ � 2 )� 2 + � 2t � (�) 1 + t(�) 2 -φ � � � 2 = 0 ( 49
)
The domain D and the function t(s) have to be adjusted according to the boundary condition:

• If µ < µ c , D =]0, 1 √ a(α-+α+) 4K( α+ α++α-)[ and t(s) = √ α -+ α + ds( � a(α -+ α + )(s + s 0 ) | α+ α++α-) (the sign -is not necessary as ds(ζ + 2K(m) | m) = -ds(ζ | m)). • If µ > µ c and φ � � < 0, then D =]0, 1 √ aα+ 2K( α++α- α+ )[ and t(s) = √ α + cs( √ aα + (s + s 0 ) | α++α- α+ ) • If µ > µ c and φ � � > 0, then D =]0, 1 √ aα+ 2K( α++α- α+ )[ too but t(s) = - √ α + cs( √ aα + (s + s 0 ) | α++α- α+ )
Indeed with such domain D, the discussion associated to Fig. 6 shows that the solution s 0 in D of the problem (49) is unique in all the cases.

Regularity of the solutions

The existence and unicity of the problem [START_REF] Cedolin | Stability of structures: elastic, inelastic, fracture and damage theories[END_REF] has been already underlined in rem-4.2. However regularity of the solution relatively to a smooth change of the boundary conditions parameters N � , φ � , M � was just interrogated in rem.4.3. According to the discussion of Fig. 6 the function t(s) is regularly varying according to the parameters α ± , a and s 0 in all the variation domain of µ, N � and φ � . It is then the case for φ(s) too. In other words φ(s) is regularly dependent to the set (µ, N � , φ � ) or (µ, N � , φ � � ). In one hand the map (N � , M � , φ � ) → (µ, N � , φ � ) is smooth and surjective, in the other hand (N 1 (s), N 3 (s), M 2 (s)) are regular function of φ(s) thanks to [START_REF] Magnusson | Behaviour of the extensible elastica solution[END_REF]. This leads to the following important result:

The solutions N 1 (s), N 3 (s), M 2 (s) (or equivalently ε 1 (s), ε 3 (s), κ 2 (s)) of the problem (17) are regularly dependent to the boundary conditions (N � , φ � , M � ). No bifurcation may occurs if the load is completely controlled (prescription of the couple intensity, force intensity and of the orientation of the force relatively to the section) in a smooth and quasi-static way at one end of the beam. Reciprocally, some bifurcation may occur only if one of these load properties are held free (standard buckling problem of beam by a dead-load) or if kinematical variable (position or orientation of the section) are controlled. This result is valid for any large transformation, any shape of the beam and any isotropic and elastic material. In term of instabilities, the preceding result shows that quasi-static instabilities cannot appear if the loads are controlled at one end in a smooth way. This is the case for follower loads. However, as this work do not invoque dynamical effect, this analysis do not allow any conjecture on dynamical-instabilities. In particular fluttering effect is over the scope of the paper.

At this stage all unknown functions and parameter are determined. It is then possible to have explicit expression of strains and shape of the beam without any approximation. This is now illustrated through a practical example.

First illustrating example

Let us consider a beam of length � = 50 and material ratio g = 5/2 supporting the following boundary conditions N � = 0.01, M � = 0.05, φ � = 3π/4. The beam is glued at s = 0 on the origin of the Cartesian frame, then θ(0) = 0 and ϕ(0) = 0.

Determination of each parameters

According to the set (N � , M � , φ � ) one obtains µ � -0.03, µ a � -0.05, µ c � 0.05 and m � 5.86 then µ a < µ < µ c and m > 1. The solution t(s) corresponds to (46) but can still be written as t(s) = √ α + cs( √ aα + (s + s 0 ) | m), if one accepts that m > 1. In Fig. 7 t(�) is presented as a function of s 0 .

Intersections with the level set tan (φ � /2) respecting the constraints (48) are highlighted by a dot. The

s 0 -100 -50 0 50 100 -4 -3 -2 -1 0 1 2 3 4 t(ℓ) tan( φℓ 2 )
s ℓ 0 possible values of s 0 form a periodic set but any value may be chosen as position of the dot and t(s) share the same periodicity. This illustrates that the definition of a restrictive domain D in (49) does not affect the generality of the resolution.

Determination of the internal forces and moments

According to [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF], ( 25) and ( 35)

N 1 (s) = N � 2t 1 + t 2 , N 3 (s) = N � 1 -t 2 1 + t 2 , M 2 (s) = -g 2t � 1 + t 2 (50)
that is easily computed as t(s) is given in (47) and because the derivation rule d dz cs(z|m) = -ns(z|m)ds(z|m) holds true for any m ∈ R [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF]. One obtains t � (s) = -√ aα + ns(ζ | m)ds(ζ | m) and then:

N 1 (s) = N � 2 √ α + cs(ζ | m) 1 + α + cs 2 (ζ | m) , N 3 (s) = N � 1 -α + cs 2 (ζ | m) 1 + α + cs 2 (ζ | m) , M 2 (s) = √ ag 2α + ns(ζ | m)ds(ζ | m) 1 + α + cs 2 (ζ | m) , ζ = √ aα + (s + s 0 ) m = α + + α - α + (51)
The internal forces and moment are illustrated in Fig. 8-left.

Determination of the rotation and placement

The rotation θ(s) of the section is directely obtained by ( 24) as φ(s) is now-determined univocally by (34) as φ(s) = 2 arctan (t(s)). Then

θ(s) = φ -2 arctan � t (s)

�

The constant φ can be chosen with a great liberty for a given set of boundary conditions as the static problem is unchanged up to a rigid rotation. For a quasi-static loading, details have been exposed in sec.3.2. As the beam is glued at s = 0 in this example, one imposes θ(0) = 0 and then φ = 2 arctan (t(0)). Result of this typical case are presented in Fig. 8-right.

The situation is more complex for the placement function ϕ(s). The two components of this vector have to be determined. In practice ϕ(s) belongs to the (e x , e z )-plane. The choice of a proper orthonormal basis is crucial. Cartesian frame can be chosen such that ϕ = ϕ x (s)e x + ϕ z (s)e z or a mobile directors frame such that ϕ = ϕ 1 (s)d 1 (s) + ϕ 3 (s)d 3 (s). However an other orthogonal frame is naturally introduced into the problem, let us denote it (e t , e y , e n ), where

e n = N � �N � � , e t = e y × e n
It is a fixed frame induced by the direction of the external force N � = N � e n . Of course this direction is known if the external force is perfectly described. In this frame the placement is expressed by ϕ = ϕ t (s)e t +ϕ n (s)e n . Note that � d 3 (s), e n = φ(s), then e n • d 3 (s) = cos (φ(s)), and e n • d 1 (s) = sin (φ(s)). The relation [START_REF] Baker | Jacobi elliptic functions and the complete solution to the bead on the hoop problem[END_REF] gives

M(s) -M(�) + (ϕ(s) -ϕ(�)) × N � = 0 that is written in this frame: M 2 (s) -M � - (ϕ t (s) -ϕ t (�))N � = 0, hence ϕ t (s) = ϕ t (�) + M 2 (s) -M � N �
Remark that M 2 (0) is now determined, hence the following expression holds too:

ϕ t (s) = ϕ t (0) + M 2 (s) -M 2 (0) N � (52)
and seems more appropriate, as in general ϕ(0) is generally imposed. In order to determine ϕ n (s) a more tricky strategy is necessary as an integration is needed. Starting with (5) that is written as

ϕ � n e n + ϕ � t e t = ε 1 (s)d 1 (s) + ε 3 (s)d 3 (s)
After projection along e n :

ϕ � n = ε 1 (s) sin (φ(s)) + ε 3 (s) cos (φ(s)) = N 1 (s) sin (φ(s)) + ( 1 g N 3 + 1) cos (φ(s)) = N � � 2t 1 + t 2 � 2 + N � g � 1 -t 2 1 + t 2 � 2 + 1 -t 2 1 + t 2
where ( 16), ( 50) and ( 35) have been used successively. Integration is then performed as

ϕ n (s) -ϕ n (0) = � s 0 ϕ � n (σ) dσ which implies ϕ n (s) = ϕ n (0) + � s 0 N � � 2t 1 + t 2 � 2 + N � g � 1 -t 2 1 + t 2 � 2 + 1 -t 2 1 + t 2 ds ( 53
)
where

t(s) = � |α + |cs( � a|α + |(s + s 0 ) | m
) is perfectly known. For this general case, no simple explicit formulation of a primitive of ϕ � n can be obtain (a simplified expression is given in sec.7 for a particular case). However, the smooth behavior of such function allows us to use a simple integration technic in order to obtain a numerical solution. In practice one uses a rectangular integration on a fine discretization of the length (Δs = �/1 000). As mentioned earlier, ϕ n (0) = ϕ t (0) = 0 is considered. The placement function is now completely determined as ϕ(s) = ϕ t (s)e t + ϕ n (s)e n . However if a Cartesian frame is privileged one has of course to remind that the angle between e z and N � satisfies θ(�) + φ � = φ (see Fig. 1). Then This result is plotted on Fig. 9 where the sections (oriented along d 1 (s)) are also presented for the sake of the clarity. Knowing that the physical radius of the section is �, these sections are plotted in a dimensionless way with a unitary radius. Hence this picture is completely dimensionless with a respected slenderness ratio. The beam (� = 50, g = 5/2) is still glued at s = 0, then θ(0) = 0 and support at the other end a pure shear-follower load. Then the boundary conditions at s = � are: N � � = 0, φ � = π/2 and M � = 0. For such type of control all the configurations presented on this section ( 7) are quasi-statically stable according to the sec-5.6. The objective of this example is to analyse the qualitative and quantitative behavior of the beam as N � increases. Some simplified asymptotic expressions are given in the case of moderate shear force N � .

Parameter analysis

For this example µ = 0. In particular, according to [START_REF] Poston | Catastrophe theory and its applications[END_REF], one gets α -= 1 then α -> 0, t(s) is given by (46). One has explicitly first order Taylor (40) can be used. In the same spirit, Taylor expansion of √ α -+ α + , � a(α -+ α + ) and α+ α++α-in terms of N � exhibits a first order approximation:

α + = ( 2g g-1 -N � )/( 2g g-1 + N � ) but if 0 < N � � 1 (
t(s) = √ 2 ds( � N � g (s + s 0 ) | 1 2 ) + O(N � )
Within this approximation, one obtains directly M 2 (s) by ( 50) and ϕ t (s) thanks to (52):

M 2 (s) � � 2gN � cn( � N � g (s + s 0 ) | 1 2 ) ϕ t (s) � ϕ t (0) + � 2g N � cn( � N � g (s + s 0 ) | 1 2 )
The expression of ϕ n is highly simplified by applying leading term approximation with respect to N � . In fact by (53) one gets:

ϕ n (s) � ϕ n (0) + � s 0 1 -t 2 1 + t 2 ds
Integration is explicit in such a case and gives

ϕ n (s) � ϕ n (0) + s -2 � g N � E( � N � g (s + s 0 ) | 1 2 )
where E(x | m) is the Jacobian epsilon function and sn(s) := sn(

� N � g (s + s 0 ) | 1 2 ) 225 7.

Qualitative and quantitative analysis

The figure [START_REF] Antman | Nonlinear problems of elasticity[END_REF] represents successive configurations as the magnitude of the transverse shear load increases. These 20 simulations have been computed with the exact formulation (in particular numerical integration of (53)). The computation cost for these 20 simulations is 0.23 s on a 2.2 GHz Intel Core i7 (without plotting curves). This analytical approach clearly allows real-time simulations. On Fig. 10 one observes qualitatively that the beam wrinkles with a periodicity controlled by the magnitude of the force. The asymptotic approach developed in the preceding section provides some tools in order to characterize this wrinkle. The periodicity P (along s) of ϕ is given by the complete elliptic integral of the first kind K(m) that is the quarter period of the Jacobian elliptic functions too. According to the approximation proposed in the preceding section, the periodicity of ϕ is

z x N ℓ = 0.15 N ℓ = 0.2 N ℓ = 0.1 N ℓ = 0.05 N ℓ = 0.01
P = 4 � g N � K( 1 2 ) � 7.41 � g N �
Hence the number of wrinkle along the beam is �/P . In one hand, the size A of the wrinkle is related to the magnitude of ϕ t : A = 2 � 2g/N � . In the other hand the spatial periodicity B (indeed P is the material periodicity) of this wrinkle may be computed as B = |ϕ n (P )ϕ n (0)|. Focusing on the leading terms, and detailing the computation, on obtain:

B = � � � � P -2 � g N � E(4K( 1 2 ) | 1 2 )
� � � � = � � � � 4 � g N � K( 1 2 ) -8 � g N � E( 1 2 )
� � � � = 4 � g N � � 2E( 1 2 ) -K( 1 2 )
� = 2π K( 1 2 ) � g N �
where the complete elliptic integral E is introduced thanks to the relation

E(4K( 1 2 ) | 1 2 ) = 4E( 1 2 ) ([ 23
]-22.16.29) and [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF]-19.7.1 has been used too. Lastly the size ratio B/A (see Fig. 11) of the wrinkle pattern is independent of any material, geometrical or loading parameter:

B A = π √ 2K( 1 2 ) ∼ 1.2
This properties of the deformed shape can been seen as a particular signature of such sollicitation by a pure-shear load.

Quasi-static stability

Let us consider a given static configuration. Let us denote the associated quantities V (for example 230 N � , θ(s), d i (s)). This quantities are supposed to be known, they are solutions of the problem [START_REF] Cedolin | Stability of structures: elastic, inelastic, fracture and damage theories[END_REF] for given boundary condition. The quasi-static stability problem consists of the analysis of the behavior of a perturbed solution V = V + δV keeping invariant some boundary conditions. The perturbation δV is infinitesimal, then all quadratic terms in δV would be neglected in the following, leading to a linear analysis of the unknowns δV. Let us focus first on the vectorial quantities N(s) = N(s) + δN(s), M(s) = M(s) + δM(s) and ε(s) = ε(s) + δε(s). The prescribed initial quantities V(s) satisfies [START_REF] Cedolin | Stability of structures: elastic, inelastic, fracture and damage theories[END_REF]. In particular N(s) = N � is a fixed vector. The perturbed quantities solve (17) too. They writes (linear δ-order):

δN � (s) = 0 δM � (s) + ε(s) × δN(s) + δε(s) × N � = 0 (54)
According to the first line δN(s) is a constant vector, so let not it δN � (= δN(s)) where δN � is indeed prescribed by boundary conditions. It is then natural to focus on

δM � (s) + ε(s) × δN � + δε(s) × N � = 0 (55)
Consideri now the components of vectorial quantities. Constitutive relations motivate the use of initial (respectively current) basis for the prescribed (respectively perturbed) configuration:

N � = N 1 (s) d 1 (s) + N 3 (s) d 3 (s) M(s) = M 2 (s) d 2 ε(s) = ε 1 (s) d 1 (s) + ε 3 (s) d 3 (s)
with

N 1 (s) = N � sin φ(s) N 3 (s) = N � cos φ(s) M 2 (s) = g θ � (s) (56)
N(s) = N 1 (s) d 1 (s) + N 3 (s) d 3 (s) M(s) = M 2 (s) d 2 ε(s) = ε 1 (s) d 1 (s) + ε 3 (s) d 3 (s)
with

N 1 (s) = ε 1 (s) N 3 (s) = g (ε 3 (s) -1) M 2 (s) = g θ � (s) (57)
Last, the perturbations are arbitrary defined on the new basis:

δN � = δN 1 (s) d 1 (s) + δN 3 (s) d 3 (s) δM(s) = δM 2 (s) d 2 δε(s) = δε 1 (s) d 1 (s) + δε 3 (s) d 3 (s)
Even if d 2 = d 2 = e y , this is not the case for the other directors that may undergo an infinitesimal rotation: θ(s) = θ(s) + δθ(s). Then d 1 (s) = cos(δθ(s)) d 1 (s)sin(δθ(s)) d 3 (s) and d 3 (s) = sin(δθ(s)) d 1 (s) + cos(δθ(s)) d 3 (s) may be approximated in this first order approach by: d 1 (s) = d 1 (s)δθ(s) d 3 (s) and d 3 (s) = d 3 (s) + δθ(s) d 1 (s) respectively. The constitutive relation for the perturbation of the moment is trivial δM 2 (s) = g δθ � (s) but this is not the case for the strains. This is obtained by considering that N i (s) may be defined by two equivalent ways:

for

N 1 (s) (N � + δN � ) • d 1 (s) = (ε(s) + δε(s)) • d 1 (s) for N 3 (s) (N � + δN � ) • d 3 (s) = g ((ε(s) + δε(s)) • d 3 (s) -1) (58)
Expanding each side and using first δ-order approximations, one obtains, respectively:

N 1 (s) -δθ(s) N 3 (s) + δN 1 (s) = ε 1 (s) -δθ(s) ε 3 (s) + δε 1 (s) N 3 (s) + δθ(s) N 1 (s) + δN 3 (s) = g (ε 3 (s) + δθ(s) ε 1 (s) + δε 3 (s) -1) (59)
Using (56), a constitutive relation depending on the initial state is obtained for the perturbation:

δε 1 (s) = cos (φ(�) -φ(s))δN 1 (�) -sin (φ(�) -φ(s))δN 3 (�) - � g-1 g N 3 (s) -1 � δθ(s) δε 3 (s) = 1 g � sin (φ(�) -φ(s))δN 1 (�) + cos (φ(�) -φ(s))δN 3 (�) � -g-1 g N 1 (s) δθ(s) (60)
where φ(�) + θ(�) = φ(s) + θ(s) (= φ) has been used. Hence strain perturbations δε i (s) are all related to the field of micro-rotation δθ(s) and controlled by both the prescribed configuration and boundary conditions imposed on the force perturbation δN � at the last section. This observation induces that (55) may be written merely in terms of a single degree of freedom δθ(s) of the perturbation with parameters controlled by the prescribed configuration and boundary conditions. Straight forward calculation gives:

δθ �� (s) + k 2 (s) δθ(s) = f (s) (61)
where k 2 (s) and f (s) are

k 2 (s) = - 1 g N 3 (s) + g -1 g 2 � N 2 3 (s) -N 2 1 (s) � f (s) = 1 g � g -1 g N 3 (s) -1 � � cos (φ(�) -φ(s))δN 1 (�) -sin (φ(�) -φ(s))δN 3 (�) � + . . . . . . g -1 g 2 N 1 (s) � sin (φ(�) -φ(s))δN 1 (�) + cos (φ(�) -φ(s))δN 3 (�) � (62)
δθ(s) is a solution of a linear, non-homogeneous, second order differential equation with non-constant coefficients. This differential equation is of the class of driven parametric oscillators for which analysis is well-documented but beyond the scope of the present paper. Boundary conditions affect only the non-homogeneous term. Indeed, the prescription of δN � is of particular physical importance. This is detailed in the next paragraph for dead and follower load cases.

Dead-load

For a dead-load on gets N � = N � then δN � = 0. The equation (62) becomes:

δθ �� (s) + k 2 (s) δθ(s) = 0 (63)
Boundary conditions at σ = 0 or � are either δθ(σ) = 0 for constant orientation of the section or δθ � (σ) = 0 if the moment kept constant.

Exemple for a pure longitudinal compression. For this initial configuration φ(s) = π and k(s) is constant:

k 2 = N � g + g -1 g 2 N 2 � ( 64
)
In that case k 2 is strictly positive, hence solutions of (63) are of the form

δθ = C 1 cos (ks) + C 2 sin (ks) (65)
for which the constants C i depend on the boundary conditions.

• For simply supported beam δθ � (0) = 0 and δθ � (�) = 0. Trivial solution C i = 0 is imposed except if k = nπ/� (with n ∈ N * ) for which δθ(s) = C 1 cos (ks) is a possible solution. A straightforward computation shows that this buckling solution occurs for

N � = g 2(g -1) �� 1 + 4(g -1)( nπ � ) 2 -1 � (66)
This non-linear relation may be simplified in first approximation as nπ � � 1. This leads to the standard Euler critical-load N � = g( nπ �

) 2 for buckling of simply supported beam.

• For clamped-hinged beam δθ(0) = 0 and δθ � (�) = 0. Again trivial solution is imposed except if k = (2n + 1)π/(2�) (with n ∈ N) for which δθ(s) = C 2 sin (ks) is a possible solution. In that case the critical load becomes:

N � = g 2(g -1) �� 1 + (g -1)( (2n + 1)π � ) 2 -1 � (67)
which corresponds in first approximations to the standard Euler critical load N � = g( (2n + 1)π 2� ) 2 for cantilever beam.

• For beam clamped at both ends, δθ(0) = 0, δθ(�) = 0, the same approach leads to (66) as for simply 250 supported beam.

Follower load

If the force acting on the last section is a follower load, then N � � = N � . However the components of the actual force at this end are unchanged: N 1 (�) = N 1 (�) and N 3 (�) = N 3 (�). Applying (58) at s = � leads to:

δN 1 (�) = N 3 (�) δθ(�), δN 3 (�) = -N 1 (�) δθ(�) . (68)
Injecting these expressions in (62), the equation (61) writes after rearrangement:

δθ �� (s) + k 2 (s) � δθ(s) -δθ(�) � = 0 (69)
where δφ(s) + δθ(s) = δφ(�) + δθ(�) has been used. For a follower-load, the angle between the normal of the last section and the external force is the same before and after perturbation, then δφ(�) = 0. In other words: δθ(s)δθ(�) = -δφ(s). According to this change of variables, the above differential equation becomes:

δφ �� (s) + k 2 (s) δφ(s) = 0 (70)
where the boundary condition at s = � is already prescribed δφ(�) = 0 even if δθ(�) is unknown (note that δθ(�) may be a boundary condition prescription).

Exemple for a pure longitudinal compression. As already mentioned k 2 is, in that case, a positive constant specified in (64). As δφ(�) = 0, solution of (70) is of the form

δφ(s) = C 1 sin (k(� -s)) (71)
However, the boundary conditions imposed on the beam still play an important role, as it is highlighted by the three following examples:

• If the couple remains null at the last end δM 2 (�) = 0 then δθ � (�) = 0 or equivalently δφ � (�) = 0. Direct calculation shows that no solution δφ(s) may respect this condition except the trivial ones: δφ(s) = 0, then δθ(s) = δθ(�). Physically speaking any perturbation of the beam induced an instability in the form of a rigid rotation. Note that this observation is valid for any constraint imposed on s = 0. In particular if the beam is clamped at the origin δθ(s) = 0: no transverse perturbation are possible.

• If the beam is simply supported at s = 0 but δθ � (�) � = 0 at the other end, the problem writes in terms of δφ: δφ � (0) = 0 and δφ(�) = 0. These boundary conditions imposes non-trivial solution if k = (2n + 1)π/(2�) (with n ∈ N). A buckling appears for a critical load presented in (67) and the mode of buckling may be written as

δθ(s) = δθ(�) -C 1 sin � (2n + 1)π 2� (� -s) �
where the constant C 1 is arbitrary.

• If δθ(0) = 0 the beam is clamped at s = 0 then δφ(0) = δθ(�) and still δφ(�) = 0. According to (71), C 1 sin (k�) = δθ(�). Let consider moderate compression (0 < k� < π) in that case the solution (71) is always possible and writes:

δφ(s) = δθ(�) sin (k(� -s)) sin (k�) , then δθ(s) = δθ(�) � 1 - sin (k(� -s)) sin (k�)

�

Then, for any infinitesimal perturbation of the orientation of the last section (and then of the follower load), a perturbed solution exists. No brutal bifurcation is observed. In conclusion the clamped beam is particularly unstable under the loading of a follower compressive force (even if this latter is infinitesimal).

Boundary problem

Examining now that the same beam (� = 50, g = 5/2) supports a pure dead-load N � = N � e z , M � = 0. Of course φ � � = 0 as M � = 0 however φ � is not prescribed. Indeed, the last boundary condition is the orientation of the section at s = 0: θ 0 = θ(0). Let us define in the same spirit M 0 := M 2 (0), φ 0 := φ(0) and 24) θ 0 = φφ 0 and according to the section (3.2) φ = 0. The questions are the following:

φ � 0 := φ � (0) (= -M 0 /g). From (
What is the configuration of the beam as θ 0 varies according to a command, for a fixed pure dead-load ?

Is this solution unique ?

In order to address these questions, one first consider the initial condition problem which has unique solution:

What is the configuration of the beam as φ � varies for a prescribed N � and such that M � = 0, φ = 0 ?

In a second step the map of the solution associated to θ 0 ∈ C → φ � is studied according to the command C of the boundary condition.

Parameter analysis

First observe that µ = 2gN � cos (φ � ) -(g -1)N 2 � cos 2 (φ � ) then µ a ≤ µ ≤ µ c and t(s) is of the form (46). Second, as φ � � = 0 we have t � (�) = 0 and according to Fig. 6 and the associated analysis:

� a(α -+ α + )(� + s 0 ) = +K, if 0 ≤ φ � < π � a(α -+ α + )(� + s 0 ) = -K, if π ≤ φ � < 2π , where K := K( α + α + + α - )
Then, for a fixed φ � , s 0 can be determined directly. In other words, t(s) is perfectly determined and then φ(s). In particular θ 0 := -2 arctan (t(0)) becomes:

θ 0 = -2 arctan � � α -+ α + ds( � a(α -+ α + )s 0 | α + α + + α - ) � = 2 arctan � � α -+ α + ds(±K + � a(α -+ α + )� | α + α + + α - ) � = ±2 arctan � √ α - cn( � a(α -+ α + )� | α+ α++α-) �
According to the preceding results and using translation rules [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF]

-22.4.(iii). The sign ± is -if φ � ∈ [0, π[ and + if φ � ∈ [π, 2π[.
Recall that a and α± depend explicitly on φ � according to (39) and (29).

Catastrophic instablities

The evolution of θ 0 versus φ � is presented in Fig. 12. The map φ � → θ 0 is bijective only for moderate magnitude N � of the dead-load. For moderate N � the angle φ � is uniquely determined for a given θ 0 hence the associated configuration of the beam is unique. For larger N � , more than one value of φ � is associated to a given value of θ 0 . This lack of uniqueness induces a non-uniqueness of the configurations for a fixed θ 0 . Supposing that the command C is "θ 0 increases from 0". The initial state is then associated to the point (φ � = 0, θ 0 = 0) for which the beam supports a pure traction. Afterward the graph in Fig. 12 has to be read from right to left and the first configurations are defined by the increasing branch initiated at (0, 0). For large N � , this branch has an inflection point (for φ � � -π/3 if N � = 0.01 see in Fig. 13-left). In order to respect the command imposed on θ 0 , the values of φ � has to present a large jump (its magnitude is higher than 3π/2 if N � = 0.01). This jump corresponds to a large and brutal change of the configuration of the beam as it is observed in Fig. 13-right. This phenomena is known as a catastrophe in instability theories ( [START_REF] Thom | Structural stability and morphogenesis[END_REF][START_REF] Poston | Catastrophe theory and its applications[END_REF][START_REF] Golubitsky | An introduction to catastrophe theory and its applications[END_REF]).

φ ℓ -2π -7 4 π -3 2 π -5 4 π -π -3 4 π -1 2 π -1 4 π 0 θ 0 0 π/4 π/2 3π
The analytical and dimensionless approach followed along the present paper allows to tackle this type of problem through explicit solution without a priori approximations or hypotheses. This methodology may be seen as a complementary approach to those more usual based on the internal energy.

Conclusion

Analytical study of large (but plane) transformation of a Timoshenko beam with linear constitutive law was conducted. Cosserat's formulation was used by modeling the beam as a curvilinear line with moving directors frame. Inspired by [START_REF] Marrec | Vibration of a timoshenko beam supporting arbitrary large pre-deformation[END_REF] and by adapting dimensionless procedures, equilibrium relations in nondimensional form were found. Even if the physical problem is defined by a mixture of kinematical (position and orientation of the crosssection) and dynamical (force and moments) boundary conditions, the formulation followed in the present paper chooses to emphasize first a Cauchy problem formalism where all strains are prescribed at one end (in contrast to most studies in this field). More precisely, at one end, the moment and force intensities are supposed to be known but the orientation of the load with respect to the normal of the section too. This approach ensures existence and unicity of the solution. It had been proved that such prescription of the moment and force (intensity and orientation relatively to the cross-section) at a given end ensures regularity of the solution with respect to the boundary conditions too. Explicit formulation of the solution is obtained thanks to Jacobian elliptic functions for which coefficients are smoothly dependents on the invariants of the problem. This methodology let us tackle several theoretical and physical problems since explicit solutions were obtained without any approximation.

Several examples are then presented to illustrate this methodology. In the case of a pure-shear follower load, load is completely prescribed at one end, and a quasi-static evolution is straightforwardly addressed. In such a case, the explicit analysis is the opportunity to exhibit an universel size ratio of a wrinkle pattern: this size-ratio is independent of the material and the geometrical properties of the beam but independent of the intensity of the load too.

For mixed boundary conditions, the perturbation of the problem is presented in a general way, which allows us to compute the perturbed solutions under various boundary conditions. This general perturbed problem written as a driven parametric oscillator for which each (space-varying) parameters is explicitly available.

In the last example, the kinematic control of a beam supporting a dead-load is presented. For such a problem, the usual perturbation methodology could be used to detect eventual instability. The present paper however focus on an alternative nonlinear method that depicts the whole equilibrium solution. Emergence of even more critical instability than fork bifurcation (as usual buckling) is observed. Indeed such so called catastrophe remains reachable in a straightforward and explicit way.

A more detailed study of energy function property and a dynamic approach could improve our understanding of planar elastic beam. In particular dynamic instabilities, such as fluttering, is one of the interesting issues that may be potentially addressed by extending the present work.

  80

Figure 1 :

 1 Figure 1: Parametrization of a current configuration of a Timoshenko beam. At a given curvilinear abscissa s, the center of mass G and the section S are given. The directors d 1 (s) and d 3 (s) of this section are obtained by a rotation θ(s) around ey. At s = � the external force N � has an angle φ(�) = φ � with the normal d 3 (�) of the last section and φ with ez.

Figure 2 :

 2 Figure 2: Representation of the two invariants (N � in green and µ in purple) in the configuration space (N 1 (s), M 2 (s), N 3 (s)). The solutions N 1 (s), M 2 (s) and N 3 (s) of the problem is along the intersection of the two surfaces.

  120

Remark 4 . 2 (

 42 Existence and unicity of the solution). According to Cauchy theorem, the ordinary differential equation (28) has a unique solution if φ � or φ �

Figure 3 :

 3 Figure 3: Variation of µ according to the parameter of the boundary conditions.

µ

  a (N ℓ ) µ c (N ℓ )

Figure 4 :

 4 Figure 4: Variation of α + (left) and log 10 (|α -|) (right) according to N � and µ. The curves µ = µa(N � ) and µ = µa(N � ) are presented too.The values of α -according to N � and µ is presented in Fig.4-right. As α -has a large domain of variation the log 10 (|α -|) is plotted and the sign of α -is specified in each domain. In practice α -> 0 if µ a < µ < µ c and α -< 0 if µ c < µ. α -= 0 if µ = µ c and for small values of µ, N � :

44) 5 . 4 .

 54 Class of solutionsParticular attention must be drawn to the values of m presented in Fig.5-left. For moderate parameters µ and N � :

Figure 5 :

 5 Figure 5: Left: Variation of log 10 (m) according to N � and µ. The curves corresponding to µa(N � ) and µc(N � ) are presented too. Right: Variation of m and 1/m according to µ for fixed N � . The approximation (45) is presented too. The vertical lines depict the position of µa(N � ) (left) and µc(N � ) (right).

  155

  the function cs(ζ | m) and ds(ζ | m) are presented for m = 1/2. These odd functions have distinct periodicity. Introducing the complete elliptic integral of the first kind K(m), the periodicity is 2K for cs(ζ | m) and 4K for ds(ζ | m). Remarks that K(m) → ∞ as m → 1. As m passes through 1 the transition between cs(ζ | m) to ds(ζ | m) through 1/ sinh (ζ) is then smooth as the non-periodic function 1/ sinh (ζ) is asymptotically considered as a periodic function with infinite period. This anodyne remark leads to a more crucial ones detailed in sec.5.6. Last, it is observed that ds(ζ | m) belongs to ] -∞, -√ 1m] ∪ [ √ 1m, +∞[ if 0 < m < 1. According to (34) and (46), this means that |φ(s)| does not belong to the interval [0, 2 arctan(

Figure 6 :

 6 Figure 6: cs(ζ | m), ds(ζ | m) and 1/ sinh (ζ) for m = 1/2. Numerically K(m) � 1.85. The horizontal dashed lines corresponds to ± √ 1m � 0.71

Figure 7 :

 7 Figure 7: Graphs of t(�) according to s 0 . The level-set tan (φ � /2) is presented too. The graph of t(s) for an admissible s 0 is presented in red (the corresponding abscissa is given in red too).

Figure 8 :

 8 Figure 8: Left: N 1 (s), N 3 (s) and M 2 (s) for the example. Right: θ(s).

  e n • e z = cos φ, e n • e x = sin φ, e t • e z =sin φ, e t • e x = cos φ. The Cartesian components ϕ x = ϕ • e x and ϕ y = ϕ • e y become: ϕ x (s) = cos ( φ) ϕ t (s) + sin ( φ) ϕ n (s), ϕ y (s) = cos ( φ) ϕ n (s)sin ( φ) ϕ t (s) .

215 7 .

 7 Pure-shear follower load

Figure 9 :

 9 Figure 9: Dimensionless deformed shape of the beam glued at s = 0 and supporting a follower load N � making an angle φ � = 3π/4 with the normal of the last section (s = �). The intensity of the force is N � = 0.01 and a bending moment M � = 0.05 is imposed on this last section too. The length of the beam is � = 50. No extra-amplification is used except for the force N � .

Figure 10 :

 10 Figure 10: Dimensionless deformed shape of the beam glued at s = 0 and supporting a follower shear load N � at s = � (φ � = π/2 and M � = 0). The properties of the beam are � = 50 and g = 5/2. The figure represents successive snap shot for N � = [0.01, 0.2] (with regular step of 0.01). Some values or N � are given and the directions of the external load are represented too in order to help the reader.

Figure 11 :

 11 Figure 11: dimensionless deformed shape of the beam glued at s = 0 and supporting a pure shear load at s = � (N � = 0.01, φ � = π/2 and M � = 0). The properties of the beam are � = 500 and g = 5/2. The wrinkles are clearly visible. For this loading B � 53.6 and A � 44.7 and P � 117.3. There is 4.25 wrinkles.

Figure 12 :

 12 Figure 12: Evolution of θ 0 according to φ � for various magnitude of the dead-load (� = 50, g = 5/2, N � = N � ez, M � = 0, φ = 0). If C is θ 0 increases linearly from 0, the graph has to be read from right to left.

Figure 13 :

 13 Figure 13: Left: evolution of θ 0 according to φ � if N � = 0.01. In orange bold the history of the φ � and θ 0 for the command C. The arrows help to read the quasi static evolution and underline the gap associated to the catastrophe. Right: the successive configurations for various values of θ 0 respecting the command C. The colors of the dots in the left figure are associated to the colors of the corresponding configurations. The purple arrow highlights the brutal transition of the configurations during the catastrophe.

marwan.hariz94@gmail.com

loic.lemarrec@univ-rennes1.fr

jlerbet@gmail.com

Acknowledgements

This work was supported in part by a grant from the Walid Joumblatt Foundation For University Studies and the Henri Lebesgue Center. The authors gratefully acknowledge the support of IRMAR -CNRS UMR 6625 for providing financial support. The authors are very grateful to Professor Lalaonirina Rakotomanana for the careful and thoughtful comments on the work.