Chedid Saade

Philippe Le Grognec
email: philippe.le_grognec@ensta-bretagne.fr

Maël Couchaux

Mohammed Hjiaj

Elastoplastic bifurcation analysis for the lateral-torsional buckling of straight beams

Keywords: Lateral-torsional buckling, Elastoplasticity, Analytical modeling, Finite element validation, Rectangular cross-section/I-beam

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The buckling problem is one of the most common failure mode of slender (or thin) structures. Buckling analyses are thus required when designing beam-type structures, so as to derive the so-called critical loads and corresponding bifurcation modes. Among possible buckling phenomena encountered in practice, the buckling of a straight beam under axial compression is certainly the most widely studied in the literature. In elasticity, the first critical load of a compressed beam was found as the solution of an eigenvalue problem by Euler at the beginning of the eighteenth century. In contrast, the pioneering works dealing with the plastic buckling problem only date back to the end of the nineteenth century. The first significant result is due to Shanley [START_REF] Shanley | Inelastic column theory[END_REF] and concerns the so-called Shanley's column under axial compressive force, which turns out to be a discrete model supposedly able to reproduce qualitatively the behavior of a beam cross-section. He provided the first satisfactory answer regarding the nature and the value of the first critical load in elastoplasticity. In both cases of a discrete or continuous structure, the first bifurcation is proved to occur at the so-called tangent modulus critical load, giving rise to an incipient unloaded zone and an increasing load during the initial post-bifurcation stage. Later, it was shown by Cimetière [START_REF] Cimetière | On the modelling and buckling of elastoplastic plates (Sur la modélisation et le flambage des plaques élastoplastiques[END_REF] that in the case of compressed rectangular plates, elastoplastic structures are likely to bifurcate Preprint submitted to International Journal of Solids and Structures November 3, 2022

J o u r n a l P r e -p r o o f

Journal Pre-proof within intervals of critical loads, by continuously modifying the unloaded zone at critical time and thus the resulting structural stiffness. Each continuous range of bifurcation points (one per mode) observed in plasticity spreads from the tangent modulus critical value to the elastic one. Whilst the bifurcation can theoretically take place at any point of the continua, it naturally occurs at the tangent modulus critical value in practice, which corresponds to the minimum buckling load and manifests itself by incipient unloading only. Moreover, when dealing with plastic buckling, different critical values may be obtained, depending on the retained plasticity theory. A large discrepancy between the results provided by the flow and deformation theories has been emphasized by many authors in the literature. In most cases, it turns out that the flow theory largely overpredicts the experimental critical values, whereas the deformation theory provides critical load values that are in much better agreement, although the latter does not include the elastic unloading possibility. This major issue, known as the plastic buckling paradox, may be due to unavoidable imperfections in experiments. Having this in mind, Liu [START_REF] Liu | FE analysis of plastic buckling of plates with initial imperfections and simulation of experiments[END_REF] analyzed the imperfection sensitivity of uniaxially compressed plates by means of finite element computations, considering both the flow and deformation theories for comparison purposes. Neale [START_REF] Neale | Effect of imperfections on the plastic buckling of rectangular plates[END_REF] carried out a theoretical analysis of the influence of the imperfections on the plastic buckling of rectangular plates using a Reissner-type variational principle. He showed the close agreement between the results obtained with the flow theory and those obtained with the deformation theory (and thus also with experimental data), provided that imperfections were conveniently accounted for in the formulation. Besides, the plastic buckling paradox has been partially explained by the fact that the flow theory induces an elastic shear modulus at onset of buckling. It is particularly detrimental, for instance, to beam torsional buckling problems, for which the torsional stiffness is not affected by the plasticity developing during uniaxial compression.

As an alternative solution, Becque [START_REF] Becque | The application of plastic flow theory to inelastic column buckling[END_REF] defined, in the context of the flow theory, a new inelastic shear modulus to be used instead of the elastic one in the torsional buckling analysis The elastic lateral-torsional buckling of beams has been widely studied over the past decades. The particular case of I-beams was intensively investigated. The first analytical expression of the elastic critical bending moment has been proposed by Timoshenko and Gere [START_REF] Timoshenko | Theory of Elastic Stability[END_REF] for I-beams resting on fork supports. Vlasov [START_REF] Vlasov | Thin-Walled Elastic Beams[END_REF] introduced the warping effect necessary for the analysis of torsion of open sections and thus for the investigation of their lateral-torsional buckling. The elastic lateral-torsional buckling phenomenon depends on several factors such

J o u r n a l P r e -p r o o f

