Ibrahim Bitar
email: ibrahim.bitar@ec-nantes.fr

Stéphane Grange

Panagiotis Kotronis

Nathan Benkemoun

A comparison of displacement-based Timoshenko multi-fiber beams finite element formulations and elasto-plastic applications

Keywords: Timoshenko, multi-fiber, displacement-based beam, shape functions

published or not. The documents may come

A comparison of displacement-based Timoshenko multi-fiber beams finite element formulations and elasto-plastic applications

Introduction

The modelling of engineering structural elements is often made using two approaches: mesoscopic (stress-strain) 2D or 3D finite element models able to capture local phenomena such as crack initiation and propagation, bond slip and aggregate interlocking (e.g. reinforced concrete beams); and global scale models (e.g. beam, plates or shell elements) where specific kinematic assumptions are adopted to simplify the global equilibrium equations and to reduce the required number of degrees of freedom.

The classical Timoshenko beam theory adopts the kinematic assumption that plane sections remain plane but not necessarily normal to the beam axis. It is often used because it can (approximately) take into account the influence of shear strains and stresses. On the contrary, the Euler-Bernoulli theory (sections remain plane and normal to the beam axis) neglects their influence and provides therefore satisfactory results for the case of slender beam structures only. Euler-Bernoulli or Timoshenko beam finite elements are either displacement-based (stiffness approach, see for example [START_REF] Nickel | Convergence of consistently derived Timoshenko beam finite elements[END_REF][START_REF] Tessler | On a hierarchy of conforming Timoshenko beam elements[END_REF][START_REF] Thomas | Finite element model for dynamic analysis of Timoshenko beam[END_REF] or force-based (flexibility approach, see for example [START_REF] Saritas | Inelastic axial-flexure-shear coupling in a mixed formulation beam finite element[END_REF][START_REF] Spacone | Fibre beam-column model for non-linear analysis of R/C frames: Part I. formulation[END_REF].

For non linear calculations, Euler-Bernoulli and Timoshenko beam finite element beams are divided in two categories: lumped beams where the material non linearity is concentrated at the element's ends; and distributed models where the inelastic deformations can spread along the element. Considering the sectional behaviour, one can distinguish the stress resultant (e.g. moment -rotation) models [START_REF] Bui | Enriched Timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF][START_REF] Dujc | Multi-scale computational model for failure analysis of metal frames that includes softening and local buckling[END_REF]Jukiće ta l . , 2013;[START_REF] Pham | Stress resultant model for ultimate load design of reinforced-concrete frames: Combined axial force and bending moment[END_REF] and multi-fiber beams, where the beam section is divided in "fibers" associated with a local strain-stress constitutive law.

During the last 20 years, various multi-fiber beam elements have been developed and implemented in research and commercial oriented finite element codes. One of the first manuscripts introducing the idea of dividing a beam section in "fibers" with a specific stress/strain relation is the book of [START_REF] Owen | Finite elements in plasticity[END_REF], followed by the works of [START_REF] Chan | Nonlinear geometric, material and time dependent analysis of reinforced concrete shells with edge beams[END_REF], [START_REF] Izzuddin | Nonlinear dynamic analysis of framed structures[END_REF], [START_REF] Scordelis | Computer models for nonlinear analysis of reinforced and prestressed concrete structures[END_REF], [START_REF] Spacone | Fibre beam-column model for non-linear analysis of R/C frames: Part I. formulation[END_REF]. Multi-fiber beams have been proven very efficient for various applications in civil engineering: the non linear analysis of beam type or bearing wall structures with non homogeneous sections (e.g. reinforced concrete) [START_REF] Grange | Numerical modelling of the seismic behaviour of a 7-story building: NEES benchmark[END_REF][START_REF] Ile | Shaking table tests of lightly RC walls: Numerical simulations[END_REF]Kotronis & Mazars, 2005;Kotronis et al., 2005), arbitrarily geometrical plane or hollow shape sections [START_REF] Desprez | Seismic vulnerability assessment of a RC structure before and after FRP retrofitting[END_REF][START_REF] Grange | The effects of soil -structure interaction on a reinforced concrete viaduct[END_REF] submitted to flexion, shear or torsion [START_REF] Mazars | Using multifiber beams to account for shear and torsion: Applications to concrete structural elements[END_REF], axial and flexure interactions [START_REF] Chan | Nonlinear geometric, material and time dependent analysis of reinforced concrete shells with edge beams[END_REF][START_REF] Izzuddin | Nonlinear dynamic analysis of framed structures[END_REF][START_REF] Pinho | Dynamic collapse testing of a full-scale four storey RC frame[END_REF][START_REF] Scordelis | Computer models for nonlinear analysis of reinforced and prestressed concrete structures[END_REF][START_REF] Spacone | Response of RC members including bond-slip effects[END_REF], soil structure interaction [START_REF] Grange | The effects of soil -structure interaction on a reinforced concrete viaduct[END_REF], vulnerability assessment [START_REF] Desprez | Seismic vulnerability assessment of a RC structure before and after FRP retrofitting[END_REF] and fiber-reinforced polymer retrofitting [START_REF] Desprez | Damage model for FRPconfined concrete columns under cyclic loading[END_REF].

The Euler-Bernoulli multi-fiber beams are used for cases where shear effects are not significant [START_REF] Spacone | Response of RC members including bond-slip effects[END_REF]. They fail therefore to account for axial, shear and flexure interactions. This is the reason why Timoshenko multi-fiber beams have been also developed: [START_REF] Garcia | Coupled model for the non-linear analysis of anisotropic sections subjected to general 3D loading[END_REF][START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF][START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF]Kotronis & Mazars, 2005;[START_REF] Mazars | Using multifiber beams to account for shear and torsion: Applications to concrete structural elements[END_REF][START_REF] Marini | Analysis of reinforced concrete elements including shear effects[END_REF][START_REF] Navarro Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF][START_REF] Petrangeli | Fiber element for cyclic bending and shear of RC structures. I: Theory[END_REF][START_REF] Ranzo | A fibre finite beam element with section shear modelling for seismic analysis of RC structures[END_REF][START_REF] Šćulac | Kinematics of layered reinforced-concrete planar beam finite elements with embedded transversal cracking[END_REF]. One has to keep in mind however that Timoshenko multi-fiber beams are unable to capture the real behaviour of structures with very small slenderness since the latter no longer behaves like a beam (see [START_REF] Mazars | A new modelling strategy for the behaviour of shear walls under dynamic loading[END_REF]. For these cases, 2D or 3D finite element models [START_REF] Ile | Nonlinear analysis of reinforced concrete shear wall under earthquake loading[END_REF] or other adequate simplified modelling strategies should be preferred [START_REF] Mazars | A new modelling strategy for the behaviour of shear walls under dynamic loading[END_REF].

Another important aspect concerns the material constitutive behaviour adopted in the context of multi-fiber beams. Different choices can be found in the literature: damage models (see for example Kotronis & Mazars, 2005), plasticity models [START_REF] Dujc | Multi-scale computational model for failure analysis of metal frames that includes softening and local buckling[END_REF], coupled plasticity and damage models (Ibrahimbegovićetal., 2007), simplified models for confined concrete columns [START_REF] Desprez | Damage model for FRPconfined concrete columns under cyclic loading[END_REF] or approaches based on the modified compression field theory introduced by Vecchio andCollins (1986, 1988).

The main purpose of the article being the influence of the different shape functions interpolating the displacement field on the element response, we focus hereafter on displacement based Timoshenko multi-fiber beams with simple uniaxial linear elastic or elasto-plastic constitutive laws. The shear strains are constant in the beam section, the shear stresses are supposed linear elastic and are multiplied by the classical shear correction factor k (Bathe, 2006;[START_REF] Cowper | The shear coefficient in Timoshenko's beam theory[END_REF]. Force-based beam elements -see for example (Ciampi &[START_REF] Ciampi | A nonlinear beam element for seismic analysis of structures[END_REF][START_REF] Spacone | Fibre beam-column model for non-linear analysis of R/C frames: Part I. formulation[END_REF] -are out of the scope of the article. Furthermore, non linear axial, shear and flexure interactions are not specifically addressed, the interested readers can among others refer to Buietal.(2014), [START_REF] Ceresa | Flexure-shear fiber beam-column elements for modeling frame structures under seismic loading state of the art[END_REF], Ceresa, Petrini, Pinho, &Sousa (2009, Saritas and[START_REF] Saritas | Inelastic axial-flexure-shear coupling in a mixed formulation beam finite element[END_REF], [START_REF] Xu | Axial-shear-flexure interaction hysteretic model for RC columns under combined actions[END_REF], [START_REF] Zona | Finite element models for nonlinear analysis of steelconcrete composite beams with partial interaction in combined bending and shear[END_REF]).

In the following, we present and compare the performance of three (3) displacementbased straight Timoshenko multi-fiber beam finite elements under the small rotations assumption. The formulations differ on the shape functions interpolating the displacements and rotations fields. The first formulation, called hereafter (Full-Linear-Independent) FLI, adopts linear polynomials while the transverse displacements and rotation fields are interpolated independently, see [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF][START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF][START_REF] Pegon | A Timoshenko simple beam element in Castem[END_REF]. Because of the low degree of the adopted polynomials, this formulation is prone to shear locking problems (see [START_REF] Beirão Da Veiga | Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods[END_REF][START_REF] Ceresa | A fibre flexure-shear model for seismic analysis of RC-framed structures[END_REF][START_REF] Crisfield | Nonlinear finite element analysis of solids and structures[END_REF][START_REF] Ville De Goyet | L'analyse statique non linéaire par la méthode des éléments finis des structures spatiales formées de poutres à section non symétrique[END_REF]Ibrahimbegović&Frey, 1993;[START_REF] Mukherjee | Convergence properties and derivative extraction of the superconvergent Timoshenko beam finite element[END_REF][START_REF] Rakowski | The interpretation of the shear locking in beam elements[END_REF][START_REF] Reddy | On locking-free shear deformable beam finite elements[END_REF][START_REF] Stolarski | Membrane locking and reduced integration for curved elements[END_REF][START_REF] Stolarski | Membrane locking and reduced integration for curved elements[END_REF][START_REF] Yunhua | Explanation and elimination of shear locking and membrane locking with field consistence approach[END_REF] if specific measures are not considered. A way to avoid this is to calculate the shear strains approximately, as proposed by [START_REF] Donea | A modified representation of transverse shear in C 0 quadrilateral plate elements[END_REF] and adopted in [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF] and [START_REF] Pegon | A Timoshenko simple beam element in Castem[END_REF]. In the second formulation, referred hereafter as (Full-Cubic-Quadratic-Material) FCQM, proposed by [START_REF] Friedman | An improved two-node Timoshenko beam finite element[END_REF] for Timoshenko beams (further bibliography on Timoshenko beam finite elements can be found in [START_REF] Kiendl | Single-variable formulations and isogeometric discretizations for shear deformable beams[END_REF], [START_REF] Litewka | The exact thick arch finite element[END_REF], [START_REF] Nickel | Convergence of consistently derived Timoshenko beam finite elements[END_REF], Papa Dukiće ta l .( 2014), [START_REF] Tessler | On a hierarchy of conforming Timoshenko beam elements[END_REF] and [START_REF] Thomas | Finite element model for dynamic analysis of Timoshenko beam[END_REF] and adopted in Kotronis and Mazars (2005) and [START_REF] Mazars | Using multifiber beams to account for shear and torsion: Applications to concrete structural elements[END_REF] in a multi-fiber context, the degree of the shape functions is of order three (3) for the transverse displacements and two (2) for the rotations. Displacements and rotations are now interdependent and the shear locking problem is avoided. The particularity of this formulation is that the shape functions depend on material properties, and therefore their use for non linear problems can be problematic. The third formulation, called hereafter (Full-Cubic-Quadratic) FCQ, is proposed by [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF] and uses shape functions of order three (3) for the transverse displacements and two (2) for the rotations and an additional internal node. This results to a finite element free of shear locking and to shape functions independent on the material properties. The choice to compare these formulations is motivated by the fact that the first two elements FLI and FCQM are implemented, respectively, in the general purpose finite element codes (Aster, Retrieved from http:// code-aster.org/ and Castem, Retrieved from http://www-cast3m.cea.fr/). The more recent FCQ element is chosen in order to study its performance.

The article is organised as follows: in the first section, the main equations of the Timoshenko beam theory are recalled as well as the general form of the displacement and strain fields, the stiffness matrices and the consistent nodal force vectors within a multifiber beam approach. In the second section, the specific equations for the FLI, FCQM and FCQ formulations are provided. The performance of each formulation is studied in the third section considering an elastic or elasto-plastic behaviour and monotonic loadings. Of particular interest is the last example where we study the specific case of non linear axial-bending interactions. Shape functions of order three (3) are proposed for the axial displacements of the FCQ element in order to increase its accuracy (Section 4.5.1). The article ends with general conclusions and guidelines.

Timoshenko multi-fiber beams

Kinematic assumptions

Consider the case of a 2D straight beam, with a length Łand a normal cross-section S(x), its middle axis = (0, Ł) oriented towards the x direction. The small rotations hypothesis is adopted and the classical Timoshenko assumptions are made (the cross-section S(x) remains plane but not necessary orthogonal to the deformed middle axis). The generalised displacements vector U (x) of the cross-section S(x) -G(x, 0) being a point of the section S(x) on the middle axis, see Figure 1 -takes the following form (in the following, the upper index T states for "transpose", the term "generalised" and the capital symbols U and refer to the section):

U (x) = U x (x) U y (x) z (x) T (1)
U x (x) being the longitudinal displacement, U y (x) the transverse displacement and z (x) the rotation of the section S(x) where the point G(x, 0) is situated (see Figure 1). The displacements u x (x, y), u y (x, y) of another point P(x, y) of the section S(x) ( o ro fa "fiber" f (x, y)) are evaluated using the section displacements as follows:

u x (x, y) = U x (x) -y z (x), u y (x, y) = U y (x) (2)
The strain field becomes (the symbol ′ defines hereafter the derivative with respect to x):

ε x (x, y) = ∂u x ∂ x = U ′ x (x) -y ′ z (x), γ xy (x, y) = ∂u x ∂ y + ∂u y ∂ x = U ′ y (x) -z (x) = β y (x) (3)
with ε x (x, y) the axial strain and γ xy (x, y) the transverse shear strain of the fiber f (x, y).

Constitutive laws

At each fiber f (x, y) of the section S(x) a constitutive law can be adopted (e.g. based on elasticity, on damage mechanics or plasticity) linking the axial and the shear terms. For the case of an isotropic elastic material, it takes the following form:

σ x (x, y) = E f (x, y)ε x (x, y), τ xy (x, y) = kG f (x, y)γ xy (x, y) (4)
where E f (x, y) and G f (x, y) are the Young and the Shear modulus, respectively, and k is the shear correction factor adopted in Timoshenko beam finite elements (due to the fact that the resulting shear stresses do not respect the boundary conditions on the beam section Bathe, 2006;[START_REF] Cowper | The shear coefficient in Timoshenko's beam theory[END_REF]. The values of the shear correction factor for different crosssectional shapes can be found in [START_REF] Cowper | The shear coefficient in Timoshenko's beam theory[END_REF]. A discussion on the shear correction factors can be also found in [START_REF] Dong | Much ado about shear correction factors in Timoshenko beam theory[END_REF].

In the following and for sake of clarity and simplicity the dependencies on x or x, y are omitted.

Virtual work principle

The virtual work principle becomes [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF] and [START_REF] Pegon | A Timoshenko simple beam element in Castem[END_REF]

: Ł 0 S (δε x σ x + δγ xy τ xy )dSdx -w external = 0( 5 )
where the first term is related to the work of the internal loadings and w external the work of the external loadings.

Introducing (3)in( 5), yields

Ł 0 S (δ(U ′ x -y ′ z )σ x + (δU ′ y -δ z )τ xy dSdx -w external = 0( 6 )
where δU x is the virtual axial displacement, δU y the virtual transversal displacement and δ z the virtual rotation.

The following definitions are adopted for the generalised forces:

F x = S σ x dS, F y = S τ xy dS, M z =-S yσ x dS Axial force Shear force Bending moment (7)
introducing them in (6) one obtains:

Ł 0 F x d dx δU x + F y d dx δβ y + M z d dx δ z dx -w external = 0( 8 )
where

β y = dU y dx -z = U ′ y -z (9)
Considering ( 3),( 4) and ( 7) gives:

F x = S σ x dS = S E f ε x dS = S E f dU x dx -y d z dx dS = S E f dSU ′ x - S E f ydS ′ z F y = S τ xy dS = S kG f γ xy dS = S kG f (β y )dS = S kG f dSβ y M z =- S yσ x dS =- S yE f ε x dS =- S yE f dU x dx -y d z dx dS =- S E f ydSU ′ x + S E f y 2 dS ′ z ( 10
)
and the generalised force vector F S of the section can be written as:

⎡ ⎣ F x F y M z ⎤ ⎦ F S = ⎡ ⎣ S E f dS 0 -S E f ydS 0 S kG f dS 0 -S E f ydS 0 S E f y 2 dS ⎤ ⎦ K S . ⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ = K S . ⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ , (11)
where U ′ x ,β y , ′ z are the generalised strains and K S is the multi-fiber section stiffness matrix, independent of the neutral axis position [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF].

Finally, considering ( 8) and ( 11), one gets for the virtual work principle:

Ł 0 δ U ′ x β y ′ z K S U ′ x β y ′ z T dx -w external = 0 (12)
The extension of the previous equations in 3D is straightforward and follows the same steps [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF].

Finite element discretisation

The beam of length Ł is discretised into n finite element beams e =[x i ; x j ] of length L e = x j -x i and external nodes i and j. The generalised displacement vector is approximated by an equation of the form U = NU e as follows:

⎡ ⎣ U x U y z ⎤ ⎦ = ⎡ ⎣ N 1 N 2 N 3 N 4 N 5 N 6 N 7 N 8 N 9 N 10 N 11 N 12 N 13 N 14 N 15 N 16 N 17 N 18 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 13
)
where U e is a vector containing the external nodal displacements of the element e and N is the matrix of the shape functions depending on x.

The generalised strains are evaluated as:

⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ = ⎡ ⎣ U ′ x U ′ y -z ′ z ⎤ ⎦ = ⎡ ⎣ B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8 B 9 B 10 B 11 B 12 B 13 B 14 B 15 B 16 B 17 B 18 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎣ B axial B shear B flexion ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 14
)
where B is a matrix containing the derivatives with respect to x of the shape functions.

Introducing ( 14) in the virtual work expression (12) provides the expression of the stiffness matrix K element of the element e [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF]:

K element = L e 0 B T K S Bdx (15)
while the assembly of the latter gives the total stiffness matrix of the beam structure:

K tot = n A e=1 L e 0 B T K S Bdx (16)
The calculation of the internal nodal forces due to a state of (internal) stresses follows the same procedure:

F int,element = L e 0 B T F S dx =⇒ F int,tot = n A e=1 L e 0 B T F S dx (17)
The consistent nodal force vector of the finite element beam e due to a distributed vertical load p(x) (Figure 2), is again evaluated using the virtual work principle. The external virtual work in (12) for one (1) element e takes the following form: with N T shear = N 7 N 8 N 9 N 10 N 11 N 12 and the consistent element nodal force vector that becomes:

w external = L e 0 δU y (x) p(x)dx = L e 0 δ(N shear U e ) T p(x)dx = δU T e L e 0 N T shear p(x)dx = δU T e × F ext,element (18)
F ext,element = L e 0 N T shear p(x)dx = L e 0 N 7 N 8 N 9 N 10 N 11 N 12 T p(x)dx (19)

FLI, FCQM and FCQ formulations

We present hereafter the shape functions, the stiffness matrices and the consistent nodal force vectors for the FLI, FCQM and FCQ formulations.

FLI formulation

Linear independent shape functions are adopted for the generalised displacement fields [START_REF] Guedes | A fibre/Timoshenko beam element in Castem[END_REF][START_REF] Pegon | A Timoshenko simple beam element in Castem[END_REF]:

⎡ ⎣ U x U y z ⎤ ⎦ = ⎡ ⎣ N 1 00N 4 00 0 N 8 00 N 11 0 00N 15 00N 18 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (20)
where

⎧ ⎪ ⎨ ⎪ ⎩ N 1 = N 8 = N 15 = 1 - x L e N 4 = N 11 = N 18 = x L e (21)
The adopted generalised strain field has the following form:

⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ = ⎡ ⎢ ⎣ -1 L e 00 1 L e 00 0 -1 L e -1 2 0 1 L e -1 2 00 -1 L e 00 1 L e ⎤ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 22
)
where following the proposal of [START_REF] Donea | A modified representation of transverse shear in C 0 quadrilateral plate elements[END_REF], the expression of β y is modified by eliminating the linear terms in the shape functions N 15 and N 18 (β y component) in order to avoid the shear locking problem. The section and the element stiffness matrices are evaluated according to ( 11) and ( 15). For the case of an homogeneous elastic section (a classical Timoshenko elastic beam) the element stiffness becomes:

K element = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ES L e 00 -ES
-EI L e + kSGL e 4 0 -kSG 2 EI L e + kSGL e 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 23
)
where E and G the Young modulus and the Shear modulus, I and S the quadratic moment and the section surface, respectively. The consistent element nodal force for a n-order polynomial p(x) = qx n , using ( 19) and ( 21), is found as:

F ext,element = qL n+1 e (n + 1)(n + 2) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 n + 1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 24
)

