Grégoire Corre

| Arthur Lebée

| Karam Sab

| Mohammed Khalil Ferradi

| Xavier Cespedes

Arthur Lebée
email: arthur.lebee@enpc.fr

The Asymptotic Expansion Load Decomposition elasto-plastic beam model

Keywords: Model reduction, Higher-order beam model, Asymptotic expansion, Plasticity, Structural analysis

A new higher-order elasto-plastic beam model is derived and implemented in this paper. The reduced kinematic approximation is based on a higher-order elastic beam model using the asymptotic expansion method. This model introduces new degrees of freedom associated to arbitrary loads as well as eigenstrains applied to the beam.

In order to capture the effect of plasticity on the structure, the present elasto-plastic model considers the plastic strain as an eigenstrain imposed on the structure and new degrees of freedom are added on the fly into the kinematics during the incrementaliterative process. The radial return algorithm of 2 plastic flow is used. Because of the constant evolution of the beam kinematics, the Newton-Raphson algorithm for satisfying the global equilibrium is modified. An application to a cantilever beam loaded at its free extremity is presented and compared to a 3D reference solution.

The beam model shows satisfying results even at a local scale and for a computation time significantly reduced.

INTRODUCTION

Beam models combine a significant simplicity in the modeling of structures and a strong time efficiency in the computation of the results. Hence, the development of linear elastic beam models has been widely investigated. Since the earliest model of Euler-Bernoulli, a large number of elastic beam models has been suggested. The main issue in elasticity lies in the definition of the kinematics of the model which entirely determines its efficiency. Higher-order beam models offer extended kinematics able to describe more accurately local phenomena in beam structures. There are numerous way to build such kinematics assuming a priori a variable separation between the longitudinal coordinate of the beam and the cross-sectional coordinates. For instance, Proper Orthogonal Decomposition, reduced basis approaches [START_REF] Miled | A priori hyper-reduction method for coupled viscoelastic-viscoplastic composites[END_REF] or Proper Generalized Decomposition 2,3,4 may be used depending on the existence of a priori error estimators and the wish to update the kinematics during the computation.

It turns out that the formal asymptotic expansion of the 3D beam problem with respect to the inverse of the slenderness of the beam provides such a basis which may be derived a priori for any given beam cross-section. This approach, suggested early [START_REF] Bernadette | A Galerkin spectral approximation in linearized beam theory[END_REF] , was recently implemented in the case of linear elastic beams submitted to arbitrary loads as well as eigenstrains [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF][START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] . Two noticeable observations were made by Miara and Trabucho [START_REF] Bernadette | A Galerkin spectral approximation in linearized beam theory[END_REF] . First, the formal asymptotic expansion delivers a free family of kinematic enrichment which is dense in the space of the 3D solution. This means that going sufficiently high in the expansion allows arbitrary refinement of the 3D solution. Second, the truncation of this family ensures that the corresponding beam model is asymptotically consistent except at the boundary. This means that the kinematic enrichment delivered by the formal asymptotic expansion is optimal in terms of approximation error far from the extremities of the beam.

Introducing elasto-plastic behavior is more complex. The inherent non-linearity of plasticity and the incremental nature of plastic analysis makes the definition of a relevant kinematics more difficult. Two main approaches are followed when solving an elasto-plastic beam problem: 1D elasto-plastic beam model based on a priori cross-section analysis and 3D elasto-plastic beam models based on a 3D beam kinematics.

The first natural approach is to express the plastic flow in terms of generalized beam variables and to solve an elasto-plastic 1D problem. This requires the elasto-plastic analysis of the cross-section for pure or combined generalized stresses and the derivation of the corresponding yield surface. The cross-section analysis may be incremental or based on limit analysis but assumes a uniform distribution of generalized stresses in the longitudinal direction: normal force, shear forces, bending moments and torque. In this direction, closed-form solutions were first devised and numerical approximation of cross-section analysis were implemented later. Indeed, the elastic problem of pure-torsion was early solved by Saint-Venant and the plastic analysis of the torsion of a beam was sketched. Nadai [START_REF] Nadai | Plasticity: A mechanics of the plastic state of matter[END_REF] was the first to suggest a solution for the elasto-plastic problem and to calculate a plastic torque thanks to the sand-heap analogy. Then, closed-form solutions of the plastic torque have then been developed for the torsion of beams with common cross-sections: Christopherson [START_REF] Christopherson | A theoretical investigation of plastic torsion in an I-beam[END_REF] solved the torsion of I-beams, Sokolovsky [START_REF] Sokolovsky | Theory of plasticity[END_REF] obtained a solution for beams with oval sections and Smith and Sidebottom [START_REF] Smith | Inelastic behaviour of load-carrying members[END_REF] for prismatic bars with rectangular sections. Closed form solutions have also been obtained for bending analysis. Combined generalized stress state were also investigated [START_REF] Hill | On the plastic distortion of solid bars by combined bending and twisting[END_REF][START_REF] Boulton | Plastic twisting and bending of an I-beam in which the warp is restricted[END_REF] . A key difficulty is the derivation of a yield surface directly function of the beam generalized stress taking into account correctly their possible interactions as well as hardening. There were recent improvements in this direction, approximating the yield surface with facets or ellipsoids [START_REF] Bleyer | Yield surface approximation for lower and upper bound yield design of 3D composite frame structures[END_REF] . Once the yield surface is defined, there remains to compute the elasto-plastic response of the beam, either with closed form solutions [START_REF] Boris | Analytical solutions in elasto-plastic bending of beams with rectangular cross section[END_REF] , limit analysis [START_REF] Olsen | Rigid plastic analysis of plane frame structures[END_REF] or by means of finite element approximations [START_REF] Argyris | Finite element analysis of two-and three-dimensional elasto-plastic frames-the natural approach[END_REF][START_REF] Papadrakakis | A computationally efficient method for the limit elasto plastic analysis of space frames[END_REF] . This approach has the advantage to present fast computation time, since only a 1D elasto-plastic problem needs to be solved. However, its accuracy remains limited by the beam theory assumptions. First, it cannot handle local phenomena related to the distribution of the applied load as well as to the boundary conditions. Second, it provides only an averaged description of the actual stress in the beam.

In order to improve the accuracy of the beam model, the second classical approach consists in setting a beam kinematics expressing the 3D-displacement field in a separate form between the cross-sectional coordinates and the longitudinal coordinate. This kinematics may be defined a priori or may evolve during the incremental procedure. For a fixed increment of the generalized displacements, the corresponding 3D stress is computed and the yield criterion is expressed locally. A local algorithm such as the radial return is processed on the whole body to compute the local plastic state of the beam. This locally admissible stress state is integrated on each cross-section yielding the corresponding longitudinal distribution of the beam generalized stresses. Finally, the beam global equilibrium is ensured with a standard Newton-Raphson procedure. This approach was compared with purely 1D approach by Gendy and Saleeb [START_REF] Gendy | Generalized yield surface representations in the elasto-plastic three-dimensional analysis of frames[END_REF] . The 3D approach appeared to be much more accurate and closer to the full 3D solution for a reasonably higher computation time than a 1D approach. Many numerical models therefore adopt the continuumbased description in terms of 3D stress components in order to benefit from its accuracy. The main difficulty lies again in the definition of a relevant kinematics able to describe the displacement related to plastic flow.