Journal Pre-proof as the bending moment distribution [START_REF] Serna | Equivalent uniform moment factors for lateraltorsional buckling of steel members[END_REF], the relative position/orientation of the applied force with respect to the shear center, the boundary conditions [START_REF] Piotrowski | Lateral-torsional buckling of steel beams elastically restrained at the support nodes[END_REF] and the cross-section symmetries [START_REF] Wagner | Torsion and buckling of open sections[END_REF][START_REF] Mohri | Déversement des poutres en I sous chargements asymétriques[END_REF]. The elastic critical bending moment is generally calculated with the three factor formula used in most design codes [START_REF] Afnor | NF EN 1993-1-1/NA, Eurocode 3 -Design of steel structures -Part 1-1: General rules and rules for buildings (French National Annex to EN 1993-1-1[END_REF][START_REF]CEN, EN 1999-1-1, Eurocode 9 -Design of aluminium structures -Part 1-1: General structural rules[END_REF]. It can be derived from differential equations [START_REF] Timoshenko | Theory of Elastic Stability[END_REF] or from the potential energy using the Rayleigh-Ritz method [START_REF] Timoshenko | Theory of Elastic Stability[END_REF][START_REF] Ritz | Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik[END_REF]. Numerous finite element models have also been developed so as to ease its practical calculation [15]. However, this elastic critical bending moment is often not sufficient to design profiles against this buckling phenomenon as initial imperfections and residual stresses may accelerate the process. The current design codes [START_REF]CEN, EN 1993-1-1, Eurocode 3 -Design of steel structures -Part 1-1: General rules and rules for buildings[END_REF] integrate these effects using classical buckling curves [START_REF] Rondal | Formulation d'Ayrton-Perry pour le flambement des barres métalliques[END_REF]. Taras

and Greiner [START_REF] Taras | New design curves for lateral-torsional buckling -Proposal based on a consistent derivation[END_REF] developed a more consistent approach for the lateral-torsional buckling analysis integrating explicitly the initial imperfections in an Ayrton-Perry formulation [START_REF] Ayrton | On struts[END_REF].

This approach, based on an elastic behavior of the member, is not adequate when the applied bending moment induces yielding of the cross-section.

Albeit the elastic lateral-torsional buckling phenomenon is well known and closed-form solutions of the elastic critical values are available for many cross-section shapes, the case of an elastoplastic beam has rarely been addressed so far. As in the case of axial compression, when considering sufficiently short beams (with low slenderness), the lateral-torsional buckling phenomenon is prone to arise after the onset of plasticity, and the use of elastic critical bending moments is no more consistent since the yielding of the cross-section should also be considered when calculating the bifurcation critical value. The determination of elastoplastic critical moments can be thus of great importance for design purposes. Neal [START_REF] Neal | The lateral instability of yielded mild steel beams of rectangular cross-section[END_REF] may be considered to have pioneered the lateral instability analysis of inelastic beams. He examined the lateral-torsional buckling (under pure or simple bending) of partially yielded beams with a rectangular cross-section, both from a theoretical and experimental point of view. He specially demonstrated and verified experimentally that the elastic torsional rigidity should be retained in the evaluation of the critical bending moment even in the elastoplastic regime.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Based on this pioneering work, the case of an I-beam was investigated by Galambos [START_REF] Galambos | Inelastic lateral buckling of beams[END_REF] a few years later from an engineering point of view. The influence of initial residual stresses was included in the analysis and not only the lateral (bending) stiffness but also the torsional stiffness was reduced due to yielding. Even later, Trahair and Kitipornchai [START_REF] Trahair | Buckling of inelastic I-beams under uniform moment[END_REF][START_REF] Kitipornchai | Inelastic buckling of simply supported steel I-beams[END_REF] also analyzed the lateral-torsional buckling of inelastic I-beams, from a theoretical and experimental point of view, still considering the effect of residual stresses. They showed that the use of a null tangent modulus and an elastic shear modulus did not greatly affect the critical moments of the beam in most cases. They also highlighted the possible significant influence of initial geometric imperfections.

In line with these works, this paper is devoted to the derivation of analytical solutions for the critical moment of elastoplastic lateral-torsional buckling of straight beams with various cross-section shapes. In Section 2, a general theoretical formulation of the lateral-torsional elastoplastic buckling problem is developed. It is based on the 3D elastoplastic bifurcation theory, using the J 2 flow theory considering the von Mises yield criterion with linear isotropic hardening. This 3D approach has already been applied successfully to elastic/plastic Timoshenko beams [START_REF] Le Grognec | On the plastic bifurcation and post-bifurcation of axially compressed beams[END_REF], two-layer beams with partial interaction [START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF][START_REF] Grognec | Plastic bifurcation analysis of a two-layer shear-deformable beam-column with partial interaction[END_REF], or sandwich beams [START_REF] Douville | Exact analytical solutions for the local and global buckling of sandwich beam-columns under various loadings[END_REF][START_REF] Le Grognec | Elastoplastic buckling and post-buckling analysis of sandwich columns[END_REF] under axial compression/pure bending, but also to the lateral buckling of offshore pipelines including twisting effects [START_REF] Grognec | Investigation of the torsional effects on the lateral buckling[END_REF]. In Section 3, two particular cases are investigated:

a beam with a thin rectangular cross-section and an I-beam. In some particular conditions, closed-form solutions are derived for the elastoplastic critical bending moments. Next, in Section 4, an in-house shell finite element program is employed to validate the previous analytical solutions. Some special features of the finite element tool are first recalled and numerical computations are performed for the two specific geometries defined above. Numerical computations are also carried out using a beam finite element formulation in Abaqus software, for comparison purposes. The analytical and numerical results are generally proved to be in very good accordance. Lastly, general conclusions and prospects are proposed in J o u r n a l P r e -p r o o f

Journal Pre-proof Section 5.

Analytical formulation of the elastoplastic lateral-torsional buckling problem

Problem definition

Let us consider a straight beam with an arbitrary uniform cross-section. For simplicity purposes, we will restrict ourselves to the case of a bi-symmetric cross-section, for which the center of torsion and the center of gravity coincide. The geometry of the beam is defined by its length L along the x-axis (with endpoints located at x = 0 and x = L) and its total depth 2h along the y-axis (with the neutral axis defined at y = 0). The cross-section area, the flexural first and second moments of area with respect to axes y and z, the torsional moment of area and the warping moment will be defined in due time. The material is assumed to be homogeneous, displaying a linear elastic regime followed by a plastic regime driven by the von Mises yield criterion with a linear isotropic hardening. In the following developments, use will be made of the Young's modulus E, the shear modulus G = E 2(1+ν) (involving the Poisson's ratio ν), the yield stress σ y and the tangent elastoplastic modulus E T = EH E+H (where H stands for the (constant) hardening modulus).

The beam is subjected to two opposite bending moments M about the z-axis acting at the beam ends in such a way that it undergoes a pure bending state with a uniform bending moment along the length. This particular loading will lead to the lateral-torsional buckling of the beam in the elastic or elastoplastic regime, depending on the geometric and material parameters. Simply-supported boundary conditions are assumed so as to provide a closed-form expression for the buckling mode. The critical load and the bifurcation mode will be derived here from a general 3D approach which has been already applied to other beam buckling problems [START_REF] Le Grognec | On the plastic bifurcation and post-bifurcation of axially compressed beams[END_REF][START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF][START_REF] Grognec | Plastic bifurcation analysis of a two-layer shear-deformable beam-column with partial interaction[END_REF][START_REF] Douville | Exact analytical solutions for the local and global buckling of sandwich beam-columns under various loadings[END_REF][START_REF] Le Grognec | Elastoplastic buckling and post-buckling analysis of sandwich columns[END_REF][START_REF] Grognec | Investigation of the torsional effects on the lateral buckling[END_REF].

J o u r n a l P r e -p r o o f

Journal Pre-proof

Theoretical formulation

The theory is developed first in a 3D framework (the beam is initially seen as a 3D body) using a total Lagrangian formulation. From a general point of view, the critical loading Λ c and the bifurcation mode X of a 3D body, occupying Ω in the reference configuration, can be classically obtained by solving the following bifurcation equation [30]:

∀ δu, Ω ∇ T δu : K i (Λ c ) : ∇X dΩ = 0 (1)
where δu can be considered as a kinematically admissible test function or, more physically, as the virtual variation of the unknown displacement field u.

In Equation ( 1), K i represents the fourth-order nominal tangent tensor and can be written as follows:

K i = ∂Π ∂F = F. ∂Σ ∂E .F T + (I.Σ) T = F.D i .F T + (I.Σ) T (2)
In the above equation, E denotes the Green strain tensor and Σ the second Kirchhoff stress tensor (symmetric); F is the deformation gradient and Π = F.Σ the first Kirchhoff stress tensor (non-symmetric); I represents the fourth-order unit tensor (I ijkl = δ il δ kj ) and the superscript T the transposition of a second-order tensor and the major transposition of a fourth-order tensor ( A T ijkl = A klij ), respectively. Finally, D i is the fourth-order material tangent tensor reflecting the 3D constitutive law (with i = e or p, depending on whether the behavior is elastic or plastic, respectively).

In the present case, the plastic zone at critical time (namely the region which plastifies at the buckling moment) does not extend to the whole structure, as it does for a column under uniform axial compressive force. However, the minimum critical value encountered in practice will also correspond here to the situation for which the plastic zone corresponding to the bifurcated solution is equal to that of the fundamental solution. In brief, the bifurcation still takes place at the tangent modulus critical load with incipient unloading.