FCQM formulation

Cubic and quadratic interdependent polynomials, functions of the material properties, are used for the transverse and rotational displacements and the formulation is free of shear locking [START_REF] Friedman | An improved two-node Timoshenko beam finite element[END_REF]. The generalised displacement and strain fields are:

⎡ ⎣ U x U y z ⎤ ⎦ = ⎡ ⎣ N 1 00 N 4 00 0 N 8 N 9 0 N 11 N 12 0 N 14 N 15 0 N 17 N 18 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 25
) ⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ = ⎡ ⎣ N ′ 1 00 N ′ 4 00 0 N ′ 8 -N 14 N ′ 9 -N 15 0 N ′ 11 -N 17 N ′ 12 -N 18 0 N ′ 14 N ′ 15 0 N ′ 17 N ′ 18 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 26
)
and the shape functions:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N 1 = 1 - x L e N 4 = x L e N 8 = 1 1 + φ 2 x L e 3 -3 x L e 2 -φ x L e + 1 + φ N 9 = L e 1 + φ x L e 3 -2 + φ 2 x L e 2 + 1 + φ 2 x L e N 11 =- 1 1 + φ 2 x L e 3 -3 x L e 2 -φ x L e N 12 = L e 1 + φ x L e 3 -1 - φ 2 x L e 2 - φ 2 x L e N 14 = 6 (1 + φ)L x L e 2 - x L e N 15 = 1 1 + φ 3 x L e 2 -(4 + φ) x L e + 1 + φ N 17 =- 6 (1 + φ)L x L e 2 - x L e N 18 = 1 1 + φ 3 x L e 2 -(2 -φ) x L e , (27)
with φ the ratio of the beam bending to shear stiffness provided by (ν the Poisson's ratio):

φ = 12 L 2 e EI kGA = 24 L 2 e I kA (1 + ν). (28)
For the case on an elastic homogeneous material, the element stiffness matrix and the consistent element nodal force vector for a n-order polynomial p(x) = qx n become:

K element = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ES L e 00 -ES L e 00 0 12EI (1+φ)L 3 e 6EI (1+φ)L 2 e 0 -12EI (1+φ)L 3 e 6EI (1+φ)L 2 e 0 6EI (1+φ)L 2 e (4+φ)EI (1+φ)L e 0 -6EI (1+φ)L 2 e (2-φ)EI (1+φ)L e -ES L e 00 ES L e 00 0 -12EI (1+φ)L 3 e -6EI (1+φ)L 2 e 0 12EI (1+φ)L 3 e -6EI (1+φ)L 2 e 0 6EI (1+φ)L 2 e (2-φ)EI (1+φ)L e 0 -6EI (1+φ)L 2 e (4+φ)EI (1+φ)L e ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (29) F ext,element = qL n+1 e (1 + φ)(D) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 6n + 12 + φ(n 2 + 7n + 12) L e (n + 1)[2 + φ 2 (n + 4)] 0 (n + 1)[(n + 2)(n + 6) + φ(n + 3)(n + 4)] -L e (n + 1)[n + 2 + φ 2 (n + 4)] ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 30
)