Most approaches where the kinematics is fixed a priori rely on the ones already used in linear elasticity such as Euler-Bernoulli, Timoshenko kinematics or even Saint-Venant solution, eventually with non-linear geometric corrections. For instance, Bathe and Chaudhary 20 suggested to introduce the Saint-Venant warping function into the kinematics in order to compute the elastoplastic torsion of a rectangular beam. Once the kinematics is defined, there remains to choose the number of integration points in the cross-section in order to compute precisely the local plastic flow. Multiplying integration points improves the accuracy of the results at the price of a higher computation time of the cross-sections integrals. This is the spirit of multi-fiber beam models (see for instance [START_REF] Mazars | Using multifiber beams to account for shear and torsion. Applications to concrete structural elements[END_REF] ).

Another direction is to enrich arbitrarily the section kinematics with degrees of freedom not necessarily related to classical cross-section displacements. An early attempt was made by Bathe and Wiener 22 who performed the elastic-plastic analysis of I-beams in bending and torsion composed of three simple beam elements. This concept was formalized extensively by Carrera et al. [START_REF] Erasmo | Beam Structures[END_REF] and co-workers.

Because plastic flow may not be easily known a priori a natural improvement of the preceding methods is to update the beam kinematics during the load increments. This is the direction followed by Baba and Kajita 24 who suggested a method in which a warping mode is determined according to the plastic state of each cross-section and which was recently updated by Tsiatas and Babouskos [START_REF] George | Elastic-plastic analysis of functionally graded bars under torsional loading[END_REF] . However, in this approach, it is necessary to compute a 2D elasto-plastic cross-section problem, which remains computationally costly.

In this paper, the linear higher-order beam model based on the formal asymptotic expansion [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF][START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] is extended to elasto-plasticity in the small strains framework. This is achieved as follows. First, a higher-order kinematics is computed a priori for the considered section and applied load. Second, during the incremental procedure, this basis is updated with few displacement modes related to the plastic flow which occurs in the beam. More precisely, the plastic strain in some chosen cross-sections is considered as an eigenstrain load and used for computing the corresponding section displacement following the formal asymptotic expansion derived by Corre [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] . This approach presents two major advantages. First, it does not require additional elasto-plastic computations in the cross-section. Second, the number of beam degrees of freedom remains very limited (about 20) thanks to the sparsity of the kinematics. Indeed, the kinematics is enriched with very few degrees of freedom (up to 10) related to the plastic flow observed during the computation. From the optimality result proved in [START_REF] Bernadette | A Galerkin spectral approximation in linearized beam theory[END_REF] , this approach is expected to be more efficient than arbitrary kinematic refinements. Note that, contrary to Nonuniform Transformation Field Analysis [START_REF] Michel | Nonuniform transformation field analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the nonuniform transformation field analysis[END_REF][START_REF] Roussette | Nonuniform transformation field analysis of elastic-viscoplastic composites[END_REF][START_REF] Felix | Reduced basis hybrid computational homogenization based on a mixed incremental formulation[END_REF][START_REF] Felix | GPU accelerated computational homogenization based on a variational approach in a reduced basis framework[END_REF] where a basis of plastic strains is introduced with the corresponding plastic multipliers, in the present approach, displacement plastic modes are added to the total 3D displacement approximation and plasticity is processed at each integration point of the 3D body.

The paper is organized as follows. The formulation of the higher-order elasto-plastic beam model is first presented in Section 2: the definition of the kinematics thanks to the asymptotic expansion method is briefly recalled. The adaptation of this higherorder beam model to the framework of plasticity is then presented. Section 3 is dedicated to the numerical discretization and the description of the iterative-incremental plasticity algorithm. A radial return algorithm is used locally and an adaptation of the Newton-Raphson procedure is suggested to satisfy the global equilibrium. An application of the model to a cantilever beam is conducted is Section 4 and the influence of some parameters is investigated.

THE ELASTO-PLASTIC BEAM MODEL

The elasto-plastic boundary value problem

FIGURE 1 The beam configuration

We consider a beam occupying the prismatic domain Ω (Figure 1 ) with a length and a cross-sectional typical size ℎ. The boundary Ω is the union of the lateral surface Ω and the two end sections  ± (clamped). The longitudinal coordinate is 3 and the section coordinates are 1 and 2 denoted as1 , the corresponding reference frame is denoted ( , 1 , 2 , 3 ) where is an arbitrary point of the plane 3 = 0.

The constitutive material of the beam is only function of the section coordinates and invariant in the longitudinal direction. Without limitation, the fourth order elastic stiffness tensor ( ) is assumed isotropic.

Let [0, ] ⊂ ℝ + be the time interval of interest of the problem. The displacement of the beam is defined by the function

∶ Ω × [0, ] → ℝ 3 (1)
and the total linearized strain tensor is the symmetric gradient of . The total strain splits into an elastic part e and a plastic part p :

( , ) = e ( , ) + p ( , ) (2)
We consider an external body force ( , ) defined on Ω×]0, ], and a surface traction ( , ) defined on Ω ×]0, ] loading the beam. The evolution is elasto-plastic, quasi-static and under small deformation. The corresponding 3D elasto-plastic boundary value problem writes as:

div + = 0 on Ω = ∶ ( -p ) on Ω = ∇ on Ω ⋅ = on Ω = 0 on  ± ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ × [0, ] , (3)
where is the outer normal to Ω , ∇ is the symmetric part of the 3D gradient operator div is the 3D divergence operator and the stress tensor. The flow rule is not recalled in equations (3) but is detailed in Appendix A.1. The 3D beam is clamped at both extremities. Other boundary conditions may be applied, depending on the approximation of the total displacement, and are detailed in the following section.

The common way to solve an elastoplastic problem is to use an incremental-iterative procedure. For each load increment, a local algorithm ensures that the local stress satisfies the elasto-plastic constitutive law. A global algorithm ensures that the body is globally in a statically admissible state. The standard algorithms used for this procedure are recalled in Appendix A.

The higher-order beam model

Solving the boundary value problem (3) with a 3D mesh of the body implies a large number of elements when the beam becomes slender and the cross-sectional resolution must be preserved. In addition to poor numerical conditioning, the assembly of the stiffness matrix and the resolution of the global balance equation quickly become time consuming. For a beam, the longitudinal dimension is larger than the two other dimensions, characterized by the length ℎ (Figure 1 ). Therefore, a more time-efficient model may be obtained, taking this geometrical feature into account.