J o u r n a l P r e -p r o o f

Journal Pre-proof

In the elastic regime, the material is linear and isotropic, and D e is the classical elasticity operator whose components in an orthonormal basis are D e ijkl = λδ ij δ kl + µ (δ ik δ jl + δ il δ kj ), where δ ij is the Kronecker symbol, and λ = Eν (1+ν)(1-2ν) and µ = G = E 2(1+ν) are the Lamé constants.

In the plastic regime, the plastic threshold is defined by the von Mises yield function with a linear isotropic hardening:

f (Σ, A) = 3 2 Σ d : Σ d -σ y -A A = Hp (3)
where Σ d denotes the deviatoric part of Σ and p the equivalent plastic strain. The material tangent elastoplastic tensor reads then:

D p = ∂Σ ∂E = D e - D e : ∂f ∂Σ ⊗ ∂f ∂Σ : D e H + ∂f ∂Σ : D e : ∂f ∂Σ (4)
where the tensor product ⊗ of two second-order tensors S and T is defined by

(S ⊗ T) ijkl = S ij T kl .
More explicit expressions of the above tensors shall now be derived by exploiting the uniaxial stress state in the beam at hand. Indeed, in the pre-critical state, the beam is under pure bending and thus the nominal (first Kirchhoff) stress tensor comprises a single non-zero component, namely the longitudinal tensile/compressive stress Π xx .

Then, the following assumption of small pre-critical deformations, which is usually satisfied in practice, can be made:

∇u 1 (5)
Thus, the stress tensor Σ writes:

Σ = F -1 .Π ≈ Π (6)
and is also limited to a single component Σ xx .

The nominal tangent tensor in Equation ( 2) is also simplified as follows:

K i ≈ ∂Σ ∂E + (I.Σ) T = D i + Σ xx e i ⊗ e x ⊗ e x ⊗ e i (7)
where the implicit summation convention on repeated indices is used with i = x, y, z.

In the case of plasticity, the material tangent tensor in Equation ( 4) takes the following explicit form:

D p = D e - µ 2 H + 3µ (I -3e x ⊗ e x ) ⊗ (I -3e x ⊗ e x ) (8)
related to the uniaxial stress state.

The components of D p in the orthonormal basis (e x , e y , e z ) are:

D p xxxx = λ + 2µ -4µ 2 H+3µ D p yyzz = λ -µ 2 H+3µ D p yyyy = λ + 2µ -µ 2 H+3µ D p xxzz = λ + 2µ 2 H+3µ D p zzzz = λ + 2µ -µ 2 H+3µ D p xxyy = λ + 2µ 2 H+3µ D p xyxy = D p xzxz = D p yzyz = µ (9)
The other components are either zero or derived from Equation (9) using both major and minor symmetries of tensor D p (D p ijkl = D p klij = D p jikl = D p ijlk ). Furthermore, when dealing with 1D beam-like models, ad hoc assumptions are usually added in order to enforce some specific stress state in the body. Namely, the transverse normal material stresses are assumed to be zero: Σ yy = Σ zz = 0. Taking into account these assumptions, one may replace the tensor D i in Equation ( 7) with the reduced tensor C i defined as:

C i ijkl = D i ijkl + D i ijyy (D i yyzz D i zzkl -D i zzzz D i yykl )+D i ijzz (D i zzyy D i yykl -D i yyyy D i zzkl ) D i yyyy D i zzzz -D i yyzz D i zzyy (i, j) = (y, y), (z, z) (k, l) = (y, y), (z, z) (10)
J o u r n a l P r e -p r o o f

Journal Pre-proof

It can be readily checked that tensor C i has the major and both minor symmetries. In the sequel, only the following reduced moduli (and their equivalents obtained by major or minor symmetries) will be needed:

C e xxxx = E C p xxxx = E T C e xyxy = C e xzxz = C e yzyz = µ C p xyxy = C p xzxz = C p yzyz = µ (11)
Considering the above simplifications, the bifurcation equation ( 1) reads:

∀ δu, Ω ∇ T δu : C i + Σ c xx e i ⊗ e x ⊗ e x ⊗ e i : ∇X dΩ = 0 ( 12
)
where Σ c xx stands for the uniaxial stress at critical time.

Let us now consider the bending-twisting response of the beam. As regards the bending problem, the beam is assumed to be sufficiently long so that transverse shear effects can be neglected. The Euler-Bernoulli kinematics is thus retained, involving only three scalar displacement fields u(x), v(x) and w(x), respectively the axial and transverse displacements along y-and z-axes of the centroid axis of the beam (the flexural rotations are directly related to the transverse displacement fields since the cross-sections are supposed to remain orthogonal to the centroid axis of the beam after deformation). As for the twisting problem, the cross-section torsional rotation θ(x) (about x-axis) is introduced together with the warping function ω(y, z) which reflects the distribution of axial displacement throughout the section due to torsion. Once the beam buckles from the fundamental solution to a benttwisted shape, the expressions for the bifurcation mode X and the displacement variation δu are both chosen according to the following kinematics, bringing together the displacement

J o u r n a l P r e -p r o o f

Journal Pre-proof fields related to the flexural and torsional effects:

X = U -yV, x -zW, x +ωΘ, x V -zΘ W + yΘ δu = δu -yδv, x -zδw, x +ωδθ, x δv -zδθ δw + yδθ (13)
The bifurcation mode gradient then writes:

∇X =       U, x -yV, xx -zW, xx +ωΘ, xx -V, x +ω, y Θ, x -W, x +ω, z Θ, x V, x -zΘ, x 0 -Θ W, x +yΘ, x Θ 0       (14)
and the gradient of the displacement variation takes a similar form.

The bifurcation equation becomes thus:

∀ δu, δv, δw, δθ, Ω [E i (U, x -yV, xx -zW, xx +ωΘ, xx ) (δu, x -yδv, xx -zδw, xx +ωδθ, xx ) +µ (-V, x +ω, y Θ, x ) (-δv, x +ω, y δθ, x ) + µ (V, x -zΘ, x ) (δv, x -zδθ, x ) +µ (-V, x +ω, y Θ, x ) (δv, x -zδθ, x ) + µ (V, x -zΘ, x ) (-δv, x +ω, y δθ, x ) +µ (-W, x +ω, z Θ, x ) (-δw, x +ω, z δθ, x ) + µ (W, x +yΘ, x ) (δw, x +yδθ, x ) +µ (-W, x +ω, z Θ, x ) (δw, x +yδθ, x ) + µ (W, x +yΘ, x ) (-δw, x +ω, z δθ, x ) +Σ c xx (U, x -yV, xx -zW, xx +ωΘ, xx ) (δu, x -yδv, xx -zδw, xx +ωδθ, xx ) +Σ c xx (V, x -zΘ, x ) (δv, x -zδθ, x ) + Σ c xx (W, x +yΘ, x ) (δw, x +yδθ, x )] dΩ = 0 ( 15
)
where E i = E or E T , depending on whether the local behavior is elastic or elastoplastic, respectively.