FCQ formulation

Cubic functions are used to interpolate the transverse displacements and quadratic for the rotations. The element is free of shear locking and uses an additional internal node [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF]. [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF] proved also that the FCQM formulation can be derived from the present formulation and that one FCQ element is able to predict the exact tip displacements for any complex loading (shear/flexion) subjected to an homogeneous elastic beam (see also the numerical examples of Sections 4.2 and 4.3). The latter property is related to the well-known Tong's proof [START_REF] Tong | A variational principle and the convergence of a finite-element method based on assumed stress distribution[END_REF]. Indeed, the interdependent nature and the high order of the selected polynomials [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF] satisfy the differential equations of equilibrium associated with the homogeneous form of the Timoshenko's beam theory [START_REF] Timoshenko | On the correction for shear of the differential equation for transverse vibrations of prismatic bars[END_REF]. The reader is referred to [START_REF] Zienkiewicz | The finite element method: The basics[END_REF]f o r more details.

The nodal displacement field takes the following form:

U e = U xi U yi zi U 1 yk k U 2 yk U xj U yj zj T , (31)
where U 1 yk , k and U 2 yk are the degrees of freedom of the internal node 'k' (with no specific physical meaning).

The generalised displacement field is:

⎡ ⎣ U x U y z ⎤ ⎦ = ⎡ ⎣ N 1 00000 N 7 00 0 N 11 0 N 13 0 N 15 0 N 17 0 00N 21 0 N 23 000N 27 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U 1 yk k U 2 yk U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (32)
where

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N 1 = 1 - x L e N 7 = x L e N 11 = 1 - x L e 2 1 + 2 x L e N 13 = 2 1 - x L e 2 x L e N 15 =-2 x L e 2 1 - x L e N 17 = x L e 2 3 -2 x L e N 21 = 1 - x L e 1 -3 x L e N 23 = 1 -1 -2 x L e 2 N 27 =- x L e 2 -3 x L e (33)
This gives for the generalised strain field:

⎡ ⎣ U ′ x β y ′ z ⎤ ⎦ = ⎡ ⎣ N ′ 1 00000 N ′ 7 00 0 N ′ 11 -N 21 N ′ 13 -N 23 N ′ 15 0 N ′ 17 -N 27 00 N ′ 21 0 N ′ 23 000 N ′ 27 ⎤ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U xi U yi zi U 1 yk k U 2 yk U xj U yj zj ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 34
)
For the case on an elastic homogeneous material, the element stiffness matrix and the consistent element nodal force vector for a n-order polynomial p(x) = qx n become:

K element = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ES Le
+ 4 EI Le ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 35
)
F ext,element = qL n+1 e D ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 6(n + 2) 0 4(n + 1) 0 -2(n + 1)(n + 2) 0 (n + 1)(n + 2)(n + 6) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (36)
where

D = (n + 1)(n + 2)(n + 3)(n + 4).
The introduction of internal degrees of freedom leads to 9 × 9 matrices (rather than 6 × 6) and to 9 × 1 vectors. To use this element to a general purpose finite element code, one has to calculate the internal degrees of freedom using two approaches: either the internal degrees of freedom are treated as external degrees of freedom and are sent to the global solver of the finite element code either they are calculated locally (inside the element subroutine) using the static condensation method (see [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF] for more details and the analytical expressions of the condensed matrices and vectors).

Comparison of the FLI, FCQM and FCQ formulations

Preliminaries

In order to study the accuracy of the different formulations, the analytical solution for the case of an homogeneous Timoshenko cantilever beam structure of length Ł subjected to vertical punctual force F y at the free end (x = Ł) and to a polynomial p(x) = qx n distributed vertical force takes the following form (for more details see [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF] Figure 3: For the numerical examples presented hereafter, numerical integrations are performed in each section and along the axis of the multi-fiber beams. Each section is discretised using 64 triangular finite elements (T 6), where each T 6 has 6 nodes and 3 Gauss integration points. In this way, the section stiffness matrices and the section force vectors are exactly integrated for the case of an elastic material behaviour [START_REF] Bitar | Un nouvel élément fini poutre pour la modélisation des structures sous sollicitations complexes[END_REF]) (e.g. the S E f y 2 dS component in the section stiffness matrix (11) requires at least 3 integration points as the interpolation functions for T 6 are of order 2). The numerical integration of Equation ( 15) along the beam x axis is performed according to the Gauss integration scheme using one (1) Gauss's point for the FLI, two (2) for the FCQM (the components of the matrix B Equation ( 26) are linear) and three (3) for the FCQ formulation (the derivatives of the transverse displacement interpolation functions being quadratic (33), the exact numerical integration of 4) dx (15) requires at least three Gauss integration points).

z = F y 2EI x 2 -2xŁ + 1 EI qŁ n+2 n + 2 x - qŁ n+1 2(n + 1) x 2 + q (n + 1)(n + 2)(n + 3) (37) U y = F y 2EI x 3 3 -x 2 Ł - F y x kGS + q kGS(n + 1) Ł n+1 x - x n+2 n + 2 + q EI Ł n+2 2(n + 2) x 2 - Ł n+1 6(n + 1) x 3 + x n+4 (n + 1)(n + 2)(n + 3)(n + 4) (38)
L e 0 B T K S B O(
Different numerical examples chosen to illustrate the performance of the FLI, FCQM and FCQ beam formulations are Figure 4: 4.2 A cantilever elastic Timoshenko beam subjected to a transverse tip displacement. 4.3 A cantilever elastic Timoshenko beam subjected to a polynomial p(x) = qx n vertical force. 4.4 A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement. 4.5 A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force.

Remark: In the FCQM formulation, the shape functions depend on the material properties ( 27),( 28). The performance of the multi-fiber Timoshenko beam for non linear calculations can therefore be undpredictable [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF]. This is the reason why in the following, the FCQM formulation is not tested with the elasto-plastic applications 4.4 and 4.5.

A cantilever elastic Timoshenko beam subjected to a transverse tip displacement

The geometrical characteristics of the beam structure and the material parameters are given in Figure 5 and Table 1.

A displacement v Ł is applied at x = Ł. The resulting analytical shear force at x = 0is given as (see also Equation ( 38) considering q = 0 and x = Ł) (Association française de normalisation): and the corresponding moment:

F ana y = v Ł Ł 3 3EI + Ł kGS (39)
M ana z = Ł × F ana y (40)
Comparison of the performances of the three formulations is provided in Tables 2 and3 where n is the number of the multi-fiber Timoshenko beam elements used for the spatial discretisation and R E (F y ) the relative error on shear forces defined as:

R E (F y ) = F ana y -F num y F ana y (41)
with F num y the result of the numerical calculations. One (1) FCQM or FCQ element provides the exact solution in terms of forces and moments, whereas the FLI formulation presents a relative error R E (F y ) of 32%, see Table 2. This is due to the fact that its shape functions corresponding to the transverse displacements and rotations are linear, whereas the analytical solutions of z and U y are a second-and a third-order polynomial, respectively (see Equations ( 37), (38) considering q = 0). By increasing the number of elements, the accuracy of the FLI formulation however quickly improves (see Table 3). Table 4. A cantilever elastic Timoshenko beam subjected to a polynomial p(x) =-20x 2 vertical force. Results of the three formulations considering one (1) multi-fiber Timoshenko beam element.

F y (KN) M z (KNm)v ( 10 -3 )θ ( 10 -3 ) R E (θ) (%)
n v(10 -3 ) R E (v)(%) θ(10 -3 ) R E (θ)(%) 1 -0.

A cantilever elastic Timoshenko beam subjected to a polynomial p(x)

= qx n vertical force A 2-order distributed load p(x) =-20x 2 is applied on the previous cantilever beam. The vertical displacements and rotations at the free end for the different finite element formulations and for different discretisations are given in Tables 4 and5. The analytical values are calculated from (37) and (38) considering F y = 0, q =-20 and n = 2, see also [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

Again, one (1) FCQM or FCQ multi-fiber Timoshenko beam provides the exact results. Because of the adopted linear shape functions, one (1) FLI element gives a relative error on rotations R E (θ ) (defined in a similar way as (41)) of 25%. Since the calculation is forced controlled, numerical results on shear forces and moments are exact for all formulations. The accuracy of the FLI formulation improves rapidly with the increase of the number of elements; the relative error drops from 25% (one element) to 1.7% (four elements), Table 5.

A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement

We consider hereafter an elastic perfectly plastic material [START_REF] Simo | Computational inelasticity[END_REF] assuming that only the normal component of the axial stress σ x can enter to plasticity. No interaction between the shear stresses and the normal stresses at the material constitutive law level is considered. The stress elastic limit f y is taken equal to 450 MPa, while the other parameters are provided in Table 1. The plastic moment is (Eurocode 2 De Normalisation, 2004): and the plastic shear force:

M pl = W pl f y = bd 2 4 f y (42)

FLI FCQ

nF y (KN) M z (KNm) R E (F y ) (%) F y (KN) M z (KNm) R E (F y ) (%)
F yp = M pl Ł (43)
These two values are obtained considering that all the fibers in the section are plastified. This is an asymptotic section behaviour state that cannot be reached with the adopted kinematic assumption of plane sections, since the axial strain at the neutral axis equals zero and thus the corresponding fibers cannot plastified. (One should notice that the shear correction factor k = 5/6 derived from elasticity is adopted in this example, since the shear behaviour is considered elastic).

The cantilever beam structure is subjected to an imposed displacement v Ł at its free end (x = Ł). Comparison of the performances of the two formulations for different discretisations are shown in Tables 6 and7. As before, one (1) FLI element does not provide good results. However, increasing the number of elements greatly improves its performance. Results are better for the FCQ formulation that provides the smallest error among the two formulations.

Figure 6 presents the moment diagrams along the beam for the two formulations (FLI and FCQ) and for different number of elements. The moment nodal values are plotted and linked with straight lines. By increasing the number of multi-fiber beam elements, the moment diagrams of the two formulations coincide. Conclusions are similar in Figure 7 showing the evolution of the shear force F y with respect to the vertical tip displacement v Ł . Notice that the graphs (6), ( 7), but also (8), ( 9) and (10) illustrate the extrapolated nodal values computed from the values at the integration points.

In order to illustrate the mesh dependency problems that can occur when a constitutive law without an internal length parameter is adopted [START_REF] Armero | An analysis of strain localization and wave propagation in plastic models of beams at failure[END_REF], Figures 8 and9 present the axial strain ε x (3) in an extreme fiber of a section and the section rotation z (13) along the length of the beam for the two formulations and two discretisations. In all cases a plastic hinge, dependent on the size of the elements, is found concentrated in the first multi-fiber beam element. Different approaches have been proposed in the literature to correctly simulate strain localisation problems: introducing rate dependency [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF], non-local theories [START_REF] Bažant | Softening in Reinforced Concrete Beams and Frames[END_REF], higher order media [START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: Localization studies[END_REF] Figure 6. A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement -moment diagrams along the beams for 2 formulations considering n multi-fiber Timoshenko beam elements.

Figure 7. A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement -shear forces vs. vertical displacements for 2 formulations considering n multi-fiber Timoshenko beam elements. and the so-called enhanced approaches (Armero &[START_REF] Armero | An analysis of strain localization and wave propagation in plastic models of beams at failure[END_REF][START_REF] Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF] among others. In [START_REF] Armero | An analysis of strain localization and wave propagation in plastic models of beams at failure[END_REF] and in order to regularise the solution, the authors proposed to incorporate the localised dissipative mechanisms in the form of strong discontinuities in the generalised displacements. Localisation issues were also treated within the context of multi-fiber beams, the readers are referred to [START_REF] Coleman | Localization issues in force-based frame elements[END_REF] for force-based, and to Jukićetal.(2014) and [START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF] for displacement-based multi-fiber elements. Nevertheless, strain localiSation issues are out of the scope of the present article.

One should notice however that mesh dependency problems do not appear if an elastoplastic constitutive law with a positive isotropic hardening K H is adopted. The previous example is revisited using K H = 20 MPa and results are shown in Figure 10. One can see that the solution is not spuriously localised at the left end.

A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force

The purpose of this last example is to study the ability of the proposed multi-fiber beam element formulations to reproduce the interaction between normal forces and bending moments (see [START_REF] Casaux | Modélisation tridimensionnelle du comportement sismique d'ouvrages en béton armé: développement de méthodes simplifiées [3D modeling of the seismic behavior of reinforced concrete structures: development of simplified methods[END_REF][START_REF] Casaux-Ginestet | On multi-fiber beam model with coupled axial and bending response[END_REF][START_REF] Casaux-Ginestet | On multi-fiber beam model with coupled axial and bending response[END_REF]. We consider hereafter the simple beam of Figure 11 having an homogeneous cross-section and an elastic perfectly plastic material. The geometric characteristics of the beam and the material properties are given in Table 8. At first, an horizontal increasing force is applied at its right end till a maximal value F max . Then, and while the axial force remains constant, an increasing rotation θ is applied at its left end.

As illustrated in Figure 11, the beam structure is discretised with two elements, such that L element1 = 0.8 m and L element2 = 1.2 m. 3 Gauss points are used for the numerical Figure 8. A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement -axial strain in an extreme fiber vs. beam length. Figure 9. A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement -section rotation vs. beam length. Remark: As already mentioned, only one Gauss point is necessary for the FLI numerical integration. Three Gauss points are however used in this particular example in order to compare with the results at the same sections of the FCQ element.

According to the equilibrium equations, since d F x dx = 0 the axial force should remain constant. The numerical expression of the normal force is F

x (x) = S E f dSU ′ x (x) - S E f ydS ′ z (x) (
10). Linear shape functions are used to interpolate the axial displacements Figure 13. A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force -normal forces at the three sections within the first element for the finite element formulations FLI and FCQ.

U x (x) (21) and ( 33) leading to a constant value for U

′ x (x) = U xj -U xi L e
along the element and this for the two formulations. The FLI formulation adopts linear shape functions to interpolate z (x) and therefore

′ z (x) = zj -zi L e remains also constant. A constant value is thus found for F x (x) = S E f dS U xj -U xi L e -S E f ydS zj -zi L e
along the FLI element as shown in Figure 13.

Results are however different for the FCQ formulation. The shape functions used to interpolate z (x) are quadratic. The expressions of ′ z (x) can be derived from the corresponding Equation (34) as follows:

′ z (x) = N ′ 21 zi + N ′ 23 i + N ′ 27 zj = 2 L e -2 + 3x L e zi + 4 L e 1 - 2x L e i + 2 L e -1 + 2x L e zj (44)
Equation ( 44) shows that ′ z (x) is linear along the element. This is the reason why the evolution of F x (x) in Figure 13 is different, depending on the element section. For example, in Section 2 which is situated at x = L e 2 we have

′ z L e 2 = zj -zi L e
, and therefore the same result with the FLI formulation is obtained. As shown in Figure 14, increasing the number of the FCQ elements improves the results.

Table 9 presents a brief comparison of the three beam formulations: " * " states for poor results, " * * " for results close to the analytical solution, " * * * " for results that coincide with the analytical solution.

A more satisfactory solution is proposed in the next section based on the kinematic enhancement of the axial displacement field.

Kinematic enhancement of the axial displacement field

Following a similar idea as in the work of [START_REF] Casaux | Modélisation tridimensionnelle du comportement sismique d'ouvrages en béton armé: développement de méthodes simplifiées [3D modeling of the seismic behavior of reinforced concrete structures: development of simplified methods[END_REF] for multi-fiber Euler-Bernoulli beams, a specific novel kinematic enhancement is proposed for the axial displacements in order to improve the ability of the FCQ Timoshenko beam to deal with non linear axial-bending interactions. Cubic shape functions can be used for the interpolation of the x L e ;

U x = N xi N 1 xk N 2 xk N xj ⎡ ⎢ ⎢ ⎣ U xi U 1 xk U 2 xk U xj ⎤ ⎥ ⎥ ⎦ , (45)
where U 1 xk and U 2 xk are the additional internal degrees of freedom that as before can solved at the global or at the element level (static condensation). Figure 16 depicts the response of the enhanced FCQ formulation for the previous example using only two elements. It is obvious that the non linear axial-bending interactions are now better reproduced.

Conclusion

In this paper, a comparison between three different Timoshenko multi-fiber beam formulations was presented: the FLI formulation with linear shape functions, the FCQM formulation with higher order shape functions dependent on the material properties and the FCQ formulation with higher order shape functions and additional internal degrees of freedom. Numerical examples showed that the FLI formulation does not provide accurate results when only one element is adopted. Results are however improved using a finer discretisation. The FCQM formulation is not suitable for non linear calculations but its performance is for elastic linear cases. Finally, the FCQ formulation provides the best results with the smaller number of finite elements. In the end, a novel kinematic enhancement of the axial displacement field of the FCQ element is proposed, in order to increase its ability to simulate non linear axial-bending problems. Further applications are necessary to study the performance of the multi-fiber Timosehnko beams for non linear axialshear-bending interactions. The authors are currently working on improving the multi-fiber Timoshenko FCQ beam introducing the embedded discontinuity approach at the sectional and fiber level to better reproduce failure of reinforced concrete structure.

Figure 1 .

 1 Figure 1. Kinematics of a Timoshenko beam.

Figure 2 .

 2 Figure 2. A finite element beam e subjected to a vertical distributed load p(x).

Figure 3 .

 3 Figure 3. Cantilever beam under distributed load and vertical punctual force.

Figure 4 .

 4 Figure 4. Gauss's points distribution in a section.

Figure 5 .

 5 Figure 5. Cantilever beam.

Figure 10 .

 10 Figure 10. A cantilever elastic plastic (with positive hardening) Timoshenko beam subjected to a transverse tip displacement -section rotation vs. beam length.

Figure 12 .

 12 Figure 12. A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force.

Figure 14 .

 14 Figure14. A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force -Normal force at the 3 sections within the first element for the FCQ formulation using different discretisations.

Figure 15 .

 15 Figure 15. New axial shape functions illustration.

Figure 16 .

 16 Figure16. A simple elastic perfectly plastic Timoshenko beam subjected to a rotation and an axial force -Normal force at the 3 sections within the first element for the enhanced FCQ formulation.

Table 1 .

 1 Cantilever beam -geometrical characteristics and material properties.

	Geometry			Material properties
	Ł1 .53 m	Young modulus	E	210 GPa
	b	0.25 m	Poisson's ratio	ν	0.3
	d	0.25 m	Shear correction factor	k	5 6

Table 2 .

 2 A cantilever elastic Timoshenko beam subjected to a transverse tip displacement. Results of the three formulations considering one (1) multi-fiber Timoshenko beam element.

		n	v Ł (m)	F y (KN)	M z (KNm)	R E (F y )(%)
	Analytical	1	0.1	5609.1	8581.9	-
	FLI	1	0.1	7428.3	11365.33 2
	FCQM	1	0.1	5609.1	8581.90
	FCQ	1	0.1	5609.1	8581.90

Table 3 .

 3 A cantilever elastic Timoshenko beam subjected to a transverse tip displacement. Results of the FLI formulation considering n multi-fiber Timoshenko beam elements.

	nF y (KN)	M z (KNm)	R E (F y )(%)	nF y (KN)	M z (KNm)	R E (F y )(%)
	4	5696.3	8715.3	1.5	20	5612.2	8587.2	0.053
	8	5630.6	8614.9	0.37	31	5610.5	8584.1	0.025
	16	5614.5	8590.1	0.087	51	5609.6	8582.73	0.009

Table 5 .

 5 A cantilever elastic Timoshenko beam subjected to a polynomial p(x) =-20x 2 vertical force. Results of the FLI formulation considering n multi-fiber Timoshenko beam elements.

	Analytical	23.9	27.4	-0.2776	-0.2453	-
	FLI	23.9	27.4	-0.2411	-0.3066	25
	FCQM	23.9	27.4	-0.2776	-0.2453	0
	FCQ	23.9	27.4	-0.2776	-0.2453	0

Table 6 .

 6 A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement. Results of two formulations considering one multi-fiber Timoshenko beam.

		n	v Ł (m)	F y (KN)	M z (KNm)	R E (F y )(%)
	Analytical	1	0.1	1148.9	1757.8	-
	FLI	1	0.1	2263	3462	97
	FCQ	1	0.1	1581	2417	37.6

Table 7 .

 7 A cantilever elastic perfectly plastic Timoshenko beam subjected to a transverse tip displacement. Results of 2 formulations considering n multi-fiber Timoshenko beam elements.

Table 8 .

 8 Simple beam -geometrical and material characteristics.

	Length	Ł	2 m	Section	S	0.2m × 0.2m
	Young modulus	E	200000 MPa	Hardening modulus	K	1%
	Stress elastic limit	f y	500 MPa	Shear correction factor	k

integration along the length of each finite element beam (in order to compare the section response at three different positions within the element). Each section is discretised with 16 fibers. The maximum axial load is taken equal to F max = 19 MN Figure12(b) while the history of the imposed rotation θ is shown in Figure12(a). Figure13presents the evolution of the normal force at the three sections within the first element. The three sections are located at the positions of the three Gauss points used for the numerical integration. The performance of the finite element formulations FLI and FCQ is studied hereafter.

interaction (F x & M) * * * * * * *

Acknowledgements

The authors would like to thank the PIA-SINAPS@ project (Sésisme et Installation Nucléaire: Améliorer et Pérenniser la Sureté), approved and funded by the National Agency of Research (ANR) following the RSNR 2012 call for projects on future investments post-Fukoshima. (SINAPS@ ANR-11-RSNR-0022).

Funding

This work was support by National Agency of Research (ANR) following the RSNR 2012 call for projects on future investments post-Fukoshima [SINAPS@ ANR-11-RSNR-0022].

Disclosure statement

No potential conflict of interest was reported by the authors.


Fabrice Barbe

Rodolphe Parisot

Samuel Forest
email: samuel.forest@mat.ensmp.fr

Georges Cailletaud

Calibrating a homogenized polycrystal model from large scale FE computations of polycrystalline aggregates

Homogenization-based estimates of the effective behaviour of elastoviscoplastic polycrystals usually resort to approximations like that of an inclusion embedded in the homogeneous equivalent medium in the case of the self-consistent scheme. We propose here to replace this idealized intermediate problem by the actual computation of a polycrystalline representative volume element. Computations of a cube made of 200 zinc grains are presented and used to calibrate the scale transition rule of a simplified homogenized polycrystal model. A coupled identification procedure is presented to determine both intragranular material parameters and scale transition parameters from experimental tension, creep and relaxation tests.

INTRODUCTION TO COMPUTATIONAL HOMOGENIZATION METHODS

Homogenization methods are now a well-established tool of the mechanics of heterogeneous materials to predict the effective properties from the knowledge of the constituents, at least in the linear regime. In elastoplastic and purely viscoplastic cases, bounds and estimates are available, at least for simple loading paths [START_REF] Suquet | Continuum Micromechanics[END_REF]. Estimates like the well-known self-consistent scheme provide the values of the mean stress and strain in each individual phase. In a metal polycrystal, a phase is an ensemble of grains having the same crystallographic orientation up to a given solid angle. According to the self-consistent scheme, the state of each phase is evaluated using the solution of the problem of an inclusion of this phase embedded in the wanted homogeneous equivalent medium (HEM). This latter problem admits no exact analytical solution in the non-linear case and one usually takes a homogeneous tangent or secant operator for the matrix as an approximation [2]. In [3], the authors resort to finite element (FE) computations to solve the inclusion problem and tune some parameters of the homogenized model. The problem of the inclusion embedded in the unknown HEM can also be seen as an implicit problem for the effective properties and an inverse approach can be used to determine the properties of the HEM. In [4] a simple J 2 -plasicity constitutive framework is adopted for the HEM and the finite element method is used iteratively to calibrate the parameters and fulfill the self-consistency condition. A similar inverse problem is solved in [START_REF] Forest | IUTAM Symposium on micromechanics of plasticity and damage of multiphase materials[END_REF] and [START_REF] Pilvin | Multiaxial Fatigue and Design[END_REF] but a more complex form of the constitutive equations is retained for the HEM including explicit scale transition rules [START_REF] Pilvin | Mai[END_REF]. Again, inclusion problems are solved using the FE method to calibrate the parameters involved in the scale transition rule. A scale transition rule is a relation giving the mean stress (or strain) in a phase as a function of the mean macroscopic stress (or strain) and other internal variables. The use of simplified explicit transition rules is illustrated in section 4.

Currently available computing capabilities make it now possible to abandon the approximation of the grain as an inclusion in an infinite matrix and to compute instead a volume element of polycrystal containing enough grains to regard it as representative of the polycrystal. Such computations of polycrystalline representative volume elements can be performed in the three-dimensional case relying on parallel computing [START_REF] Quilici | [END_REF]9]. This approach provides the full heterogeneous stress and strain fields inside the grains including the grain boundary effects. The mean stress and strain can be computed over the whole volume or over each phase and can be compared to the prediction of simplified homogenization models. In particular such computations can be regarded as reference calculations to assess the quality of new estimations of bounds from homogenization theory.

In the present work we propose to calibrate simplified homogenization models directly from the prediction of computations of polycrystalline aggregates. Simplified models are still necessary since the computations of polycrystalline aggregates are lengthy and require parallel computing. They should however include no explicit yield surface at the macroscopic level and instead incorporate simplified scale transition rules in order to be able to be able to account for large distorsion of the effective yield surface. We apply this methodology to the case of an elastoviscoplastic material for which accurate estimates remain seldom.

MECHANICAL PROPERTIES OF BULK ZINC

The elastoviscoplastic behaviour of a zinc alloy used as a coating material in galvanized steel sheets has been investigated at room temperature [10]. It has been possible to produce bulk zinc ingots that have been rolled and from which plate samples have been machined. A heat treatment has been chosen to produce a grain size of about 40µm. The material is anisotropic due to the typical zinc rolling texture and only tests in the rolling direction are considered here. Tensile tests at different strain rates, creep and relaxation tests have been performed; some results are given in Fig. 1 together with simulations that will be explained in section 5. Because of its low melting temperature, zinc displays a strongly viscoplastic behaviour at room temperature. Because of the alloying elements, the material is found to be much harder than pure zinc.

Zinc has a closed-packed hexagonal crystallographic structure. Two main slip systems have been found to be activated during the deformation of bulk zinc, namely basal slip on the basal plane ¡ 0001¢ in direction a £ 1120 ¤ (3 slip systems) and pyramidal Π 2 slip on plane ¥ 1122¦ in direction c § a £ 1123 ¤ (6 slip systems). These systems are included in the constitutive modelling. The constitutive equations describing the intragranular behaviour of zinc are the following [START_REF] Quilici | [END_REF] :

ε © ε ¨e § ε ¨p σ © C ¨¨: ε ¨e ε ¨p © N ∑ s 1 γs ¡ m s n s § n s § m s ¢ 2 (1) γs © Max 0 ! τ s " x s " τ s c K s # n s $ sign ¡ τ s " x s ¢ τ s © σ ¨: ¡ m s n s ¢ (2)
The analysis is restricted to small strains. The total strain is split into elastic and plastic parts and the stress is computed from Hooke's law [START_REF] Suquet | Continuum Micromechanics[END_REF]. A Schmid law is adopted. The slip rate γs is given by a Norton type equation with a threshold τ s c . τ s denotes the resolved shear stress on slip system s. The N slip systems include both basal and pyramidal ones. Two hardening variables are introduced : non linear isotropic hardening τ s c and kinematic hardening x s according to :

τ s c © τ s 0 c § N ∑ r 1 Q s h rs % 1 " exp& b r γ r cum ' γs cum © γs (3) x s © c s α s αs © γs " d s α s γs " x s M s # m s sign ¡ α s ¢ (4)
The internal variable x s is necessary to better account for the relaxation tests and its evolution law (4) contains dynamic and static recovery terms. The isotropic hardening evolution law (3) describes the non linear evolution of the critical resolved shear stress as a function of cumulated slip γ cum and involves an interaction matrix : basal π 2

( Q s h rs) © basal π 2 0 1 1 1 1 1 1 1 1 2 Q b 0 1 2 h b 1 h b 2 3 43 53 h b 2 h b 1 3 43 53 h b 2 h b 2 . . . 6 87 9 Q b h b π 2 Q π 2 h b π 2 Q π 2 0 1 2 h π 2 1 h π 2 2
3 53 53

h π 2 2 h π 2 1 3 53 53 h π 2 2 h π 2 2
. . . 6 87 9 6 87 7 7 7 7 7 7 7

9

(5) Each slip @ system familly has its own set of material parameters K n Q b c d M m, but in order to reduce the number of unknowns, a constraint on the ratio of the respective parameters is imposed in the following identification procedure :

Q bas Q π 2 © τ bas 0c τ π 2 0c © K bas K π 2 © τ bas 0c τ π 2 0c © c bas c π 2 © 1 5 (6)
n b d M m are the same for all slip systems. Furthermore, we take [10] :

h b 1 © h b 2 2 © h π 2 1 © h π 2 2 2A 2 © h π 2 2 2A 6A a.
B DC FE HG PI RQ S TC FU WV YX abX X ac ed f hg B iC FE pG PI RQ eS TC FU V qX arX WX X ac ed f hg s ut av ws yx DC E s yV S X abX X ac ed f hg s ut v ws yx DC FE s yV S X arX WX X ac ed f hg r

X abX X rX X rX X arX W X abX X rX ac X c d X c e X c eX WX f X WX X X X b.
¦ § © ª « 8¬ u¬ 8¬ ¬ u¬ 8¬ ® 8¬ u¬ 8¬ 8¬ 8¬ 8¬ °¬ 8¬ 8¬ ± u¬ 8¬ 8¬ ² s¬ 8¬ 8¬ ³ i¬ 8¬ 8¬ ¬ ³ °¬ ³ D² s¬ ³ i¬ 8¬ « 8¬ ® 8¬ °¬ ² s¬ ¬ Figure 1.
Viscoplastic behaviour of bulk zinc : tensile tests at two strain rates (a); creep test at 80 MPa followed by a jump to 100MPa (b); relaxation tests at three different strain levels (c).

COMPUTATION OF POLYCRYSTALLINE AGGREGATES

The three-dimensional distribution of the grains in a sample can be observed experimentally for instance in the case of multicrystalline specimens containing a small number of large grains [11]. In the general case it still remains almost out of experimental reach. This is why the polycrystalline volume element is generated numerically as a distribution of Voronoi polyhedra [9]. A multiphase element technique is used to superpose a finite element mesh on the image of the microstructure : the same crystallographic orientation is attributed to all integration point inside a given grain. In the present work, an isotropic random set of orientations has been randomly distributed among 200 grains. Such a grain distribution and the corresponding mesh are given in Fig. 2a. In the present work, we do not try to represent the anisotropic overall behaviour of the considered material. Furthermore, homogeneous isotropic elasticity is considered, which is of course actually not the case for zinc grains. The observed macroscopic value of Young's modulus in the rolling direction (130000MPa) is retained. The constitutive equations (1) to [START_REF] Forest | IUTAM Symposium on micromechanics of plasticity and damage of multiphase materials[END_REF] are inte ´grated at each Gauss point. It is essential to have enough elements in each grain in order to have a converged solution in the sense of mesh refinement. This requires at least 1000 integration points per grain. This is why we resort to parallel computing for the computation of polycrystalline aggregates [11]. The RVE computed in this work is a 18x18x18 quadratic brick cube containing 200 grains and divided into 6 subdomains, each subdomain being computed on a single processor. a. It is clear that the number of grains remains too small for this volume to be really representative. However, in the case of f.c.c. elastoplastic crystals, it has been found sufficient to reach a stabilized response in tension @ when increasing the number of grains [8]. This has not been checked in the case of h.c.p. crystal yet. Moreover some effects due to the application of boundary conditions can be expected. Again in the case of f.c.c. crystals, it has been shown that the boundary layer effect extends over about one grain, provided a fine enough mesh is used, and even less in the case of free surfaces. In the case of 200 grains, boundary layer effects on the overall curves can therefore not be excluded.

Two tensile tests and a shear test have been computed on this aggregate, for a given set of material parameters. For tension, the vertical displacement is prescribed on two opposite faces of the cube, the remaining faces being free of forces. This corresponds to a mean prescribed axial deformation from 0 to 1%. The distribution of basal slip and pyramidal slip activity is given in Fig. 2b and 2c for mean strain E © 0A 5%. A very strong heterogeneity is obtained with pyramidal slip occuring mainly at the grain boundaries where plastic incompatibilities develop. This is corroborated by the metallographic observations of Fig. 3 which shows slip lines near some grain boundaries that have been identified as traces of pyramidal slip planes [10]. This is due to the fact that the initial critical resolved shear stress for pyramidal slip is much larger than for basal slip. Note that the observations have been made at a free surface of a specimen and that the view of Fig. 2 also corresponds to free lateral surfaces. The computation gives a similar effect in the core of the aggregate. The mean axial stress Σ vs. mean axial strain E curves are given in Fig. 5, for two different tensile strain rates, together with other simulation results to be explained in section 5.

EFFECTIVE POLYCRYSTAL MODEL

Homogenized polycrystal models are available to estimate the mean stress-strain in the phases of the polycrystal in the case of elastoplastic, purely viscoplastic and even elastoviscoplastic grain behaviour, for instance according to the self-consistent scheme [START_REF] Suquet | Continuum Micromechanics[END_REF]2]. In the polycrystal model proposed in [START_REF] Pilvin | Multiaxial Fatigue and Design[END_REF][START_REF] Pilvin | Mai[END_REF] an explicit scale transition rule is introduced linking the mean local stress σ ¨g in grain g to the macroscopic stress Σ ¨, the mean viscoplastic deformation E ¨p and the mean local viscoplastic deformation ε ¨pg according to :

E ¨© E ¨e § E ¨p Σ ¨© C ¨¨: E ¨e Σ ¨© N g ∑ g 1 f g σ ¨g Ė ¨p © N g ∑ g 1 f g ε ¨pg (7) σ ¨g © Σ ¨ § C % B ¨" β ¨g' B ¨© N g ∑ g 1 f g β ¨g β ¨g © ε ¨pg " Dε pg eq β ¨g (8)
where homogeneous elasticity is assumed. The volume fraction of phase g is denoted by f g , N g is the number of discretized orientations (phases). Non linear accommodation variables β ¨g are introduced, that gradually ¶ depart from ε ¨pg according to the non linear evolution rule in [START_REF] Quilici | [END_REF]. For specific values of parameters C and D, classical formulation like Sachs, Lin or Kröner polycrystal models are retrieved. In [START_REF] Forest | IUTAM Symposium on micromechanics of plasticity and damage of multiphase materials[END_REF][START_REF] Pilvin | Mai[END_REF], the scale transition parameters C and D have been tuned to mimic the response of the self-consistent scheme. C varies then from 0 to the order of magnitude of the shear modulus. Once the stress σ ¨g is computed in each grain, the single crystal constitutive equations ( 1) to ( 5) are used to compute the viscoplastic strain rate ε ¨pg . In the next section, we will calibrate them to account for the actual response of the aggregate. The question is whether the form of the equations ( 7) and ( 8), which is not far from standard homogenization models for elastoplastic polycrystals, will still be suited to describe the elastoviscoplastic behaviour of a zinc aggregate.

OPTIMIZATION PROCEDURE

Two group of parameters arose in the course of the previous sections : material parameters describing the intragranular behaviour and scale transition parameters appearing in the scale transition rule [START_REF] Quilici | [END_REF]. The second group can be identified using the response of a polycrystalline aggregate with given material parameters. The homogenized model must then be able to reproduce the macroscopic behaviour given in section 2. The identification procedure here requires simulation with both models (aggregate and homogenized polycrystal) and the computation of two cost functions : the first one • 1 measures the difference between the mean response of the aggregate and the corresponding simulation with the homogenized model, whereas the second one • 2 penalizes the difference between this simulation and the experimental curve. An optimization procedure enables one then to find the best set of material and scale transition parameters fulfilling both requirements. In the process of comparing the RVE computation and the homogenized estimate, one should not only check that the macroscopic response is corrected but also that the mean stress/strain in each phase are correctly reproduced by the homogenized model.

A computation with 6 processors of the considered aggregate lasts 5 days whereas a simulation with the homogenized polycrystal model needs a few minutes. This is why a slightly different strategy is adopted represented in Fig. 4. In fact we decouple the two identification procedures. A first trial for C D is used to optimize the material parameters to minimize • 2 . Once an acceptable fit is obtained, one checks the validity of C D by computing a polycrystalline aggregate with the new material parameters, which gives C 1 . A new estimation of C D is then proposed and one goes back to the first iterative procedure. In the present work, we have carried out these three stages only once. A tensile test at 0A 001s & 1 has been used to identify C D, an acceptable fit has been obtained (Fig. 5). As a validation, a tensile test at 0A 0001s & 1 has been simulated with the aggregate and the mean stress-strain curve has been compared to the simulation of the homogenized model. Figure 5 shows that the viscous effect is significantly underestimated by the homogenization model.

The bunch of curves of Fig. 6a and 6b gives the mean value of axial stress vs. axial strain for all 200 grains for both models. It appears that stress heterogenity is much more pronounced in the aggregate than predicted ¸by the homogenization model. In fact, the number of grains in the aggregate should be large enough to have several representatives of the same given orientations, to be able to compute a mean value over several grains, which is the actual meaning of the phase in the homogenization model. This has been done for f.c.c. crystals and the result is closer to the homogenized solution [9]. In the present case of h.c.p. crystals, four specific grains are compared for both simulations in Fig. 7a. Finally a shear test has been performed on the aggregate using homogeneous deformation conditions at the boundary, at 0A 0001s & 1 .

The mean shear stress-strain curve is given in Fig. 7b and compared to the prediction of the homogenized polycrystal model, which underestimates the attained shear stress levels.

CONCLUSIONS AND PROSPECTS

At the end of this first attempt of a coupled identification, the description of the macroscopic tests remains rather poor : creep deformation is overestimated and the relaxation is too rapid and too small (Fig. 1). The difficulty of the homogenized model to reproduce the different reponses of the aggregate may have two origins :

¹ the chosen aggregate is too small and a larger number of grains should be considered; ¹ the scale transition rule ( 8) is not sufficient to describe elastoviscoplastic materials : improvements involving 1 to 3 more scale transition parameters are proposed in [START_REF] Pilvin | Mai[END_REF][START_REF] Quilici | [END_REF] that should be tested.

Finally, a better description of the experimental results could be achieved by relieving the constraint on material parameters [START_REF] Pilvin | Multiaxial Fatigue and Design[END_REF]. Note that the real situation of the zinc alloy requires to take anisotropic local elasticity and the crystallographic texture into account. This will probably require bigger aggregates and a more sophisticated scale transition rule. The methodology will remain the same and the present work is a first attempt in this direction which remains to be significantly improved. Absolute grain size effects are not accounted for in this work, so that the found values of materials and scale transition parameters are bound to one specific grain size. However the method could be extended to size-sensitive crystal plasticity models like [12]. Comparison between the equivalent stress-strain curves obtained by both models for 4 grains (left); overall shear stress-strain curve for the aggregate and prediction of the homogenized model.

  g ih rj k ml rn Do ph bq sr t vu iw vw vx yo pw vg io {z s| i} 5| s| u n o ph bj w g 5| s| s| 8| |

Figure 2 .

 2 Figure 2. Distribution of 200 grains (a); equivalent basal slip distribution for a mean axial deformation of 0µ 004 (b, the red color corresponds to an amount of slip greater than 0µ 01); equivalent pyramidal Π 2 slip (c, the red color corresponds to an amount of slip greater than 0µ 003).

Figure 3 .

 3 Figure 3. Experimental evidence of basal and pyramidal Π 2 slip in zinc grains (grain size about 40µm).

Figure 4 .

 4 Figure 4. Optimization schemes for the identification of material and scale transition parameters.

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Simulation of tension tests on the polycrystalline aggregate and with the homogenized polycrystalline model.