Assuming a separation of the cross-sectional coordinates ( 1 , 2 ) and the longitudinal coordinate 3 , a dimensional reduction of a 3D model into a beam model may be obtained from the following expression of the displacement:

= mod ∑ =1 ( ) ( 3 ) (4)
where are mod 3D-displacement modes defined on the 2D cross-section of the beam and are the beam generalized displacements. The 3D-displacement is therefore decomposed into a "2D+1D" form.

Let us recall that the principle of virtual work of boundary value problem (3) writes as:

( , ̂ ) = ∫ Ω [ ] ∶ [ ̂ ]dΩ -∫ Ω ⋅ ̂ -∫ Ω ⋅ ̂ d = 0, ∀ ̂ ∈ ( ± , ). (5)
where, ( ± , ) is the set of 3D-displacement fields vanishing on  ± . Hence, restricting the displacements and ̂ to the approximation (4) into equation ( 5) and integrating the constitutive law for each cross-section leads to an expression of  in terms of the 1D kinematic unknowns and ̂ . The corresponding 1D global equilibrium equation may be solved and the total displacement can then be reconstructed thanks to equation (4).

Whereas, the dimensional reduction described previously is fairly classical in beam modeling, the choice of the displacement modes is critical for the efficiency of the final algorithm. Indeed, the kinematics can be arbitrarily assumed. It can also be established by means of an automated procedure, for any fixed cross-section, which is detailed in the following.

The Asymptotic Expansion Load Decomposition beam model

The higher-order elastic beam model developed in Ferradi et al. [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] is based on the asymptotic expansion method which offers a systematic procedure for enriching the kinematics of the beam in the form of equation ( 4) in the framework of linear elasticity. In 6 , the kinematics is composed of two kinds of modes. The first collection of modes are the 12 modes of the Saint-Venant's solution. The collection of modes they form is denoted by S-V . These modes are specific to the geometry of the cross-section. The second collection of modes comprises modes both specific to the geometry of the cross section and to the applied load. For a given applied force ( and in ( 3)), the model presented by Ferradi et al. enriches the kinematics with specific additional modes. This additional basis of force modes is denoted by .

Using the developments made in Ferradi et al. [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] , the computation of the modes of the Saint-Venant's solution and modes specific to the force applied on the structure are briefly sketched in this section. The boundary value problem (3) is considered with p = 0 and is therefore linear.

Scaling and expansion

Noticing that the ratio = ℎ is small for the geometry of a beam, a change of coordinates is operated as follows:

1 , 2 , 3 = ℎ 1 , ℎ 2 , 3 (6)
The scaled section is denoted  0 , and  0 is its boundary. Then, the method is based on two main assumptions. First, the load applied on the structure is assumed as products of a single longitudinal function and cross-sectional distributions:

= 1 2 ̃ ( ) ( 3 ), 3 = 1 ̃ 3 ( ) ( 3 ), = 3 ̃ ( ) ( 3 ), 3 = 2 ̃ 3 ( ) ( 3 ). (7)
Second, the variables are expressed as power series of the scaling ratio (asymptotic expansion):

= 0 + 1 + 2 2 + ... , = 0 + 1 + 2 2 + ..., = 0 + 1 + 2 2 + ..., (8)
and introduced in (3). The powers of are then identified: for each power ∈ ℕ, each compatibility equations, boundary conditions and constitutive equations for and equilibrium equations for -1 yield an auxiliary problem on the cross-section which splits in two uncoupled 2D boundary value problems.

Transverse displacement

First, the in-section displacement problems (transverse mode)  are gathered for ≥ 0:

 +1 ∶ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ , + -1 3,3 + 3 ̃ = 0 on  0 = 2 + , 33 = 2 33 + on  0 2 = +1 , + +1 , , 33 = 3,3 on  0 = 3 ̃ on  0 (9)
where -1 = 0 and 3 = 1 if = 3 and 3 = 0 else. For a simply connected cross-section, this 2D boundary value problem on the displacement +1 is well-posed if the applied load is globally self-equilibrating for translations and a rotation in the plane of the cross-section:

∫  0 -1 3,3 + 3 ̃ d + ∫  0 3 ̃ d = 0 and ∫  0 -1 3,3 + 3 ̃ d + ∫  0 3 ̃ d = 0 (10)
where 11 = 22 = 0, 12 = +1, 21 = -1. The solution is thus defined up to a rigid motion of the section in its plane.

Longitudinal displacement

Second, the longitudinal displacement problems (warping mode)  are obtained for ≥ 0:

 +1 ∶ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 3 , + -1 33,3 + 2 ̃ 3 = 0 on  0 3 = 2 3 + 3 on  0 2 3 = +1 3, + ,3 on  0 3 = 2 ̃ 3 on  0 (11)
Again, for a simply connected cross-section, this 2D boundary value problem on the displacement +1 3 is well-posed if the load applied is globally self-equilibrating for the longitudinal translation:

∫  0 -1 33,3 + 2 ̃ 3 d + ∫  0 2 ̃ 3 d = 0 ( 12
)
The solution is defined up to a longitudinal displacement.

Kinematic approximation

The successive resolutions of  +1 and  +1 for each ≥ 0 yields a collection of displacement modes 0≤ . In problems  +1 and  +1 , the computation of the modes specific to the force applied is linearly dependent on the longitudinal function and its higher gradients in the longitudinal direction. As a result, problems ( 9) and ( 11) may be solved for a unit on a single cross-section and only the cross-sectional functions ̃ , ̃ 3 , ̃ and ̃ 3 need to be known. Once orthonormalized thanks to a Grahm-Schmidt procedure,with respect to the 2 -norm defined as follows:

‖ ‖ ‖ ‖ = ⎛ ⎜ ⎜ ⎝ ∫ 3 ∑ =1 2 d ⎞ ⎟ ⎟ ⎠ 1∕2 (13)
these modes yield the kinematic approximation of the model of the form (4). For a fixed load which does not generate a boundary layers, Miara and Trabucho 5 proved a higher-order convergence result with this kinematic approximation.

Boundary conditions

In the 3D boundary value problem (3), the beam extremities were assumed fully clamped. Considering the approximated kinematics (4), this is achieved enforcing = 0 at extremities. Other boundary conditions may be applied. Indeed, the first six Saint-Venant modes corresponds to the rigid motion of the section. Restraining only these degrees of freedom is actually the boundary condition classically used in structural mechanics: warping and transverse displacements are let free. Note that, the computation of the modes is only made possible by assuming variable separation of the applied load: it is decomposed and expressed as the product of a function of the cross-sectional coordinates and a function of the longitudinal coordinate. Therefore this model is here named the Asymptotic Expansion Load Decomposition beam model (AELD-beam model).

The AELD extended to the case of eigenstrains

This higher-order beam model has been extended in Corre et al. [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] to the case of eigenstrains. Indeed, this enables the model to deal with various situations such as creep, thermal loads or prestressed loads. Since the total plastic strain may be considered as an eigenstrain, it appears relevant to enrich the total kinematics with the corresponding displacement modes. Considering a prescribed and fixed eigenstrain p , expressed as the product of functions of the cross-sectional coordinates and three functions of the longitudinal coordinate:

p 33 = ̃ 33 ( ) 1 ( 3 ), p 3 = ̃ 3 ( ) 2 ( 3 ) and p = ̃ ( ) 3 ( 3 ), (14)
new modes are computed thanks to the same systematic procedure as before -fully detailed in [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] -and are added to the kinematics of the model. The basis of modes specific to p is denoted p . The distinction between 33 , 3 and is motivated by the fact that 33 and are related to traction and bending whereas 3 is related to torsion at leading order in the asymptotic expansion. As for the modes of , the modes specific to the distribution of the plastic strain in a cross-section are linearly dependent on the longitudinal functions .

Adaptation of the AELD-beam model to the elasto-plastic behavior

The AELD-beam model introduced in the previous section has proven its efficiency for linear elastic materials. The model is now adapted to the elasto-plastic behavior.

We consider the elasto-plastic boundary value problem expressed in equation (3). The first collection of modes to take into account is the basis of Saint Venant modes -described previously. The load applied on the beam then generates an additional collection of force modes .

The introduction of degrees of freedom related to the plastic strain in the beam is necessary to correctly describe the effect of plasticity in the total displacement approximation. Therefore, the plastic strain computed at a given iteration of the global algorithm is taken into account for enriching the kinematics of the following iteration. This is possible using the procedure described in the previous section for a fixed p distribution in the cross-section. Considering now the whole beam, p is not longitudinally uniform. Hence, several chosen cross-sections can be used for taking snapshots of the plastic strain in order to sufficiently enrich the kinematics of the model. These new plastic modes are computed and added to the kinematics on the fly. The basis of modes specific to a plastic strain p is denoted by p . Finally, the kinematics of the model is evolving during the Newton-Raphson procedure and is the union of the basis S-V , and p . This union of basis is orthonormalized to form the total basis, and the total number of modes mod used in the kinematics is the number of modes in the total basis.

THE ELASTO-PLASTIC ALGORITHM

The implementation of the general framework introduced in the previous section is now detailed. This requires first the definition of the numerical approximation of the 3D body. Then, the incremental resolution of the elasto-plastic problem is adapted so that processing the local constitutive equations remains standard whereas the global equilibrium iterations are performed with the reduced basis. Hence, the local algorithm remains defined by Algorithm 2. Major changes are made at the global level of the algorithm.

Numerical approximation of the higher order beam model

The approach suggested in the previous section requires the definition of a 3D mesh of the beam composed of cross-sections meshes positioned along the longitudinal direction (Figure 2 ). Indeed, these cross-sections will be the domain of integration of the constitutive law in the principle of virtual work (5).

FIGURE 2

Discretization of a square beam

Longitudinal discretization

A longitudinal discretization of the beam is defined for the functions ( 3 ) introduced in equation ( 4). As in [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] , we choose the same collection of NURBS basis functions for each generalized displacement ( 3 ):

( 3 ) = NURBS ∑ =1 ( 3 ) , ( 15
)
where ( 3 ) are the NURBS interpolation function and , are the corresponding degrees of freedom. Note that, contrary to conventional finite element interpolation, , is not the displacement at a given node, except at the extremities: = 1 or = NURBS . The number of NURBS interpolation functions is NURBS = knot + order -1, where knot is the number of knots used for the definition of the NURBS and order is their interpolation order. A set of longitudinal integration points are also defined for the integration of the interpolation functions. It is natural to place the cross-section meshes at the positions of these longitudinal integration points. This set of > NURBS longitudinal positions is denoted 1 , ..., .

Cross-section discretization

The cross-section mesh used for the computation of the modes is the same as the one used in Corre et al. [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] : the modes are computed by means of quadratic Lagrange triangle elements:

( ) = sec ∑ =1 ( ) , ( 16
)
where sec is the number of nodes in the section, ( ) are Lagrange interpolation functions and , are the nodal values of the displacement. The number of elements in the section is denoted and the number of Gauss points by element is denoted .

= × is therefore the total number of Gauss points in the section.

As already mentioned, the values of the incremental fields Δ , Δ p and Δ are computed in the cross-sections meshes. For simplicity, the discretization of each cross-section is the same as the one defined for the computation of the modes and the state variables are evaluated at the Gauss points of the quadratic triangle elements.

Plastic-mode cross-section

During the computation, sections where the incremental plastic strain is not zero are gathered in ℙ :

ℙ = ∈ 1, ..., ∕∃ ∈ 1, ..., , Δ p ( , ) ≠ 0 (17)
where denotes the transverse coordinates of the th Gauss point of the section. All cross-sections in ℙ could be used for the computation of the plastic modes: for each plastic strain distribution in each cross section, one or several modes could be computed. However, it would excessively increase the number of generalized displacement degrees of freedom dof and also increase the computation time dedicated to the corresponding modes. In order to limit the number of plastic modes to a few, only one cross-section called plastic-mode cross-section is chosen for taking snapshots of the plastic strain distribution. As a first approach, this choice is based on an educated guess and will be automatized in the future. Alternative approaches based on some projection of the collection ℙ of plastic strains could also be considered.

Adaptation of the Newton-Raphson procedure

Formulation of the tangent stiffness of the beam model

For a given basis of mod displacement modes, the numerical approximation of the total displacement may be written as follows:

( ) = mod ∑ =1 NURBS ∑ =1 sec ∑ =1 ( ) ( 3 ) , , = NURBS ∑ =1 sec ∑ =1 ( ) ( 3 ) , ( 18
)
where the 3D = NURBS × sec local displacement degrees of freedom are:

, = mod ∑ =1 , , . (19)
Considering the second form of ( 18), it appears that the 3D strain ( ) and consequently the local constitutive equation may be directly computed from , without the need to specify the basis . Hence, provided the finite element solution is stored as the collection of the , , the local integration of the constitutive equations as well as the computation of the local elasto-plastic tangent stiffness ep remain unchanged. Hence this framework appears as a generalization of multi-fiber models to any kind of basis . Accordingly, { } denotes the finite element vector of , . { } and {Δ } denotes the corresponding iterations and increments.

For a fixed basis , injecting the numerical approximation of the kinematics (18) into the principle of virtual work (5) leads to the expression of the residual expressed in terms of the increment of the dof = mod × NURBS generalized displacement degrees of freedom Δ , and the corresponding test degrees of freedom. The standard Newton-Raphson procedure is used in order to cancel this residual (see Appendix A.2.2) which leads to the following reduced equilibrium equation:

ep { } = { } , ( 20
)
where ep ∈ (ℝ dof ) 2 is the global stiffness matrix of the beam model expressed in basis , computed with the local elastoplastic tangent stiffness ep +1 , { } is the finite element vector of the generalized degrees of freedom , and { } ∈ ℝ dof is the residual vector. Because, the tangent stiffness as well as the residual both depend on the choice of the basis , they need to be updated each time the basis is changed. Significant computational time is gained because dof which sets the size of the tangent stiffness is much smaller than the rather large 3D which is required for a sufficiently detailed description of the fields in the cross-section. In practice about mode = 20 modes are used whereas in approaches only based on variable separation such as the Carrera Unified Formulation mode = sec since the basis functions are directly the interpolation functions of the cross-section mesh.

Description of the global algorithm

The modified global algorithm is presented in Algorithm 1 and corresponds to the following procedure.

The basis of modes is first initialized as described in [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] and is composed of the 12 modes of the Saint-Venant solution S-V and fAE modes associated to the applied load . This collection of modes is then orthonormalized:

0 = S-V ∪ ⟂ (21)
where (•) ⟂ means that the basis is orthonormalized. During the orthormalization procedure some modes may appear as redundant. This is the case if several load distributions in the cross-section as well as eigenstrain distributions are mixed. Hence, redundant modes are discarded with a 10 -3 tolerance. While increments do not generate plastic flow, the basis remains unchanged. The global tangent stiffness corresponds to the elastic one and each increment is solved in one iteration.

Let assume that increment + 1 generates a non-vanishing plastic strain increment Δ p +1 at the first iteration = 1. In this case, before starting the second iteration, the basis of modes is enriched and orthonormalized with pAE modes computed from the plastic strain distribution Δ p +1 observed in the plastic-mode cross-section. This also requires the update of the residual. It has been noticed from experience that plastic modes computed at subsequent iterations of the increment were very similar. Therefore the basis used at iteration = 2 is kept until the convergence of the increment is reached. However, the converged plastic strain of the increment Δ p +1 may have changed. Hence, at the first iteration = 1 of the following increment + 2, the basis is updated, replacing only plastic modes with new ones. Again, at the second iteration = 2 the basis is updated and then remains fixed until the convergence of the increment. This choice of updating the plastic modes only at the first two iterations of the increment remains valid as long as the load increments are not too important.

APPLICATION TO A CANTILEVER BEAM

Study of alternative Newton-Raphson methods

Before exposing the performance of the beam model on a I-beam, we first investigated on the possible alternatives to the standard Newton-Raphson method. The most time-consuming step of Algorithm 1 is the assembly of the consistent elasto-plastic stiffness matrix. It is assembled by an integration operated both cross-sectionally and longitudinally. The update of the basis of modes and the update of the consistent elasto-plastic moduli ep imply a new computation of the stiffness matrix.

A simplified method commonly used in standard 3D plasticity consists in approximating at each iteration the consistent elastoplastic stiffness matrix by the elastic stiffness matrix. This method, called the modified Newton-Raphson method [START_REF] Zienkiewicz | The finite element method[END_REF] , naturally implies more iterations within a load increment, but each iteration is computed faster since it avoids the update of the stiffness matrix. For the present higher-order beam model, it means that we always consider the elastic moduli instead of the consistent elasto-plastic moduli ep for the assembly of the global stiffness. However the stiffness matrix must still be computed each time the basis of modes changes.

A third solution can be formulated in between the Newton-Raphson and the modified Newton-Raphson method. The elastic moduli is updated at the first iteration of each increment, but is kept constant during the whole increment. The update is therefore operated only once. This method, called the quasi Newton-Rapshon's method, should provide time performances in between the performances of two first methods.

The three methods were used on a simple cantilever beam loaded at its end in order to assess their respective performance. As expected, the most time-efficient method is the standard Newton-Raphson method, followed by the quasi Newton-Raphson method and the modified Newton-Raphson method. A higher average time per iteration is observed for the Newton-Raphson method, but the gain in terms of convergence is sufficiently important to yield the shortest computational time. The standard Newton-Raphson method will therefore be the method considered in this paper.

Cantilever beam loaded at its free extremity

To illustrate the efficiency of the model presented, we consider a steel beam clamped at one end and loaded on its free end. The beam chosen is a wide flange beam HE600M. This section is class 1 in Eurocode 3, meaning that the beam reaches its limit of elasticity with no risk of local buckling. The geometry of the 6 m long beam is detailed in Figure 3 . A load is applied with eccentricity at the top edge of the free end of the beam, as represented in Figures 3 and4 end if

16: if = 1 then 17: ref = ‖ ‖ { } ‖ ‖ 18:
end if

19:

Assemble the consistent elasto-plastic stiffness matrix We consider the following values for the Young's modulus the Poisson's ratio, the strain hardening modulus and the yield modulus:

= 210 GPa, = 0.3, = 0.02 , 0 = 235 MPa (22) All the computations are performed on a processor i7-4510U (2 cores at 2.00 GHz).

Higher-order beam solution 0

The model is first computed with a set of parameters chosen with an educated guess. This solution is called 0 . Some sensitivity studies are carried out later in the following sections.

The section of the solution 0 is meshed with 399 quadratic triangle Lagrange elements, as shown in Figure 3 .

The NURBS basis functions in the longitudinal direction are defined by the following knot vector: NURBS = {0, 0.125, 0.25, 0.5, 1, 2, 3, 4, 5, 6} and knot = 10. We consider second-order NURBS: NURBS = 11. Using Simpson's integration, the total number of integration points is defined by the relation:

= 1 + 2 × ⌊ order +1 2 ⌋ × knot -1 (23)
where ⌊ ⌋ denotes the integer value of . The = 19 integration points of longitudinal axis are the knots NURBS where an integration point is introduced in the middle of each interval. The interpolation functions are represented in Figure 5 .

The corresponding integration sections are represented in Figure 6 . The mesh is refined close to the clamped extremity since plasticity is expected to occur mainly at this location.

FIGURE 5

Second-order NURBS basis functions used for the longitudinal interpolation of the element

FIGURE 6 mesh composed of 19 integration sections

The plastic-strain cross-section is placed at 3 = 0.25 m at the 5 th integration point. The number of force modes is fAE = 4 , and of plastic modes is pAE = 9 (this choice is based on experience).

During the computation, the maximum number of modes in the basis is 22. Indeed, in the orthormalization procedure, redundant modes are discarded. The number of interpolation shape functions being 11, the maximum number of degrees of freedom during the computation is therefore dof = 242. This number could be reduced by associating the plastic modes only to the interpolation functions with non-zero values where plasticity has been detected. But at this time, the modes of the basis are considered all along the beam element.

3D Reference solution

A 3D reference solution is computed with the finite element software Code_Aster. The beam is meshed by extruding a crosssection with 430 triangles along the longitudinal axis. The longitudinal discretization is the same as for the beam model and forms 7740 prismatic elements in total. The prismatic elements are interpolated with quadratic functions.

Comparison with the 3D reference solution

The main features of the 3D solution and solution 0 are presented in Table 1 . The deformed structure obtained with solution 0 is presented in Figure 7 . The torque due to the eccentricity of the load induces a longitudinal rotation of the cross-section, and the transverse part of the load induces a bending of the beam. In order to compare both solutions, the deflection at point A placed on the free section of the beam ( 3 = 6 ) is represented in Figure 8 during the 10 time steps of the study.

FIGURE 8 Deflection at point A

In view of the force-displacement curves shown in Figure 8 , the results obtained with the beam solution 0 are consistent with the reference solution. The beam solution satisfactorily captures the plastic branch despite a low kinematic hardening ( = 0.02 ). The curve of 0 is slightly above the curve of ref : for = 0.8 m, solution 0 associates a force 2, 42% higher than the force obtained by ref .

In order to assess the accuracy of the beam solution when compared to the reference solution the six components of the plastic strain computed by the solution 0 and by the reference solution at 3 = 0, 5 m for the given displacement of point A = 0.8 m are presented in Figure 9 . All the variables presented for a fixed displacement of point A are obtained by linear interpolations between the increments defined in Section 4.2 and represented by dots in Figure 8 .

We also define the following 2 -estimator to compare the beam model and the reference model on a cross-section S:

2 p = ⎡ ⎢ ⎢ ⎢ ⎣ ∫ p,3D - p 2 d ∫ p,3D 2 d ⎤ ⎥ ⎥ ⎥ ⎦ 1∕2 (24)
where p,3D is the plastic strain computed by the reference solution.

The axial plastic strain presented in Figures 9 a to 9 c shows that nearly all parts of the section have reached the elastic limit. As expected, the eccentricity of the load on the free extremity of the beam creates a slightly uneven progression of plasticity in the section. Therefore the highest values of the plastic strain components are observed at the top left and the bottom left of the section where the absolute values of the stresses are largest. The plastic strain computed by 0 is slightly lower than the plastic strain computed by ref for each component. A late detection of plasticity due to the longitudinal refinement is suspected to originate this phenomenon. The influence of the longitudinal mesh refinement is investigated in the next section. The values of the 2 -estimator defined previously are for the axial strains 2 The non-axial components presented in Figures 9 d to 9 f seem less satisfying but their amplitude is about 10 times lower than the axial components. Thus, the beam solution presented here shows satisfying results with a good comparison with the 3D solution. The solution 0 has been defined with a set of parameters: mesh refinement, NURBS order, expansion order, etc. The following sections investigate the influence of these parameters on the results and discuss possible ways to get more accurate results.

Mesh refinement

Solution 0 shows satisfying results in terms of displacement. Its longitudinal mesh has been chosen arbitrarily and it is now interesting to investigate the influence of the refinement of the longitudinal mesh on the results. All solutions considered here are interpolated with second-order NURBS, meshes can therefore be characterized by the knot vector used for the definition of the NURBS, as explained in Section 4.3. A reference mesh 0 is characterized by the following knot vector: {0, 1, 2, 3, 4, 5, 6}. More refined meshes are defined by adding new knots in the interval [0, 1]. These additional knots are placed at the position 3 = 1∕2 . Thus, the knot vector of mesh is 0, 1 2 , 1 2 -1 , ..., 1 2 , 1, 2, 3, 4, 5, 6 . For mesh 0 , the plastic-mode cross-section is placed at integration point 3 = 0.5 m. The same case study as before is performed for = 1..6. The parameters of the six meshes investigated are gathered in Table 2 and the deflection computed for each mesh is presented in Figure 10 . The refinement of the longitudinal mesh close to the clamped end means an improvement of the results. We define the relative distance between the reference solution and solution m for a given by:

( ) ( ) = | | | | | ( ) - ref ( ) ref ( ) | | | | | (25)
The relative distance ( ) is computed for each mesh for = 0.4 m and = 0.8 m. The results are presented in Figure 11 with the CPU time needed by the computation of each solution . ( = 0.8) cannot be computed for ≤ 1 since the computed displacement does not reach this value (see in Figure 10 ). The two curves associated to ( = 0.4) and ( = 0.8) confirms the results of Figure 10 : the more the mesh is refined, the more the results are accurate. Solutions m4 and m5 reach

Interpolation functions

The interpolation functions have a role in the quality of the results. The functions used in the present paper are NURBS as described in Corre et al. [START_REF] Grégoire | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] . The degree of interpolation of the NURBS has an incidence on the model: the higher is the degree of interpolation, the more numerous are the functions of interpolation and the broader is their support. The parameters used for the 4 solutions studied in this section are the same as the ones for solution 0 and the interpolation degree of NURBS takes values from 1 to 4. The main parameters of the 4 solutions studied here are gathered in Table 3 and the deflection at point A computed for each solution is presented in Figure 12 . The relative distances ( = 0.4) and ( = 0.8) and the CPU time required by solution are presented in Figure 13 . Solution p1 is notably less satisfying than the other solutions: the relative distance ( = 0.4) is large (12%) and the force-displacement curve does not reach = 0.8 m. These poor results were expected since locking occurs in the element for this interpolation order which corresponds to Lagrange linear elements. The increase of the order of the interpolation functions comes with an improvement in the results: while (0.4) = 2.56% and (0.8) = 2.42% for 0 , the relative distance is improved in p3 with (0.4) = 1.09% and (0.4) = 1.13%. Best values are reached with p5 : (0.4) = 0.93% and (0.4) = 0.93%. The computational cost of the solutions increases with the interpolation order. This is mainly due to the increase of the number of degrees of freedom and of the integration sections induced by a higher interpolation order. The CPU time curve presented in Figure 13 shows levels of computation time: p1 and 0 need approximately the same CPU time, and so do p3 and p4 . This is explained by the increase of the number of integration sections: this number is determined by equation ( 23) and is the same for the interpolation order 2 + 1 and 2 + 2. The small increase of CPU time between = 1 and = 2 and between = 3 and = 4 is only due to the increase of degrees of freedom. In view of the results, solutions p3 and p3 can be an alternative to 0 : the increase of the computational cost is notable but the improvement of the results is significant.

Asymptotic expansion order for the computation of the plastic modes

The number of plastic modes used in solution 0 has been set to pAE = 9. In this section, we investigate on the impact of pAE on the results. The higher is pAE , the richer should be the kinematics and the better should be the solution. Based on solution 0 , 4 other solutions are studied with pAE taking values from 0 to 12. For each solution, the new modes correspond to going one order higher in the asymptotic expansion detailed in Section 2.2.1. The main parameters of these solutions are gathered in Table 4 . The relative distance for = 0.8 m computed by each solution and the relative total time of each solution are presented in Figure 14 . The force-displacement curves at point A are not shown here since they are too close to afford a good comparison of the solutions. As expected, the highest distance from the reference solution is obtained with the solution without plastic modes: ( = 0.4) = 3.13% and ( = 0.8) = 3.05%. Its computational cost is the smallest since no time is spent on extending the kinematics with plastic modes. The increase of pAE brings a reduction of the distance from the reference solution, and comes with an increase of the computational cost. n4 0 is not better than n3 0 (respectively 2.56% versus 2.57% for = 0.4 m and 2.42% for both solutions for = 0.8 m). This means that the plastic modes computed at the fourth expansion order do not bring additional information when compared to the basis computed at the third order. However, it is interesting to note that the solution without plastic modes quickly converges and yield satisfying results. This is explained by the fact that the displacements induced by the plastic strain are already well described by the modes of the Saint-Venant's solution and the force modes in this particular case study. Based in Figure 14 , the value pAE = 9 seems relevant in the present case.

Loading -unloading example

In order to illustrate the ability of the model to deal with non-monotonic loadings, we consider the same problem loaded from 0 MN to 2.5 MN in 10 time steps and add 10 additional steps to unload the structure from 2.5 MN to 0 MN. We consider the solution 0 previously presented and the same 3D reference solution. The deflection of point A is presented in Figure 15 .

As expected, the unloading occuring between time steps 11 and 20 is elastic. The gap between the beam solution and the 3D solution induces a difference between the two residual displacements observed at 20 . For the solution 0 considered here, the relative distance between the two solutions at 20 is of 8.51%. This relative distance originates from the gap observed at 10 and is logically constant between 11 and 20 since the unloading occurs according to the elastic stiffness. A residual displacement closer to the 3D reference could be obtained with more refined longitudinal meshes and with higher interpolation order for the interpolation functions as exposed before. But in regards with the very low computational time offered by the beam solution, the results obtained can be considered satisfying for engineering applications.

We have here exposed the ability of the model to consider non-monotonic loading. The model can also handle multiple load cases and different time evolutions. The kinematics of the solution must therefore be enriched with force modes according to each load and plastic modes potentially computed from different plastic-mode cross-sections and the basis of modes needs to be updated accordingly.

CONCLUSION

A new higher-order elasto-plastic beam model has been presented. The model is based on an enrichment of the reduced kinematics on the fly during the incremental-iterative process. The kinematics of the model, represented by 2D-displacement modes, is updated at each increment according to the plastic state of the beam. This higher-order elastic beam model does not need any a priori knowledge on the solution of the problem to extend its kinematics.

The method has been applied to a cantilever beam loaded at its free extremity by a force applied with eccentricity. The beam model required a computation about 100 times shorter than the reference 3D solution computed on Code_Aster. Since the kinematics of the element is adapted to the forces applied and to the plastic state of the structure, this beam model is able to capture all the deformations induced by the load considered: bending, torsion or shear forces. The numerical differences between the beam solution and the reference solution for the example presented are low, and could still be lowered. Indeed, the computation of the displacement modes is optimal only far from boundary conditions. At the clamped extremity the local stress varies rapidly. The model could therefore be improved if displacement modes specific to the boundary conditions were computed and added to the kinematics of the model. Moreover, the gains in computation time are already very interesting and could be even larger by thanks parallelization, in particular for the computation of the stiffness matrix.

The elasto-plastic beam model has been presented with isotropic material with a 2 yield criterion. Its adaption to different yield criteria is easy. The extension of the model to anisotropic material or to more complex material like reinforced concrete could be the next step of the development of this elasto-plastic beam model.

A.2 Standard 3D plasticity algorithm

Considering the problem introduced in Section A, we recall here the procedure classically used to deal with an elasto-plastic 3D model.

The equations governing the beam evolution are gathered in equations ( 3) and (A.1). We assume that ( , 0) = ( , 0) = p ( , 0) = ( , 0) = ( , 0) = ( , 0) = 0. The time range [0, ] of the study is uniformly discretized into time steps: 0 , 1 , ..., with 0 = 0 and = . The objective of the plasticity algorithm is to compute the state of the structure = , , p , , for each time step = . The state +1 is computed from the state of the previous step and from the input +1 and +1 . The approach is iterative: the discretized equilibrium equations are solved at a global level, generating an incremental strain Δ . This strain yields new state variables { , p , }, by integrating the local constitutive equations. The global balance equation is then tested with the new stress. The iteration process is continued until the global balance is satisfied.

A.2.1 The radial return algorithm

The local plasticity algorithm exposed in this section is the well-known radial return method. The radial return algorithm was first introduced by Wilkins 33 , and was applied to 2 flow theory. The linear isotropic and kinematic hardening have been considered and introduced into the radial return algorithm by Krieg and Key [START_REF] Krieg | Implementation of a Time Dependent Plasticity Theory into Structural Computer Programs[END_REF] .

We consider the time step , +1 , and assume the state to be known. The algorithm solves the following problem: determine the state variables ep +1 , +1 , p +1 and +1 at time +1 with the knowledge of , p , and the strain increment Δ at time .

The equations ( 3) and (A.1) are discretized and expressed at = +1 :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ +1 = + tr Δ + 2 Δ -Δ p , ‖ ‖ +1 ‖ ‖ - + Δ ≤ 0, Δ p = Δ √ 3 2 +1 , Δ ≥ 0, Δ ‖ ‖ +1 ‖ ‖ -( + Δ ) = 0 (A10)
where Δ is the deviatoric part of Δ . The radial return algorithm is a well documented procedure. It is summarized here in Algorithm 2. Initialize the tangent stiffness matrix = + 1

1

 1 . The force is applied on the length Algorithm Beam global algorithm 1: Initialize state variables: 0 = 0 , 0 , p 0 , 0 , 0 2: Compute the bases S-V and ← see Section 2.2.1 3: Assemble and orthonormalize the initial basis of modes = 0 = S-V ∪

FIGURE 3

 3 FIGURE 3Dimensions (mm) of the HE600M section, mesh and applied load

3 1 2 FIGURE 4

 24 FIGURE 4 3D representation of the HE600M cantilever beam loaded at its end

FIGURE 9 FIGURE 10 FIGURE 11

 91011 FIGURE 9 Plastic strain computed by 0 and ref close to the clamped extremity at 3 = 0.5 m for = 0.8 m

FIGURE 12

 12 FIGURE 12Deflection at point A for solutions p1 to p5

2 FIGURE 15

 215 FIGURE 14and computational cost for solutions n0 0 to n4 0

Algorithm 2 1 :Algorithm 3

 213 Radial Compute trial elastic stresses trial +1 = + 2 Δ and trial +1 = + tr(Δ ) + 2 Δ end if The consistent elasto-plastic tangent moduli The consistent elasto-plastic tangent moduli is the discrete counterpart of the continuum elasto-plastic tangent moduli defined in equation (A9). The notion of consistent tangent moduli was presented in Simo and Taylor 35 , and originates in Hughes and Approximation of the global procedure Assuming a 3D discretization of the beam, the Newton-Raphson procedure is approximated with finite elements. In all the following, the finite element matrices are denoted with the notation [•], and the finite element vectors are denoted with the notation {•}. The approximation of equation (A21) yields ̂ ] dΩ = { ̂ } ep,( ) ( ) , (A22) where ep is the global tangent stiffness matrix computed from the local elasto-plastic tangent moduli ep +1 . The first member of equation (A19) is approximated as follows: -(Δ ( ) , ̂ ) = { ̂ } ( ) (A23) ( ) is the residual force, defined as the sum of the external and the internal forces: resolution of equation (A25) yields the displacement correction ( ) . The global algorithm is presented in Algorithms 3. Standard global algorithm 1: Initialize state variables: 0 = 0 , 0 ,

ep 6 :

 6 Assemble the residual{ } = extint , ref = ‖{ }‖,

TABLE 1

 1 Main features of the 3D solution and solution 0

		Solution ref	Solution 0 (2D+1D)
	type of elements	15-nodes prisms 6-nodes triangle + 19 longitudinal nodes
	number of elements	7740	399 + 1
	number of DOF	74517	3562 + 242
	CPU computation time	818 s	11 s

FIGURE 7

Deformed shape of the beam after the 10 load increments (solution 0 )

TABLE 2

 2 Main parameters of solutions m0 to m5

	Solution	m0	m1	m2	0	m4	m5
	plastic modes number pAE	9	9	9	9	9	9
	force modes number fAE	4	4	4	4	4	4
	NURBS order	2	2	2	2	2	2
	knots	7	8	9	10	11	12
	integration sections	13	15	17	19	21	23

TABLE 3

 3 Main parameters of solutions p1 to p5

	Solution	p1	0	p3	p4	p5
	plastic modes number pAE	9	9	9	9	9
	force modes number fAE	4	4	4	4	4
	NURBS order	1	2	3	4	5
	knots	10	10	10	10	10
	integration sections	19	19	37	37	55

TABLE 4

 4 Main parameters of solutions n0 to n4

	Solution	n0	n1	n2	0	n4
	plastic modes number pAE	0	3	6	9	12
	force modes number fAE	4	4	4	4	4
	NURBS order	2	2	2	2	2
	knots	10	10	10	10	10
	integration sections	19	19	19	19	19

In the following, Greek indices , , = 1,

denote cross-sectional dimensions and Latin indices , , , = 1, 2,

3, all three dimensions. Einstein summation convention on repeated indices is used.

APPENDIX A STANDARD 3D PLASTICITY

A.1 2 flow theory with isotropic hardening

We consider the beam presented in Section 2.1.

The state of the beam at a time is defined by the variables , p , and , where is an internal variable. In stress-space, the space of plastically compatible states is defined by:

where is the space of statically compatible stresses, and ∶ × ℝ + → ℝ is the yield criterion function. The state is elastic for any ( , ) ∈ × ℝ + such that ( , ) < 0. The boundary of is defined by

and is called the yield surface.

We consider an associated flow rule, therefore the plastic strain flow follows the normality flow rule:

where is a non-negative function called the consistency parameter, and ̇• = • ∕ . We can then express the Kuhn-Tucker conditions, describing the evolution of in :

The Huber-Von Mises [START_REF] Mises | Mechanik der festen Körper im plastisch-deformablen Zustand[END_REF] or 2 criterion is defined by the yield function:

where is the deviatoric part of and ‖•‖ = √ • ∶ •. We consider a linear isotropic hardening defined by:

where 0 is the yield stress, is the plastic modulus and ∶ Ω × ]0, ] → ℝ + is the equivalent plastic strain, expressed by:

The criterion used here ensures the convexity of and a smooth boundary . The Prandlt-Reuss equations then write as:

Thanks to equation (A3), the plastic flow is orthogonal to the yield surface since is the normal to the yield surface. Using equations (A4) and (A.1), the rate of change of can be expressed in terms of the total strain rate ̇ : ̇ = ep ∶ ̇ . Here ep is the elasto-plastic tangent moduli, given by:

where is the bulk modulus, is the shear modulus, is the fourth-order symmetric unit tensor, and is the second-order unit tensor.

Taylor [START_REF] Thomas | Unconditionally stable algorithms for quasi-static elasto/visco-plastic finite element analysis[END_REF] . Its expression is given by: ep +1 = -+1 .

(A11)

Here +1 = 2 Δ p Δ is a plastic correction to the elastic tensor . Its expression is established as follows:

A.2.2 The global algorithm

The objective of the plasticity algorithm is to yield the state +1 with the knowledge of the state . The local integration of equations ( 3) and (A.1) ensures that the stress computed with the algorithm 2 is plastically admissible, meaning that is inside the elasto-plastic domain . The global algorithm will ensure that is statically admissible, meaning that the global equilibrium is verified. If the evolution is plastic, both conditions are not reached simultaneously, motivating the iterative procedure. The procedure presented here is based on the Newton-Raphson algorithm.

Formulation of the balance equations and the Newton-Raphson algorithm

The balance equation writes as the weak form of the local equilibrium equations, here expressed as the time step = +1 :

where ( ± , ) is the space of displacements kinematically admissible for a zero displacement imposed on  ± . The local Algorithm 2 is represented by the function  of the variables Δ , , p and .

The main input of the algorithm is the strain increment Δ . This motivates the choice of the displacement increment Δ as the main unknown of the problem. We write:

For convenience, we only keep the first variable in the expression of stress at time +1 :

We must find Δ ∈ ( ± , ) such as the global balance is ensured, meaning that the residual is zero:

where

The iterative process consists in finding Δ satisfying equation (A17) by using a Newton-Raphson method: we iteratively correct Δ ( ) thanks to the linearized equation:

where ( ) = Δ ( +1) -Δ ( ) is the correction brought to Δ ( ) . The correction must satisfy ( ) ∈ ( ± , ). The convergence is reached for such that:

where is a scalar setting the convergence tolerance. The increment is then updated (Δ , ̂ ) = (Δ ( ) , ̂ ). Using equation (A18), the second member of equation (A19) can be written as:

The consistent elasto-plastic tangent moduli ep,( ) +1 is computed from the strain increment Δ ( ) thanks to the radial return algorithm exposed in Algorithm 2.