Under an axial compressive force, plasticity occurs instantaneously and evolves uniformly over the whole structure whereas under a bending moment, plastic zones appear first at the top and bottom faces of the cross-section and spread progressively towards the neutral axis

J o u r n a l P r e -p r o o f

Journal Pre-proof of the beam. At the sought critical point (which will be assumed to correspond to the tangent modulus critical value), the beam will then be partially plasticized with an elastic central zone and two plastic zones located at the cross-section extremities, separated from each other by transition heights. In the above equation, the uniaxial stress field is thus not uniform and, due to plasticity, it cannot even be simply related to the bending moment.

Let us consider only the upper half of the beam (y > 0) by symmetry and let c (> 0) be the y-coordinate of the transition between the elastic and plastic zones, measured from the neutral axis (see Figure 2). Under pure bending and the assumption of small pre-critical deformations, the axial strain writes simply:

E xx = -yK ( 16
)
where K is the uniform curvature, which is supposed to be positive along the length of the beam. Owing to the uniaxial stress state, the longitudinal stress can be deduced easily:

Σ xx = E E xx = -yEK = -y c σ y in elasticity (0 < y < c) Σ xx = -σ y + E T E xx + σy E = -σ y 1 + E T E y c -1 in plasticity (c < y < h) (17)
Let us now define the reduced cross-section area A(ξ), first and second flexural moments of area Z(ξ), I y (ξ) and I z (ξ), and warping moment I ω (ξ), as follows: Thanks to the above definitions, the bending moment can be expressed in the following concise manner:

A(ξ) = 2 ξ 0 dS Z(ξ) = 2 ξ 0 y dS I y (ξ) = 2 ξ 0 z 2 dS I z (ξ) = 2 ξ 0 y 2 dS I ω (ξ) = 2 ξ 0 ω 2 dS
M = 2 h 0 -yΣ xx dS = 2 c 0 y 2 σy c dS + 2 h c yσ y 1 -E T E dS + 2 h c y 2 σy c E T E dS = σy c I z (c) + σ y 1 -E T E (Z(h) -Z(c)) + σy c E T E (I z (h) -I z (c)) (19)
The bifurcation equation (15) can also be integrated over the cross-section using the same definitions [START_REF] Taras | New design curves for lateral-torsional buckling -Proposal based on a consistent derivation[END_REF]. During this step, negligible higher-order terms will be eliminated (assuming that Σ c xx E) in such a way that the bifurcation equation can be simplified into:

∀ δu, δv, δw, δθ, L 0 [(EA(c) + E T (A(h) -A(c))) U, x δu, x + (EI z (c) + E T (I z (h) -I z (c))) V, xx δv, xx + (EI y (c) + E T (I y (h) -I y (c))) W, xx δw, xx + (EI ω (c) + E T (I ω (h) -I ω (c))) Θ, xx δθ, xx +µJΘ, x δθ, x -σy c I z (c) (Θ, x δw, x +W, x δθ, x ) -σ y 1 -E T E (Z(h) -Z(c)) (Θ, x δw, x +W, x δθ, x ) -σy c E T E (I z (h) -I z (c)) (Θ, x δw, x +W, x δθ, x ) dx = 0 ( 20
)
where J = S (y + ω, z ) 2 + (zω, y ) 2 dS corresponds to the torsional inertia of the crosssection.

Finally, integrating by parts Equation ( 20) with respect to x yields four local differential equations for the components U , V , W and Θ of the eigenmode:

(EA) * U, xx = 0 (EI z ) * V, xxxx = 0 (EI y ) * W, xxxx +M cr Θ, xx = 0 (EI ω ) * Θ, xxxx -GJΘ, xx +M cr W, xx = 0 (21)
where:

(EA) * = EA(c) + E T (A(h) -A(c)) (EI z ) * = EI z (c) + E T (I z (h) -I z (c)) (EI y ) * = EI y (c) + E T (I y (h) -I y (c)) (EI ω ) * = EI ω (c) + E T (I ω (h) -I ω (c)) (22)
and M cr is the critical bending moment given by Equation [START_REF] Ayrton | On struts[END_REF] for the critical value of c which will act further as the bifurcation parameter.

The natural boundary conditions are also provided by the integration by parts, together with the local differential equations. Let us first define the displacement boundary conditions fulfilled by the eigenmode. The beam is simply-supported at both ends, so that V (0) =

J o u r n a l P r e -p r o o f

Journal Pre-proof

W (0) = V (L) = W (L) = 0.
The torsional rotation is also prevented at both ends: Θ(0) = Θ(L) = 0. Last, the longitudinal displacement is fixed at an arbitrary point (say U (0) = 0) so as to prevent rigid body modes. Assuming a kinematically admissible virtual displacement δu, the conditions δv(0) = δw(0) = δv(L) = δw(L) = 0, δθ(0) = δθ(L) = 0 and δu(0) = 0 must be taken into account during the integration process leading to the remaining stress boundary conditions:

(EA) * U, x (L) = 0 (EI z ) * V, xx (0) = (EI z ) * V, xx (L) = 0 (EI y ) * W, xx (0) = (EI y ) * W, xx (L) = 0 (EI ω ) * Θ, xx (0) = (EI ω ) * Θ, xx (L) = 0 (23)

Solution procedure

The first two uncoupled equations in [START_REF] Galambos | Inelastic lateral buckling of beams[END_REF] together with the corresponding boundary conditions give rise to the following trivial solutions: U (x) = V (x) = 0. The last two equations in [START_REF] Galambos | Inelastic lateral buckling of beams[END_REF] are coupled since they both involve the lateral displacement field W and the rotation Θ. The sought buckling modes are thus prone to combine bending and torsion.

Nevertheless, these two equations can be combined into a single differential equation only involving Θ:

(EI ω ) * Θ, xxxxxx -GJΘ, xxxx - M 2 cr (EI y ) * Θ, xx = 0 (24)
At this time, by solving Equation ( 24) together with the corresponding boundary conditions, it can be shown that the first mode takes necessarily the following form:

Θ(x) = Θ 0 sin πx L (25)
and that W (x) is also found to be a sine function (with a similar form).

In the following case studies, the plastic zones within the cross-section will be such that

J o u r n a l P r e -p r o o f

Journal Pre-proof the warping inertia (EI ω ) * will be negligible compared to the remaining terms. Equation [START_REF] Le Grognec | On the plastic bifurcation and post-bifurcation of axially compressed beams[END_REF] can thus be rewritten as follows:

M cr = π L GJ(EI y ) * (26)
Equation ( 26) provides the expression of the sought critical bending moment. In elasticity, the same developments can be made. By replacing the elastoplastic tangent modulus E T with the elastic modulus E, one obtains the classical expression for the critical bending moment of an elastic beam in lateral-torsional buckling (as far as the warping inertia is still negligible):

M e cr = π L GJEI y ( 27
)
where the second moment of area of the whole section I y (h) has been simply renamed here I y , for convenience purposes.

In elastoplasticity, Equation ( 26) is implicit, as M cr depends on the transition coordinate c as well as (EI y ) * . This equation will be solved, analytically or numerically, so as to give the critical value of c and, subsequently, the critical value of the bending moment M cr in the elastoplastic regime.

Application cases

Thin rectangular cross-section

Let us first consider the case of a beam with a thin rectangular cross-section of depth 2h measured along y-axis and width b measured along z-axis. It is assumed that b h which corresponds to the configuration of a thin-walled open section whose torsional inertia is classically estimated through the following expression:

J = 2b 3 h 3 ( 28
)
J o u r n a l P r e -p r o o f

Journal Pre-proof

The useful partial moments of area and the elastoplastic flexural stiffness (EI y ) * take on the following values:

Z(ξ) = bξ 2 I y (ξ) = b 3 ξ 6 I z (ξ) = 2bξ 3 3 (EI y ) * = (Ec + E T (h -c)) b 3 6 ( 29
)
whereas the warping inertia is neglected here due to the thinness of the beam.

To facilitate further analysis, two dimensionless parameters are introduced:

a = c h
(relative coordinate of the transition between elastic and plastic zones)

η = E T E
(ratio between the tangent elastoplastic and elastic moduli) (30)

Finally, inserting the expression [START_REF] Ayrton | On struts[END_REF] of the critical moment into Equation ( 26), one obtains a sixth-order polynomial equation for the critical value of a:

(η-1) 2 a 6 -6(η-1) 2 a 4 +(4η+β)(η-1)a 3 +(9η 2 -(β+18)η+9)a 2 -12η(η-1)a+4η 2 = 0 (31) with β = π 2 EGb 4 L 2 h 2 σ 2 y .
In the general case, this equation must be solved numerically. If the geometric and material parameters are such that the lateral-torsional buckling occurs after the onset of plasticity, it can be found that Equation [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF] has only one real root between 0 and 1. Using this particular root, the value of the critical moment is then obtained in a straightforward manner.

In the case of a relatively small tangent modulus (namely for η 1), the influence of hardening can be neglected, as it will be verified in the sequel. One can therefore consider

J o u r n a l P r e -p r o o f

Journal Pre-proof that η = 0 and Equation ( 31) becomes:

a 4 -6a 2 -βa + 9 = 0 (32)
Such a quartic equation can be analytically solved, see Appendix A. Without going into details, closed-form solutions are given below for the relative transition coordinate a and the critical moment M cr for this particular case:

a = α 3 - √ α 3 (-α 3 +12α+2β) 2α 2 M cr = πb 3 h 3L √ EGa = M e cr √ a (33)
with:

α = γ 2 + 32 γ + 4 γ = 3 4β 2 + 512 + 4 β 4 + 256β 2 (34)
In both cases, one can observe that the heterogeneous elastoplastic behavior (and stress state) clearly complicates the calculations and leads to critical moments that depend not only on the tangent modulus but also on the yield stress, among other parameters.

I-beam

Let us now consider the case of a beam with a I-shaped cross-section whose dimensions are represented in Figure 3. The depth and thickness of the vertical web are respectively 2h and b. The horizontal flanges are similar (so as to ensure the bi-symmetry of the beam), with a width w and a thickness t. It is assumed that b, t h, w. In this configuration of thin-walled open section, the torsional inertia is classically estimated through the following expression:

J = 2 3 b 3 h + wt 3 (35)
At the critical time, the beam is supposed to be partially plasticized assuming that the

Z(ξ) = bξ 2 for ξ ∈ [0, h[ and Z(h) = bh 2 + 2wth I y (ξ) = b 3 ξ 6 for ξ ∈ [0, h[ and I y (h) = 1 6 (b 3 h + w 3 t) I z (ξ) = 2bξ 3 3 for ξ ∈ [0, h[ and I z (h) = 2bh 3 3 + 2wth 2 (EI y ) * = Ecb 3 6 + E T (b 3 (h-c)+w 3 t) 6 (36)
In this new configuration, the warping inertia may not be negligible due to the warping displacement in the flanges. Nevertheless, as the flanges are supposed to behave plastically at the critical time, with a tangent modulus much lower than the elastic one, it will practically have no influence on the critical moment.

Finally, inserting the expression [START_REF] Ayrton | On struts[END_REF] of the critical moment into Equation ( 26), one obtains again a sixth-order polynomial equation for the critical value of a = c/h, which reads:

-L 2 b 2 h 4 [(η -1)a 3 -3(η -1)a + 2η] 2 + 12L 2 btwh 3 [(η -1)a -η] [(η -1)a 3 -3(η -1)a + 2η] -h 2 [(η -1)a -η] [π 2 a 2 b 6 ρ + 36L 2 t 2 w 2 (η -1)a -36L 2 t 2 w 2 η] -π 2 a 2 b 3 twhρ [t 2 (η -1)a -η(t 2 + w 2 )] + π 2 a 2 w 4 t 4 ρη = 0 (37)
with η = E T /E like in the previous case, and ρ = EG/σ 2 y .

In the general case, this equation must be solved numerically and allows one to derive the critical value of the bending moment. Again, considering a null tangent modulus, Equation ( 37) can be simplified into the following quartic equation:

L 2 b 2 h 3 a 2 -3 2 -12L 2 btwh 2 a 2 -3 -h π 2 ab 6 ρ -36L 2 t 2 w 2 -π 2 ab 3 t 3 wρ = 0 (38)
The closed-form solution of this equation will not be given explicitly as it involves too cumbersome expressions. In addition, it can be shown that, for a given non zero tangent modulus, the above approximate solution is in general less accurate for an I-section than for a rectangular cross-section, due to the presence of flanges.

Numerical results and finite element validation

Finite element modeling

The previous theoretical solutions will be now validated thanks to finite element computations. For that purpose, use is made of an in-house shell finite element program. Developed within a total Lagrangian framework, the program can handle finite plasticity and geometric non-linearities (finite displacements and rotations) in order to deal with the plastic buckling and advanced post-critical behavior of general structures. The use of such an advanced program brings many advantages, among which:

• it is based on a shell finite element model which will validate the beam formulation four fundamental steps are given below:

1. At the end of each increment, it must be checked whether one has passed across one or several critical points. The detection of critical points is based on the singularity of the tangent stiffness matrix, which is factorized following the Crout formula. The critical points are determined by counting the negative pivot number.

2.

Each critical point detected has to be isolated in order to determine its nature: limit point or bifurcation point. To do this, the current arc-length is re-estimated several times using a dichotomy-like method. In the case of prescribed loads, a simple way to distinguish between a limit point and a bifurcation point is to determine the sign of the so-called current stiffness parameter introduced in [START_REF] Bergan | Solution techniques for non-linear finite element problem[END_REF], which changes when passing a limit point but remains unchanged when passing a bifurcation point.

3.

If it is a bifurcation point, the step increment size is reevaluated so as to reach a point just behind the bifurcation point.

4.

Finally, the switching on a bifurcated branch is performed using the mode injection method [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]: at the first step of a bifurcating branch, the eigenvector corresponding to the zero eigenvalue is computed and used as an alternative prediction.

For more information, the whole program is described in full details in Le Grognec and

Le van [START_REF] Le Grognec | Elastoplastic bifurcation and collapse of axially loaded cylindrical shells[END_REF].

Alternatively, finite element computations will also be performed using Abaqus software, for comparison purposes. Again, incremental calculations are needed so as to take into account the material non-linearities, and use is also made of arc-length (Riks) methods in order to handle possible limit points. In addition, imperfections are introduced so as to trigger the buckling phenomenon and follow the post-critical equilibrium path. To ease implementation, a perturbation force technique is adopted rather than other alternatives including geometric imperfections. The resulting equilibrium curves obtained with such a computation are supposed to be slightly degenerated, when compared to the idealized curves

J o u r n a l P r e -p r o o f

Journal Pre-proof stemming from perfect structures. Moreover, no bifurcation point is explicitly detected with such calculations, but the corresponding equilibrium curves classically display limit points, which can be identified as the sought critical values, failing another solution. Numerical estimates of the critical moment will be obtained using three-node (quadratic) beam finite elements with six degrees-of-freedom per node (B32 elements in Abaqus). The use of B33

elements should have been more consistent with the analytical formulation, since transverse shear effects are neglected, but such elements are known to be not recommended for torsional buckling applications.

Thin rectangular cross-section

The first case study concerns the beam with a rectangular cross-section. The dimensions retained for the cross-section are such that lateral-torsional buckling may be the preferred buckling mode for the considered loading conditions. Therefore, the width b is clearly smaller than the depth 2h, but nevertheless, the cross-section is chosen sufficiently thick in such a way that buckling occurs after the onset of plasticity (and before any other undesirable mode) for a maximum length range. The material parameters coincide with the properties of a S235 grade structural steel, except the hardening modulus which is set to a rather high value so as to prevent any localized extreme deformation during the shell finite element computations. The geometric and material parameters are listed in Tables 1 and2.

Half-depth h (mm) Width b (mm) 100 20 The beam is modeled using shell finite elements whose thickness is equal to the width b.

A quick mesh convergence study shows that 4 elements are needed along the depth direction whereas 30 elements are used along the beam length. A typical buckling mode is illustrated To reproduce the physical model, some loading and boundary conditions must be enforced at both ends. Owing to the use of shell elements, the centers of the two end sections are fixed in both transverse directions in order to meet the simply-supported boundary conditions whereas a single end section is fixed in the longitudinal direction to cancel beam rigid body modes. In addition, both end sections are fixed in the direction normal to the shell, in order to prevent torsional rotation at the beam ends. Last, a linear longitudinal shell edge load is properly applied to both ends so as to generate the equivalent of a bending moment.

Several numerical computations are performed considering beam lengths varying be-

J o u r n a l P r e -p r o o f

Journal Pre-proof tween 800 mm and 3200 mm. In each case, with such a cross-section geometry, the beam is shown to bifurcate to a lateral-torsional deformed shape while the beam has already partially plasticized. Table 3 shows the comparison between the critical bending moments obtained analytically with or without considering the hardening modulus. The corresponding positions (along the y-axis) of the transition between the elastic and plastic zones at the critical instant are also listed. In the case of perfect plasticity, the closed-form solutions ( 33) and ( 34) are employed, whereas in presence of hardening Equation ( 31) is solved numerically.

Length M cr (N.m) Transition (%) M cr (N.m) Transition (%) Error (%) The length is varied in such a way that buckling always occurs in the elastoplastic regime, but with a plastic zone kept to a minimum (when L = 3200 mm) towards a beam almost completely plasticized (when L = 800 mm). In all those conditions, the critical values are shown to be roughly the same, regardless of the fact that hardening is considered or not in the formulation (even if a high hardening modulus has been retained). The maximum relative error reaches 3% only. Accordingly, the transitions between elastic and plastic regions at the critical point are found to be very similar, whether or not hardening is taken into account.

L (mm) (E T = 0) c/h (E T = 0) (E T = 0) c/h (E T = 0) |M E T =0 cr -M E T =0 cr | M E T =0
It means that, in this configuration, the hardening modulus (as long as it is not excessively high) has no particular influence on the elastoplastic lateral-torsional buckling. It is certainly due to the fact that an elastic zone always remains in the vicinity of the beam neutral axis

J o u r n a l P r e -p r o o f

Journal Pre-proof (even in the case of a very short beam which is highly plasticized at the critical moment), which has a major influence on the value of the buckling load. In such practical conditions, it means that the closed-form solution for the elastoplastic lateral-torsional buckling moment is a very good estimation.

Next, the previous analytical values (considering E T = 0) are compared against the numerical results obtained using the shell finite element model (see Table 4). Again, the relative error between the critical bending moments does not exceed 4%. This outcome fully validates the closed-form expression.

Length M cr (N.m) M cr (N.m) Error (%) The same problem is finally solved using Abaqus and all the analytical and numerical solutions are plotted in Figure 5 Alternatively, the previous results can be displayed in a normalized graph plotting the dimensionless critical lateral-torsional buckling moment χ LT = M cr /M pl with respect to the which is equal here to:

L (mm) (E T = 0) (FE) |M E T =0 cr -M F E cr | M F E
M el = 2bh 2 σ y 3 (39)
• M pl is the maximum bending moment corresponding to the limit state of the beam in the case of perfect plasticity, which is equal here to:

M pl = bh 2 σ y = 3 2 M el (40)
In the case of a rectangular cross-section and considering the above definitions, the ratio between the elastic and plastic moments is always M el /M pl = 2/3. Accordingly, one can express the following relationship between χ LT and λ LT (when considering that the influence of hardening is negligible):

χ LT = 2 3λ 2 LT √ a (41)
where a only depends on β = 4/λ 4 LT (see Equations ( 33) and ( 34)) in elastoplasticity and is simply equal to 1 in elasticity.

I-beam

The second case of an I-beam is now considered. The dimensions of the reference crosssection are listed in Table 5. The web dimensions are similar to that of the previous rectangular cross-section whereas the flanges are a little less wider and thicker than the web.

With this configuration, several buckling mode types may arise (more than in the previous configuration). Their order of appearance highly depends on the beam geometry. The choice of a more massive web, in comparison with the flanges, permits to avoid any local or distorsional buckling mode before the sought lateral-torsional buckling phenomenon occurs.

The material parameters are the same as in the previous configuration (see Table 6), except the hardening modulus which is set to a higher value in order to prevent any localization of plastic deformations, especially in the flanges.

The beam is again modeled using shell finite elements with the appropriate thicknesses. needed for the web. For efficiency purposes, only 20 elements are finally used along the beam length as shown in Figure 9 where a typical buckling mode is illustrated. Several numerical computations are performed considering beam lengths varying between 800 mm and 1600 mm. For each case with such a cross-section, the beam is shown to bifurcate to a lateral-torsional deformed shape after the onset of plasticity (in the flanges

J o u r n a l P r e -p r o o f

Journal Pre-proof and a portion of the web). Table 7 and Figure 10 show the comparison between the analytical and numerical critical bending moments (the analytical solutions are derived from Equation [START_REF] Le Grognec | Elastoplastic bifurcation and collapse of axially loaded cylindrical shells[END_REF] including the effect of hardening). The corresponding positions (along the y-axis) of the transition between the elastic and plastic zones at the critical instant are also listed in with a hardening modulus ranging between 500 M P a and 5000 M P a. Within this range, it appears that the hardening modulus has only a limited influence on the value of the critical moment (it has been already checked that it has no influence on the critical moment of a beam with a rectangular cross-section). Nevertheless, the analytical and numerical critical moments follow the same trends, despite the differences between the two analyses. These last results complete the validation of the analytical solution for smaller hardening moduli. In contrast with the previous case, the critical bending moment of an elastic I-beam involves the warping inertia term. Instead of Equation ( 27), one must then use the following expression:

M e cr = π L GJEI y 1 + π 2 EI ω L 2 GJ ( 43
)
which turns out to be the solution of Equation ( 24), in elasticity, after inserting the buckling mode expression [START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF].

Figure 13 displays both elastic and elastoplastic (including hardening) analytical critical values in such a graph. Numerical finite element results are also added to the same graph for comparison purposes. The general trends observed before can be seen again. In contrast with the previous case, the elastoplastic critical moment exceeds here the plastic moment (with a ratio between the elastoplastic critical moment and the plastic moment of about 1.2). This result proves the more significant influence of hardening for the current configuration. In this case, one can observe a gap between the elastic and elastoplastic curves. The maximum value of λ LT considered in the previous calculations (around 0.6) corresponds to a relative transition coordinate of 94.3 % (see Table 7). For higher values of λ LT , namely between 0.6 and 1, the buckling phenomenon still occurs in the elastoplastic range, but the plastic zone at the critical time is likely to be reduced to only a portion of the flanges. This type of configuration is disregarded here since it has not been handled in the theoretical analysis.

J o u r n a l P r e -p r o o f

Journal Pre-proof (involving tangent moduli) was retained as it corresponds overall to a better option, which is likely to give solutions in line with numerical finite element results. Nevertheless, it could be interesting to replace the elastic shear modulus involved in the present formulation by an inelastic value (as recommended in [START_REF] Becque | The application of plastic flow theory to inelastic column buckling[END_REF], for example), which would require further developments but could perhaps help achieving better solutions with regard to experimental results. At the end, the resolution of the so-called bifurcation equation leads to differential governing equations which are first solved in general terms.

Then, two particular cases are investigated, namely the case of a thin rectangular crosssection and the case of an I-beam. In the first case, it is shown that the hardening modulus has almost no influence on the buckling values, at least when using geometric and material parameters relevant to practice. Under such conditions, the critical position of the transition between elastic and plastic zones along the thickness can be obtained as the positive root of a fourth-order polynomial equation, and a closed-form solution can thus be derived for the critical bending moment. In the second case, the solution is more generally deduced from the numerical resolution of a sixth-order polynomial equation, including the effect of hardening.

In both cases, the analytical solutions are validated by comparison against numerical finite element results obtained with an in-house shell finite element program. This finite element code can handle plasticity and finite transformations, and it is able to compute the elastic/plastic bifurcation loads during an incremental procedure thanks to branch-switching methods. Additional computations are performed using a beam finite element model on Abaqus software, for completeness purposes. Analytical solutions are shown to be in very good agreement with numerical results, as long as the same modes are naturally considered (the geometric and material parameters have been retained so that the lateral-torsional buckling happens in the first instance, before any other buckling phenomenon).

As a conclusion, analytical (and possibly closed-form) solutions for the critical bending moments of lateral-torsional buckling of elastoplastic beams are likely to be achieved and

Figure 1 :

 1 Figure 1: Lateral-torsional buckling of a simply-supported I-beam under pure bending

  r n a l P r e -p r o o f Journal Pre-proof

Figure 2 :

 2 Figure 2: Strain and stress profiles along the thickness of an elastoplastic beam under pure bending

Figure 3 :

 3 Figure 3: Geometry of the cross-section of the I-beam

4 .

 4 Young's modulus Poisson's ratio Hardening modulus Yield stress E (M P a) ν (-) H (M P a) σ y (M P a) In practice, an incremental calculation is performed and the buckling mode is evaluated by removing the trivial solution at the bifurcation point to the immediate postbuckling deformed shape (a suitable deformation amplification factor is applied for visibility purposes).

Figure 4 :

 4 Figure 4: Deformed modal shape of the beam with a thin rectangular cross-section obtained with the shell FE model

  for validation purposes. The commercial finite element code Abaqus does not offer the possibility to use branch-switching techniques and compute bifurcation values. Accordingly, only a maximum admissible moment can be estimated from an equilibrium curve such as the one plotted in Figure6. Despite this, and although an initial imperfection has been introduced which degenerates the problem and the response curve, the corresponding results are in good agreement with all the previous analytical and numerical solutions. Let us mention that the same transverse perturbation force of 200 N has been applied at mid-span in all cases. This probably explains the fact that the limitJ o u r n a l P r e -p r o o fJournal Pre-proof moment values are increasingly deviating downwards from the other buckling solutions as the beam is getting longer and thus more flexible. Another configuration resulting from the reduction of the thickness b to 15 mm has also been tested with Abaqus. With this thinner cross-section geometry, other failure modes are likely to occur in practice before the sought lateral-torsional buckling mode. The use of shell elements here would then lead probably to irrelevant results. Conversely, with the beam model in Abaqus, one still obtains consistent results with respect to the analytical solutions. In Figure7, only the analytical solutions without hardening are plotted together with the numerical results from Abaqus. The beam length is again varied within a certain range for which plasticity appears before the onset of buckling. The comments made for the previous configuration apply to the current one.These last results complete the validation of the analytical solution for this problem.

Figure 5 :

 5 Figure 5: Comparison between analytical and numerical critical moments of a beam with a rectangular cross-section -b = 20 mm

Figure 8

 8 Figure 8 displays both elastic (Equation (27)) and elastoplastic (Equation[START_REF] Lam | Arc-length method for passing limit points in structural calculation[END_REF] if no hardening is considered) analytical critical values in such a graph. Numerical finite element results are also added to the same graph for comparison purposes (results from Abaqus are not displayed in this figure since they do not really represent bifurcation values). It can be seen that for λ LT > 1 (for long beams), the lateral-torsional buckling naturally happens in the elastic range. Once λ LT < 1 (for short beams), the elastoplastic critical value is strictly lower than the elastic one. It is shown that the critical moment M cr barely exceeds the maximum bending moment M pl , namely χ LT is almost never higher than 1, since the hardening has no influence in this configuration.

Figure 8 :

 8 Figure 8: Normalized representation of the elastic/plastic lateral-torsional buckling moments of a beam with a rectangular cross-section

Figure 9 :

 9 Figure 9: Deformed modal shape of the I-beam obtained with the shell FE model

Figure 10 :

 10 Figure 10: Comparison between analytical and numerical critical moments of an I-beam: influence of the beam length

Figure 12 :

 12 Figure 12: Comparison between analytical and numerical critical moments of an I-beam: influence of the hardening modulus

Table 1 :

 1 Geometric parameters

Table 3 :

 3 Analytical critical bending moments of a beam with a rectangular cross-section

Table 4 :

 4 Analytical and numerical critical bending moments of a beam with a rectangular cross-section

  Two line of shell elements are employed to discretize the flanges whereas 6 elements are Length Web half-depth Web thickness Flange width Flange thickness

	L (mm)	h (mm)	b (mm)	w (mm)	t (mm)
	1000	100	20	80	12

Table 5 :

 5 Geometric parameters

	Young's modulus Poisson's ratio Hardening modulus Yield stress
	E (M P a)	ν (-)	H (M P a)	σ y (M P a)
	210000	0.3	5000	235

Table 6 :

 6 Material properties

Table 7 .

 7 Even if the critical values do not vary significantly, the transition at the critical instant moves from a highly plasticized configuration to a nearly entirely elastic web within the range of lengths considered. These results fully validate the analytical solution.

	Length M cr (N.m) Transition (%) M cr (N.m)	Error (%)
	L (mm) (analytical)	c/h	(FE)	|M analytical cr M F E -M F E cr | cr
	800	96781	23.29	96286	0.51
	1000	91732	41.04	91820	0.1
	1200	87423	60.04	87550	0.15
	1400	82911	78.14	83755	1.01
	1600	78256	94.3	79789	1.92

Table 7 :

 7 Influence of the beam length on the analytical and numerical critical bending moments of an I-beam

J o u r n a l P r e -p r o o f

Journal Pre-proof developed in this paper in a more comprehensive way (conversely, a beam finite element model would have validated the analytical expressions only, but not the general beam assumptions);

• it is provided not only with arc-length solution methods, but also with branch-switching algorithms, that enable one to deal with non-linear plastic bifurcation analyses, without resorting to any initial imperfections.

In the following numerical computations, the thin-walled beams are thus described by means of shell elements. As regards plasticity, the von Mises yield criterion with a linear isotropic hardening is adopted. A plane stress radial return mapping is employed to integrate the evolution equations and compute the plastic strains and the stresses, together with the reduced consistent elastoplastic tangent tensor. The shell surfaces are then discretized using eight-node isoparametric quadrilateral elements (involving quadratic shape functions) with reduced integration. The virtual work principle yields the discretized non-linear equations of the problem, which are solved by the iterative Newton-Raphson procedure, requiring the computation of the structural tangent stiffness matrix.

The computations are performed in an incremental way, and a quadratic arc-length method is used in order to proceed on the equilibrium branches and pass easily through possible load or displacement limit points found on the different paths [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF][START_REF] Crisfield | Non-linear finite element analysis of solids and structures[END_REF]. The standard solution scheme is modified according to Lam and Morley [START_REF] Lam | Arc-length method for passing limit points in structural calculation[END_REF], so as to cope with the complex roots which may be encountered when solving the quadratic equation in terms of the loading factor. It consists in eliminating as much as possible the component of the residual force which is orthogonal to the load direction, as often as necessary.

Branch-switching techniques are also included in the numerical procedure in order to detect the bifurcation points and bifurcate onto a given branch at any primary (or secondary) critical point. Specific methods are implemented following Riks [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF][START_REF] Riks | On formulations on path-following techniques for structural stability analysis[END_REF] and Seydel [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]. The

J o u r n a l P r e -p r o o f

Journal Pre-proof Then, the reference length (L = 1000 mm) is considered and the yield stress is varied in such a way that the beam still yields partially before buckling. Analytical solutions are compared against numerical results (see Table 8 and Figure 11). In both cases, the analytical and numerical critical bending moments almost coincide with each other (with a maximum relative error of about 2%).

Yield stress M cr (N.m) Transition (%) M cr (N.m) Error (%)

Conclusions

In this paper, an analytical study of the elastoplastic lateral-torsional buckling of perfectly straight beams has been carried out. A general formulation has been developed with the aim of obtaining analytical solutions for the critical moments of elastoplastic beams under pure bending and undergoing lateral-torsional buckling. This formulation is based on the 3D elastoplastic bifurcation theory, which has been already applied successfully to the case of elastoplastic Timoshenko beams under axial compression [START_REF] Le Grognec | On the plastic bifurcation and post-bifurcation of axially compressed beams[END_REF], two-layer beams with partial interaction [START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF][START_REF] Grognec | Plastic bifurcation analysis of a two-layer shear-deformable beam-column with partial interaction[END_REF], or sandwich beams [START_REF] Douville | Exact analytical solutions for the local and global buckling of sandwich beam-columns under various loadings[END_REF][START_REF] Le Grognec | Elastoplastic buckling and post-buckling analysis of sandwich columns[END_REF] under axial compression/pure bending, but also to the lateral buckling of offshore pipelines including twisting effects [START_REF] Grognec | Investigation of the torsional effects on the lateral buckling[END_REF].

In the present study, it is the first time that the formulation is used in the elastoplastic regime considering bending loading conditions, in such a way that, at critical instant, the beam is partially plasticized with a priori unknown elastic and plastic zones. The elastoplastic buckling problem is solved in the framework of the J 2 flow theory, assuming the von Mises yield criterion and a linear isotropic hardening. Let us mention that the flow theory

J o u r n a l P r e -p r o o f

Journal Pre-proof to give satisfactory results for any cross-sectional shape. These new elastoplastic solutions deviate significantly from the classical well-known elastic values as soon as plasticity reaches a sufficiently high level, and are thus of great importance for relatively short and thick beams.

Appendix A. Solving a quartic equation

Let us consider a quartic equation which takes the standard form:

where a, b, c and d are real. Equation (A.1) can be rewritten as follows:

J o u r n a l P r e -p r o o f

Journal Pre-proof By adding x 2 + ax 2 y + y 2 4 to both sides of the above equation, one obtains:

The objective is then to choose a value for y such that the right hand side of Equation (A.3) becomes a perfect square. This can be done by letting the discriminant of the quadratic function in x become zero, namely:

This leads to a cubic equation with respect to y:

It is well-known that a cubic equation can be analytically solved and has always at least one real root, say y 0 . Consequently, Equation (A.3) can be rewritten as:

Finally, the four roots (real and/or complex) of the quartic equation can be found by solving the two following quadratic equations for x: ☒ The authors declare that they have no known competng fnancial interests or personal relatonships that could have appeared to infuence the work reported in this paper.

☐The authors declare the following fnancial interestsppersonal relatonships which may be considered as potental competng interests: