Ibrahim Bitar
email: ibrahim.bitar@live.com

Nathan Benkemoun

Panagiotis Kotronis

Stéphane Grange

A multifiber Timoshenko beam with embedded discontinuities

An enhanced high order multifiber Timoshenko beam is introduced to simulate structural behavior up to failure. The beam is displacement based and geometrically linear, its section can be of arbitrary shape and a local constitutive law is assigned to each fiber. The Strong Discontinuity Approach is adopted to enhance the displacement field of the fibers to describe crack openings. The material behaviour at the discontinuity is characterized by a linear cohesive law linking the axial stress and the displacement jump. The variational formulation is presented in the context of the incompatible modes method and details are given on the corresponding computational procedure and the numerical integration of the constitutive laws. The simulation of the non linear behavior of a cantilever beam structure and of a reinforced concrete frame are provided to illustrate the performance of the novel enhanced high order multifiber Timoshenko beam.

Introduction

Experimental observations show that under severe static or dynamic loadings strains localize in specific zones sometimes called 'plastic hinges'. With increasing loading, discontinuities can appear leading to partial or total structural collapse. Among the different existing numerical models to study post-peak behavior and failure, beams permit to reduce the necessary number of degrees of freedom and are therefore computational very efficient.

A novel enhanced high order multifiber Timoshenko beam is introduced in this article. The beam is displacement based and geometrically linear, its section can be of arbitrary shape and a local constitutive law is assigned to each fiber. The Strong Discontinuity Approach (SDA) is adopted and the material behaviour at the discontinuity is characterized by a linear cohesive law linking the axial stress and the displacement jump. When a significant amount of fibers reaches failure, failure at the beam and structural level can be therefore 'naturally' reproduced.

Several finite element multifiber beam formulations have been developed and implemented in various Finite Element codes [START_REF] Chan | Nonlinear geometric material and time dependent analysis of reinforced concrete shells with edge beams[END_REF][START_REF] Scordelis | Computer models for nonlinear analysis of reinforced and prestressed concrete structures[END_REF][START_REF] Izzuddin | Nonlinear dynamic analysis of framed structures[END_REF][START_REF] Spacone | Fibre beam-column model for non-linear analysis of R/C frames: part I. formulation[END_REF]. Multifiber beams have proven highly effective for civil engineering applications: non-linear analysis of beams or walls with non-homogeneous sections (e.g. reinforced concrete (RC)) [START_REF] Kotronis | Simplified modelling strategies to simulate the dynamic behaviour of {R}/{C} walls[END_REF][START_REF] Kotronis | A simplified modelling strategy for R/C walls satisfying PS92 and EC8 design[END_REF][START_REF] Ile | Shaking table tests of lightly {RC} walls: numerical simulations[END_REF][START_REF] Grange | Numerical modelling of the seismic behaviour of a 7-story building: NEES benchmark[END_REF]; Soil-Structure Interaction problems [START_REF] Grange | The effects of Soil-Structure Interaction on a reinforced concrete viaduct[END_REF]; Fiber Reinforcing Polymer confinement [START_REF] Desprez | Damage model for FRP-confined concrete columns under cyclic loading[END_REF]; seismic vulnerability assessment of retrofitted RC structures [START_REF] Desprez | Seismic vulnerability assessment of a RC structure before and after FRP retrofitting[END_REF]; sections submitted to bending, shearing or torsion [START_REF] Mazars | Using multifiber beams to account for shear and torsion[END_REF], flexure-shear interaction [START_REF] Feng | An efficient fiber beam-column element considering flexure-shear interaction and anchorage bond-slip effect for cyclic analysis of RC structures[END_REF], axial and bending interactions [START_REF] Chan | Nonlinear geometric material and time dependent analysis of reinforced concrete shells with edge beams[END_REF][START_REF] Scordelis | Computer models for nonlinear analysis of reinforced and prestressed concrete structures[END_REF][START_REF] Izzuddin | Nonlinear dynamic analysis of framed structures[END_REF][START_REF] Spacone | Response of RC members including bond-slip effects[END_REF][START_REF] Pinho | Dynamic collapse testing of a full-scale four storey RC frame[END_REF]. Euler-Bernoulli multifiber beam formulations are used when the shear effects are negligible [START_REF] Spacone | Response of RC members including bond-slip effects[END_REF]. Timosheko formulations are more suitable to reproduce the interactions between axial forces, shear forces and moments [START_REF] Kotronis | Simplified modelling strategies to simulate the dynamic behaviour of {R}/{C} walls[END_REF][START_REF] Mazars | Using multifiber beams to account for shear and torsion[END_REF][START_REF] Ranzo | A fibre finite beam element with section shear modelling for seismic analysis of RC structures[END_REF][START_REF] Petrangeli | Fiber element for cyclic bending and shear of RC structures. I: theory[END_REF][START_REF] Marini | Analysis of reinforced concrete elements including shear effects[END_REF][START_REF] Garcia | Coupled model for the non-linear analysis of anisotropic sections subjected to general 3D loading. Part 1: theoretical formulation[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] Šćulac | Kinematics of layered reinforced-concrete planar beam finite elements with embedded transversal cracking[END_REF][START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF][START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

1

In order to realistically reproduce strain localization, some authors propose to replace the part of the beam where strain concentration is expected by 3D finite elements. This requires the definition of proper boundary conditions to satisfy compatibility with the adjacent beam elements [START_REF] Ian | Coupling 1d beams to 3d bodies[END_REF]. An energetic equivalence is therefore considered between the work done at the nodes of the volume interface and the single node of the beam element. Another method, less expensive in terms of calculation time, is to modify the postpeak material behavior at the fiber level [START_REF] Giry | Modélisation objective de la localisation des déformations et de la fissuration dans les structures en béton armé[END_REF].A3 Dfinite element modelling is first performed to calculate the energy dissipated during the development of the stain localization zone. The energy dissipated by the multifiber model is then considered equal by suitably modifying the constitutive law. This makes possible to derive an estimate of the crack opening in a fiber. A kinematic enhancement of the Timoshenko beam can be also adopted. In [START_REF] Mazars | Using multifiber beams to account for shear and torsion[END_REF], the authors proposed to enhance the kinematics of the multifiber section of a Timoshenko beam by a torsional warping function calculated solving of a local problem. [START_REF] Capdevielle | A multifiber beam model coupling torsional warping and damage for reinforced concrete structures[END_REF] propose to update the torsion warping function at each calculation step. A multifiber kinematic enhancement is also introduced in [START_REF] Khoder | Enhancement of a two-dimensional multifibre beam element in the case of reinforced concrete structures for taking into account the lateral confinement of concrete due to stirrups[END_REF] to take into account the effect of transverse reinforcements in reinforced concrete sections.

In order to simulate failure and to describe crack openings with classical (no multifiber) beams, several authors enhanced the kinematic using the Strong Discontinuity Approach (SDA). The material behaviour at the discontinuity is characterized by a generalized law, linking for example the bending moment and the rotational jump, which allows capturing the released fracture energy [START_REF] Armero | An analysis of strain localization and wave propagation in plastic models of beams at failure[END_REF][START_REF] Ehrlich | Finite element methods for the analysis of softening plastic hinges in beams and frames[END_REF][START_REF] Armero | Numerical modeling of softening hinges in thin Euler-Bernoulli beams[END_REF][START_REF] Dujc | Multi-scale computational model for failure analysis of metal frames that includes softening and local buckling[END_REF][START_REF] Pham | Stress resultant model for ultimate load design of reinforced-concrete frames: combined axial force and bending moment[END_REF][START_REF] Jukić | Embedded discontinuity finite element formulation for failure analysis of planar reinforced concrete beams and frames[END_REF][START_REF] Bui | Enriched Timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF][START_REF] Nikolic | Brittle and ductile failure of rocks: embedded discontinuity approach for representing mode i and mode ii failure mechanisms[END_REF][START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF]. A generalized higher order Timoshenko beam with embedded rotation discontinuity has been recently presented by the authors of this article in [START_REF] Bitar | A generalized timoshenko beam with embedded rotation discontinuity[END_REF].

Up to now however, few authors have tried to combine SDA and multifiber beams. In [START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF], the authors developed a multifiber Timoshenko beam with a discontinuity of the axial displacement at the fiber level. In [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF], the authors presented an enhanced multifiber Timoshenko beam for reinforced concrete structures. The displacement field within the fiber is enhanced by one strong discontinuity variable to describe the local failure. Concrete is modeled by means of a one-dimensional damage model coupled to a cohesive model and steel with a one-dimensional elasto-plastic model coupled to a cohesive model. In the following, the model of [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF][START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF] is referred as Full-Linear-Independent (FLI) since linear functions are used to interpolate the displacement fields at the beam level.

We present in this article the way to combine SDA (rotation discontinuity) with a higher order multifiber Timoshenko beam. The original beam formulation has shape functions of order three for the transverse displacements and two for the rotations. It is free of shear locking and one element is able to predict the exact tip displacements for any complex distributed loadings and any suitable boundary conditions [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF]. In the following, it is stated as Full-Cubic-Quadratic (FCQ) and its performance has been already compared with respect to other beam formulations found in the literature [START_REF] Pegon | A Timoshenko simple beam element in Castem[END_REF][START_REF] Guedes | A Fibre/Timoshenko beam element in Castem[END_REF][START_REF] Friedman | An improved two-node Timoshenko beam Finite Element[END_REF] in one of the previous articles of the authors [START_REF] Bitar | A comparison of displacement-based timoshenko multi-fiber beams finite element formulations and elasto-plastic applications[END_REF].

This article is structured as follows: Section 2 presents the enhanced formulation of the FCQ multifiber Timoshenko beam and describes the discontinuity kinematics. The variational formulation, the way to obtain the equilibrium equations and specific details on the determination of the enhancement functions are also provided. Section 3 introduces the constitutive laws and Section 4 the computational procedure. Two numerical applications are finally studied in Section 5 to show the performance of the novel enhanced high order multifiber beam.

2. Enhanced multifiber Timoshenko Full-Cubic-Quadratic beam 2.1. Timoshenko FCQ beam 2.1.1. Modified shape functions considering an internal rotational degree of freedom Fig. 1 illustrates a Timoshenko FCQ beam element denoted e (for simplicity reasons, presentation is provided hereafter in 2D). We note L e the length of the beam, i and j the external nodes and k the internal node.

At time t, the generalized displacements U S of a section S located at position x of the beam element axis are: U x being the longitudinal displacement, V y the transverse displacement and Θ z the rotation of the section S. d e is the nodal dis- placement vector of the FCQ formulation defined by: and V Δ yk 2 are the degrees of freedom of the internal node (with no specific physical meaning) [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

= UV V V UV d [ ΘΔ Δ Θ Δ Θ ]
Fig. 1. Timoshenko FCQ beam [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

NN , u v and N θ are the shape functions of the three displacement components defined as [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF]:

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ x x x NN NNNN NN N N N N () () () 00000
where With respect to the original formulation [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF], a slight modification is introduced hereafter re-writing the interpolation of the rotations as a function of the rotations at the nodes i and j and the rotation at the note k situated at the middle of the element = t ΘΘ ( , )

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ =- = =- ⎛ ⎝ + ⎞ ⎠ =- ⎛ ⎝ ⎞ ⎠ =- ⎛ ⎝ -⎞ ⎠ = ⎛ ⎝ - ⎞ ⎠ = ⎛ ⎝ -⎞ ⎠ ⎛ ⎝ - ⎞ ⎠ =- - =-⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎜ - ⎞ ⎠ ⎟ ( ) ( ) () () ( )
u x L u x L v x L x L v x L x L v x L x L v x L x L θ x L x L θ x L θ x L x L
zk z L 2 e
(see Fig. 2). To do this, the relation between ΔΘ zk and Θ zk is first found calculating the rotation at the middle of the element from Eq. ( 1) (expressing the shape functions N N , θθ [START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF] and

N θ 6 at = x L 2 e ): ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ == - + L tt t t t Θ 2 , Θ () Δ Θ () 1 4 [Θ ( ) Θ ( )] z e zk zk zi zj (5)
Therefore,

=+ + tt tt ΔΘ( ) Θ( ) 1 4 [Θ ( ) Θ ( )] zk zk zi zj (6)
Introducing Eq. (6) in Eq. ( 1), the rotation field can now be interpolated function of the rotations at the external notes i j , and the rotation at the middle point k:

∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ == = * * xt N x t x t Nd Θ, ( ) Θ ( ) ( )( ) z αi j k α θ zα θ e ,, (7)
with the three modified interpolation functions:

=+ = - + =+ = -+ == - N x Nx Nx N x Nx Nx Nx Nx () () () 1
() () () () () θθ θ x L x L θθ θ x L x L θθ x L x L 3 *
Eq. ( 8) satisfy all the necessary properties of a shape function [START_REF] Bitar | A generalized timoshenko beam with embedded rotation discontinuity[END_REF][START_REF] Zienkiewicz | Finite element methodthe basis[END_REF][START_REF] Bitar | Modeling the failure of reinforced concrete structures using generalized and multifiber finite element beams[END_REF]. The above modification helps to determine the enhancement functions of the multifiber Timoshenko FCQ beam element (see Section 2.3) and has no influence on the performance of the FCQ formulation [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF][START_REF] Bitar | A comparison of displacement-based timoshenko multi-fiber beams finite element formulations and elasto-plastic applications[END_REF].

The use of linear interpolation functions for the axial displacement and quadratic approximation for the rotation provides inconsistent results when N-M interactions and nonlinear material effects are considered. A solution to this problem was provided by the authors in [START_REF] Bitar | A comparison of displacement-based timoshenko multi-fiber beams finite element formulations and elasto-plastic applications[END_REF] using higher order interpolation functions for the axial displacement field for the FCQ beam. In this paper however we focus on M-T interactions and therefore the initial version of the formulation is adopted [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

Axial displacement: fiber vs. beam axis interpolation

Following the planar section hypothesis (Timoshenko theory), the displacements of a fiber of coordinates x y , at instant t are deduced from the generalised displacements (Eq. ( 1)) as: ) xx z yy [START_REF] Grange | The effects of Soil-Structure Interaction on a reinforced concrete viaduct[END_REF] The multifiber beam formulation can be presented in two equivalent: an interpolation at the fiber level (Path a) and an interpolation at the beam level axis (Path b), illustrated in Fig. 3 and Eq. [START_REF] Desprez | Damage model for FRP-confined concrete columns under cyclic loading[END_REF].

=- = uxy t Uxt y xt vxyt Vxt (, , ) (, ) Θ(, ) (, , ) (,
⇒ ⎫ ⎬ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪ ⇒ ⎫ ⎬ ⎪ ⎭ ⎪ ↘ ⇒ ⎫ ⎬ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪ ⇒ ⎫ ⎬ ⎭ ↗ ⇒⇒ ⇒⇒ ⇒ Ut Vt
K a () () Θ() () () Θ() (, ) (, ) (, ) (, ) () () Θ() () () Θ() (, ) (, ) Θ( , ) (, , ) (, , ) (, , ) (, , ) (, , ) (, , ) () (
The first formulation, Path (a) in Fig. 3a, consists in calculating the axial displacements vector yt d (, ) f at the nodes of each fiber f of ordinate y considering the hypothesis of planar sections [START_REF] Grange | The effects of Soil-Structure Interaction on a reinforced concrete viaduct[END_REF]. That works: [START_REF] Desprez | Seismic vulnerability assessment of a RC structure before and after FRP retrofitting[END_REF] Hence, the fiber can be considered as a 1D three-node bar as shown in Fig. 4.

⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ = ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ - - - ⎤ ⎦ ⎥ ⎥ ⎥ + yt
The axial displacement field u xyt (, , )

x along the fiber is thus interpolated by shape functions x N ( )

f such that: =- = * u xyt x t y x t x yt Nd N d Nd (, , ) () () () () () (, )
x u e θ ef f [START_REF] Mazars | Using multifiber beams to account for shear and torsion[END_REF] with

= x Nx Nx N x N () [ () () () ]
f fi fj fk [START_REF] Feng | An efficient fiber beam-column element considering flexure-shear interaction and anchorage bond-slip effect for cyclic analysis of RC structures[END_REF] These functions take the same expressions as * N x ( )

α θ in (8)
:

=- + = -+ = - Nx Nx N x () 1 , () , () fi x L x L fj x L x L fk x L x L 32 2
As a result, the multifiber beam element formulation is transformed into a combination of 1D three-node bar elements. On the other hand, Path (b) consists in interpolating the generalized displacement fields with the shape functions of the adopted beam formulation. Then, the axial displacement along each fiber is determined using the planar section hypothesis, Eq. ( 9), at the integration points. The displacement field is therefore no longer interpolated along each fiber but along the axis of the beam element as follows:

=- = - * * u xyt x t y x t x y x t Nd N d N N d (, , ) () () () () ( () () ) ()
x u e θ e uθ e [START_REF] Pinho | Dynamic collapse testing of a full-scale four storey RC frame[END_REF] As shown in Sections 2.2, 2.3 and 4 Path (a), Eq. ( 12), allows writing the axial displacement at the fiber level and thus to define Fig. 4. A three-node fiber.

the discontinuity kinematics and the corresponding enhancement functions. Path (b), Eq. ( 15), facilitates the writing of the variational formulation, the numerical implementation procedure and the construction of the element stiffness matrix.

Fiber kinematic enhancement

We propose hereafter a way to adapt the SDA for beams (see [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF][START_REF] Armero | An analysis of strain localization and wave propagation in plastic models of beams at failure[END_REF][START_REF] Ehrlich | Finite element methods for the analysis of softening plastic hinges in beams and frames[END_REF][START_REF] Armero | Numerical modeling of softening hinges in thin Euler-Bernoulli beams[END_REF][START_REF] Dujc | Multi-scale computational model for failure analysis of metal frames that includes softening and local buckling[END_REF][START_REF] Pham | Stress resultant model for ultimate load design of reinforced-concrete frames: combined axial force and bending moment[END_REF][START_REF] Jukić | Embedded discontinuity finite element formulation for failure analysis of planar reinforced concrete beams and frames[END_REF][START_REF] Bui | Enriched Timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF][START_REF] Bitar | A generalized timoshenko beam with embedded rotation discontinuity[END_REF][START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF]) for the higher order multifiber Timoshenko beam introduced in [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF].

Enhanced axial displacement

In order to account for several discontinuities within a fiber, the enhanced axial displacement field at the fiber level is written as follows (see Eq. ( 12)):

∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ +- = u xyt x yt H x ϕ x u t Nd ,, () , ( [ ()
() ] ( ) ) xf f d npg x d x d 1 () () d ( 16
)
with npg the number of integration points, Hx ( )

xd
the Heaviside function, ϕx ( )

d ()
the functions associated with each discontinuity (d) allowing to fulfill the compatibility conditions between the elements and u t ( )

x d ()
, the displacement discontinuities variables. Two Gauss integration points are sufficient to obtain an accurate integration of the axial forces at the fiber level since the interpolation functions of the strain field are linear. It is assumed that each integration point carries only one discontinuity (Fig. 5):

The coordinates of the integrations points and the associated integration weights in the real and reference element are:

         ⎜ ⎟⎜ ⎟ = ⎡ ⎣ - ⎤ ⎦ ⇒= ⎡ ⎣ ⎢ ⎛ ⎝ - ⎞ ⎠ ⎛ ⎝ + ⎞ ⎠ ⎤ ⎦ ⎥ xx 11 pg ref pg ele LL 1 3 1 3 Reference element 1 3 2 1 3 2 Real element ee ( 17
)
= w [START_REF] Desprez | Seismic vulnerability assessment of a RC structure before and after FRP retrofitting[END_REF] pg [START_REF] Marini | Analysis of reinforced concrete elements including shear effects[END_REF] In the proposed model in this paper, each fiber can carry one or two discontinuities. Also, it is a question of enhancing the integration point here rather than enhancing the element. This choice is made in order to respect the ultimate values of the stresses throughout the fiber, especially at the points that interest the calculation, i. e. at the integration points. More details will be addressed when elaborating the equilibrium equations.

Enhanced axial strain

The strain field at the fiber is easily derived from the axial displacement field as follows:

∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = ∂ ∂ + ⎛ ⎝ ∂ ∂ -+ ∂ ∂ ⎞ ⎠ = εx y t x xy t x ϕx ut x Hx u t Nd , , { () (, ) } { () } () { () } () xf f d npg d x d x x d 1 () () () d ( 19
)
To simplify the equations, the vector t u () x carrying the discontinuities by fiber is introduced (see Fig. 5):

= ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ t ut ut u () () () x x x (1) (2)
(20)

We also use the vector:

=+=+ δ xxxx G GGG () () () () rrrrx d (21)
divided into a regular (symbol •) and a singular part (symbol •) where x G () r is the vector of the enhancement functions associated with each discontinuity:

= x GxGx G () [ () () ] r rr (1) (2) (22)
with

=- ∂ ∂ = G x x ϕx d () { () } 1 ,2 r d d () () (23)
and δ xd the following Dirac functions vector:

= δ δxδx [ () () ] x xx (1) (2)
d dd [START_REF] Ian | Coupling 1d beams to 3d bodies[END_REF] The index r indicates that the interpolated discontinuities are real variables (see Section 2.3 for the virtual variables). Therefore, the enhanced strain filed takes the following form:

       ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ =+ + δ εx y t x y t x t t Bd Gu u ,, () ( ,) () ( ) ( ) xf f r x εx y t xx εx y t (,,) (,,) x d x ( 25
) with = ∂ ∂ == ⎡ ⎣ -+ -+ -⎤ ⎦ x x x Bx Bx B x BN () () [ () () () ]
In order to determine the enhancement functions G r d ( ) , known as the compatibility operators, two requirements should be met: (1) the introduction of displacement field discontinuities must not influence the nodal displacements so that the elements compatibility is satisfied. (2) to avoid stress locking phenomena G r d ( ) must be able to reproduce the zero hinge mode [START_REF] Ehrlich | Finite element methods for the analysis of softening plastic hinges in beams and frames[END_REF]; i.e. in the case of a fully opened discontinuity the cohesive stress at the discontinuity must vanish and the strain εx must tend to zero.

1st requirement: compatibility

The following three equations have to satisfied:

∑∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ =⇒ - =⇒ - = ∀∈ = == uy tu H ϕ u t ϕ u t dϕ 0, , ([ (0)
(0)] ( )) 0 ([0 (0)] ( )) 0 therefore [1, 2] (0) 0 xi d npg x d x d d npg d x d d 1 () () 1 () () () d (27)
∑∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ =⇒ - =⇒ - = ∀∈ = == uLy t u HL ϕ Lu t ϕ Lu t dϕ L ,, ( [ ()
() ] ( ) ) 0 ( [ 1 () ] ( ) ) 0 therefore [1, 2] ( ) 1
∑ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ =⇒ ⎡ ⎣ - ⎤ ⎦ =⇒ ∑ ⎡ ⎣ - ⎤ ⎦ == ∑ ⎡ ⎣ - ⎤ ⎦ == == = = = ( ) ( ) ( 1 () () 1 () () () d (28)
)

uy t u H ϕ u t ϕu t x x ϕu t x x xϕ xϕ ,, () () ( ) 0 1( ) ( ) 0 i f 0( ) ( ) 0 i f therefore for , ( ) 1
and for , ( ) 0 Following ( 27)-( 29), the enhancement functions ϕx ( )

x L k d npg x L d L x d d npg d L x d d d npg d L x d d d L d L 2 1 2 () 2
() 1 () 2
() 1 1 () 2
() 2 1 () 2 2
d ()
corresponding to the two discontinuities can be defined as:

=+ < => ϕx N x N x x ϕx N x x () () () s i n c e () () s i n c e fj fk L fj L (1) 1 2 (2) 2 2 e e ( 30
)
The derivatives of these functions give:

=- + = - =- = - + Gx B x B x Gx B x x ‾ () ( () () ) ‾ () ()
The values of Eq. ( 31) at the two integration points are summarized in the table below:

In Table 1, the values of G x ( )

r (1)
2 and G x ( )

r (2)
1 are not zero, indicating that the strain states at the integration points are coupled, i.e. the presence of a discontinuity at the first integration point influences the state of the discontinuity at the second integration point. This leads to a complex numerical implementation, see Section 4.

2nd requirement: zero hinge mode Eq. ( 31) should verify the zero hinge mode requirement. Three cases have to be checked: (i) only the first discontinuity is active, (ii) only the second discontinuity is active, (iii) both discontinuities are active.

The first case implies:

=- u uu x xk xi (1) and = u u xj xk . Therefore,      ̃⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ + =++ - ⎛ ⎝ ⎜ ⎜ + ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ - ⎞ ⎠ ⎟ ⎟ =+ + = = εx t x y t G x u t Bx u Bx u B x u Bx Bx u u Bx Bx B xu Bd , () , () () () () () () () (( ) ( ) ( ) ) 0 xf frx fi xi fj xj xk xk fj fk xj xi fi fj fk xi (1) (1) 0 (32)
The second case implies:

=- u uu x xj xk (2) and = u u xk xi . Therefore,      ̃⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ + =++ - ⎛ ⎝ ⎜ ⎜ - ⎞ ⎠ ⎟ ⎟ =+ + = = ε x t x y t G x u tB x u B x u B x u B x u u B xB xB x u Bd , () , () () () () () () ( () () () )
0 xf frx fi xi fj xj xk xk fj xj xk fi fj fk xi (2) (2) 0 (33)
The third case implies:

=- u uu x xj xk (2) and =- u uu x xk xi (1)
. Therefore,

     ̃⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ + + =++ - ⎛ ⎝ ⎜ ⎜ + ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ - ⎞ ⎠ ⎟ ⎟ - ⎛ ⎝ ⎜ ⎜ - ⎞ ⎠ ⎟ ⎟ =+ + = = ε x t x y t Gx ut Gx ut B x u B x u Bx u B x B x u u Bxu u Bx Bx B xu Bd , () , () () () () () () () () () () ( () () () ) 0 xf frx rx fi xi fj xj xk xk fj fk xk xi fj xj xk fi fj fk xi (1) (1) (2) (2) 0 (34)
Eqs. ( 32)- [START_REF] Bui | Enriched Timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF] show that the compatibility operators G r d ( ) satisfy the zero hinge mode requirement.

Variational formulation 2.3.1. Interpolation of the virtual fields

The upscript * • is adopted hereafter for the virtual variables. The virtual displacements are interpolated with the same shape functions as the real displacements. Using [START_REF] Pinho | Dynamic collapse testing of a full-scale four storey RC frame[END_REF], the axial virtual strains take the following expression:

=+ * * * εxyt xy t x t BdG u (, , ) (, ) () () () x f e v x ( 35
)
where G v is the enhancement function of the virtual discontinuities * u x , also known as the equilibrium operator, defined as the sum of a regular G v and a singular G v part:

=+=+ δ xxxx x G GGG () () () () () vvvvx d ( 36
)
The use of ( 15) enables to write the variational formulation using the virtual values of the nodal degrees of freedom * t d ( ) e , thus facilitating the numerical developments (see Section 4).

The enhancement functions G v associated with virtual discontinuities are not necessarily equal to G r associated with the real discontinuities (see for example [START_REF] Bitar | A generalized timoshenko beam with embedded rotation discontinuity[END_REF]). Actually, real discontinuities are interpolated following kinematic considerations (see Section 2.2) while virtual discontinuities are interpolated following static considerations (see Section 2.4). This is in order to satisfy the balance between the discontinuous and the continuous parts in the fiber, as well as the Patch test. Both interpolations in the same formulation are first proposed in [START_REF] Dvorkin | Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions[END_REF][START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: fundamentals[END_REF]. Following Timoshenko's theory, the virtual shear strain is constant and calculated as:

Table 1 Values of the kinematic operator x G () r at the two integration points.

x 1 x 2 Gx () r (1) ⎜ ⎟ -⎛ ⎝ -- ⎞ ⎠ 1 Le 2 2 3 ⎜ ⎟ ⎛ ⎝ -+ ⎞ ⎠ 1 Le 2 2 3
Gx ()

r (2) ⎜ ⎟ ⎛ ⎝ -+ ⎞ ⎠ 1 Le 2 2 3 ⎜ ⎟ ⎛ ⎝ -- ⎞ ⎠ 1 Le 2 2 3 ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ = ∂ ∂ -= ∂ ∂ -= *
∫ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ + ⎛ ⎝ ⎞ ⎠ * * W yt δ ε xyt σ xyt δ γ xyt τ xytd , ( ,,) ,, ( ,,) ,,
Ω int f x T x xy T xy f Ω f (38)
where σ x and τ xy are the axial and shear stresses at the fiber, * δε x and * δγ xy are respectively the axial and the shear virtual strains variations and Ω f is the fiber volume.

Introducing ( 35) and ( 37) in [START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF] results:

=+ * * W yt δ t yt δ yt yt dF u F (, ) () (, ) (, ) (, ) int f T int B f x T int G f , , (39)
where

∫ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ yt xy x σxyt τx y t d F BB ,[ (, ) () ] (, , ) (, , ) Ω int B f f γ T x xy f , Ω f (40)
and

∫ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ yt x σ xytd FG ,( ) , , Ω int G f v T x f , Ω f (41)

Beam level

The internal work at the multifiber beam element level is found simply by summing up the internal work of the fibers:

∑∑ = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ =+ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ * * W t W yt δ t t δ yt yt dF u F () , ( ) () ( , )
,

int e f nfib int f T int B e f n x T int G f , , f (42)
with n f the number of cracked fibers per beam element. A cracked fiber is a fiber with at least one discontinuity. Structural level Finally, the total internal work a the level of the structure is:

∑∑ ∑ == + ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ * * W tW t δ t t δ y ty t dF u F () () () () (
with n the number of elements with at least one cracked fiber.

The introduction of the last expression in the principle of virtual work gives:

∑∑ -= -+ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = * * Wt Wt δt t t δ y t y t dF F uF () () 0 () ( () () ) ( , ) , 0 int str ext str strT int B str ext str e n f n x T int G f , , f (44)
The equilibrium of the structure must be respected for all virtual displacements * δt d ()

str
as well as for any virtual displacement jump * δy t u (, )

x at the fiber level. This provides the following system of equations:

-= ∀∈ … ∀∈ … = tt en fn ey t FF F () () 0 {1, 2, , } et {1, 2, , ( )} : ( , ) 0 int B str ext str f int G f , , (45)
The first equation corresponds to the overall equilibrium of the structure. The other equations represent the local equilibrium at the level of each fiber with active discontinuities and can be solved locally at the fiber level.

yt F (, ) int G f ,
can be further developed using the decomposition [START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF] to give:

∫∫ ∫ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ + ⎞ ⎠ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ + ⎛ ⎝ ⎞ ⎠ = δσ ty t xσ x y t d x x σ x y t d xσ x y t d xy t FG G G () , () , , Ω () () , , Ω () , , Ω , , 0 int G f v Tf v T x Tf v Tf d , ΩΩ Ω ff d f (46) with ∫ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ δσ xσxytd x yt () , , Ω , , x Tf d Ω f d leading to: ∫ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ =-⎛ ⎝ ⎞ ⎠ =-⎛ ⎝ ⎞ ⎠ σ xσxytd x yt yt GC () , , Ω , , , v Tf d f Ω f (47) with yt C (, )
f the vector of cohesive stresses at discontinuities level belonging to the same fiber defined as:

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ yt Cy t Cy t C , (, ) (, ) f f f , (1)
,

(48)

2.4. Determination of x G ( ) v It remains to determine the functions = xG x G x G () [ () () ] v vv (1) (2)
. Let's assume a linear function for G x ( )

v d ()
of the form + ax b (a and b constants) and a linear expression for the fiber stress =+ σ xα x β () (α and β constants). a and b are determined by an identification procedure that ensures the following equality

∫ =- Gx σ x y t d σ x y t () (, , ) Ω ( , , ) v d Tf d Ω () f (47) for each discontinuity level = xx d . Indeed, ∫ ⎛ ⎝ ⎜ ⎜ + ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ + ⎞ ⎠ ⎟ ⎟ =- ⎛ ⎝ ⎜ ⎜ + ⎞ ⎠ ⎟ ⎟ ⇒ ⎧ ⎨ ⎩ =- =-+ ax b αx β d αx β ax bx Ω f d αβ identification LL d L L d Ω with
Therefore, the function G xx (, )

v d d ()
takes the following general form:

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎡ ⎣ ⎢ - ⎤ ⎦ ⎥ +- G xx LL xx L x L , 61 2 6 4 v d d ee d e d e () 23 2 (50)
Finally, for = xx dp g e l e (Eq. ( 17)), the equilibrium operators are: . The enhancement method with embedded discontinuities is similar to the incompatible mode method and as mentioned in the literature [START_REF] Taylor | The patch testa condition for assessing FEM convergence[END_REF] the equilibrium operator should verify the Patch test. This test, initially proposed by [START_REF] Irons | Experience with the patch test for convergence of finite elements, The mathematical foundations of the finite element method with applications to partial differential equations 557[END_REF],d e fines a convergence condition during mesh refinement (the elementary strains and stresses should stay constant). In other words, this test guarantees the ability to represent a constant state of stress per element [START_REF] Ibrahimbegović | Mécanique non linéaire des solides déformables: Formulation théorique et résolution par éléments finis[END_REF].

⎜ ⎟ ⎜ ⎟ =-⎛ ⎝ + ⎞ ⎠ =- = =- -⎛ ⎝ - ⎞ ⎠ == - Gx x Gx Gx Gx x Gx Gx () 1 3 w i t h ( ) a n d ( ) 0 () 1 3 w i t h () 0 a n d () v L L v L v v L L vv L (1) 23 1 (1) 1 2 (1) 2 (2
The patch test requires that the additional virtual work associated with the enhancement must be zero if the stress is constant along the element. This results the following condition:

∫ =- Gx d x () 1 L v d () e ( 52
)
This condition is easily verified using Eq. ( 50).

Constitutive laws

In an enhanced multifiber beam, each fiber corresponds to a specific material. Constitutive laws composed of a continuous and a cohesive part should be therefore defined at the fiber level. A damage mechanics law, an elasto-plastic law and the corresponding cohesive laws are detailed hereafter. This choice is particularly suitable for reinforced concrete structures (see Section 5), where damage mechanics is often adopted for concrete and plasticity for steel. The laws are presented in 1D as shear is considered decoupled and linear elastic (see Section 3.3).

Damage mechanics modelcohesive model

Continuous part

In damage mechanics constitutive laws, damage evolution is often driven by a strain threshold [START_REF] Mazars | A description of micro-and macroscale damage of concrete structures[END_REF][START_REF] Simo | Strain-and stress-based continuum damage models-I. Formulation[END_REF][START_REF] Simo | Strain-and stress-based continuum damage models-II. Computational aspects[END_REF][START_REF] Borderie | Phenomènes unilateraux dans un matériau endommageable: Modélisation et application à l'analyse de structures en béton[END_REF][START_REF] Jirásek | Nonlocal damage mechanics with application to concrete[END_REF][START_REF] Desprez | Damage model for FRP-confined concrete columns under cyclic loading[END_REF][START_REF] Mazars | A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings[END_REF]. A stress threshold is however chosen hereafter in order to simplify the numerical implementation. The damage model developed in [START_REF] Ibrahimbegović | Mécanique non linéaire des solides déformables: Formulation théorique et résolution par éléments finis[END_REF][START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF] is adopted. The origin of this model goes back to the work of [START_REF] Ortiz | A constitutive theory for the inelastic behavior of concrete[END_REF]. A scalar internal damage variable called compliance and noted D is in- troduced such that the 1D strain-stress relation takes the following form:

=∈ ⎡ ⎣ ⎢ ∞ ⎤ ⎦ ⎥ εD σ D E avec 1 , (53)
The stress-based damage threshold surface becomes:

=--⩽ ϕσ σ q () 0 by [START_REF] Jirásek | Nonlocal damage mechanics with application to concrete[END_REF] where σ by the elastic limit for a material that may have a different behavior in compression and tension (e.g. concrete). q the stress-like variable associated with the hardening mechanism controling the evolution of the damage threshold surface (50)

=- = ⎧ ⎨ ⎩ q Hξ H H H with for compression for tension bb b c b t (56)
H b the (positive) damage hardening modulus, ξ the strain-like internal variable that controls the hardening mechanism. For more details on the thermodynamic formulation of this model, the reader is referred to the work of [START_REF] Ibrahimbegović | Mécanique non linéaire des solides déformables: Formulation théorique et résolution par éléments finis[END_REF]. In the following, only those elements necessary for understanding the article are presented.

The free energy Ψ is written as follows:

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ =+ - Dξ ε D ε ξ H ξ Ψ, 1 2 1 2 b 1 (57)
The dissipated energy D b due to damage is (where the symbol •̇is the derivative with respect to time):

∫∫ =- = +> Dσ ε d σ D σ q ξ d ̇Ψ̇Ω 1 2 ̇̇Ω0 b ΩΩ (58)
As usual, the evolution of the internal variables are determined using the principle of maximum dissipation; among all permissible stresses, those that maximize dissipation are selected.

== == ∂ ∂ ∂ ∂ Dγ γ ξγ γ ̇̇σ ϕ σ sign σ σ ϕ q 1 () (59)
with γ ̇the Lagrange multiplier. The Kuhn-Tucker conditions and the consistency condition are:

⩾⩽ = = γ ϕ γϕ γϕ ̇00 ̇0 ̇̇0
The damage process is active when γ ̇is positive. Therefore, the value of ϕ ̇must necessarily be zero to comply with the consistency condition. This makes possible to deduce the expression of γ :

= + - - γ sign σ D ε DH ̇() ḃ1 1 (60)
Finally, the stress rate is written as function of the strain rate:

= ⎧ ⎨ ⎩ = > - + - - σ Dε γ εγ ̇̇if ̇0 ̇if ̇0 DH DH 1 b b 1 1 (61)

Cohesive part

A cohesive model is adopted to describe the behaviour at each discontinuity. The discontinuity activation criterion is again formulated in terms of stresses; a failure surface ϕ d ( ) is introduced at each integration point in order to check the activation and evolution of the discontinuity.

=- - = ϕt Ct σ qt d () () ( ()) 1 ,2 d d bu d () () () (62)
where Ct ()

d ()
is the cohesive stress determined by Eq. ( 47), written as: == Cy tσ x y tσ x y t ( , )( , , )( , , )

d dp g () (63)
Since the discontinuities d ( ) occur at the integration points pg () , the cohesive stresses Cy t (, )

d ()
should be equal to the continuous stresses σ xy t (, , ) pg . Therefore, the failure surfaces ϕ d ( ) at the integration points can be written as:

=- - ϕt σ x y t σ q x t ( ) ( , ,) ( ( ,)
)

pg pg bu pg () (64)
where σ bu is the ultimate stress, that may differ in compression and tension:

= ⎧ ⎨ ⎩ σ σ σ for compression for tension bu bu c bu t (65)
q is the stress-like variable associated with the cohesive model defined as: The cohesive behaviour is defined by the following relation: . Fig. 6 illustrates the cohesive model. The relation between the cohesive stress and the discontinuity can be also expressed as:

⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ =- ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ = ⎧ ⎨ ⎪ ⎩ ⎪ =- =- q xt S ξ
= u xt D xt C xt (, ) (, )(, )
=+ C x tS u x tσ (, ) (, )
f pg b pg bu (68)
The total free energy of a fiber with several discontinuities takes the following form:

=+ δ tt t Ψ Ψ() Ψ() () xd ( 69
)
with δ xd the Dirac function vector defined in Eq. ( 24) and t Ψ( )the free energy vector associated with the discontinuities defined as:

= ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ t t t Ψ() Ψ( ) Ψ( ) (1)
(

(70

)
where Ψ (1) and Ψ (2) are respectively the free energies associated with discontinuities 1 and 2. The expressions are determined according to the corresponding internal variables of the cohesive model:

=+ = - tD u S ξ d Ψ( ) 1 2 () () 1 2 () 1 , 2 dd d b d () () 1( ) 2 () 2 (71)
Using the strain Eq. ( 25) and the expression of total free energy (69), one can deduce the total dissipation at the fiber level as:

=+ Dy t D y t D y t (, ) (, ) (, ) tot ( 72
)
where D is the dissipation due to discontinuities:

∫ ∑ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = ⎡ ⎣ ⎢ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ - ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎤ ⎦ ⎥ =- ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = δ D yt σxytε yt d yt t yt Cu Ψ ,, , Ψ ̇,Ω ( , ) ( ) ̇, x d d ffT xx Ω 1 2 () f ( 73
)
Among the admissible internal variables that verify the failure criteria, those that maximize dissipation are selected. To do this, a Lagrange multiplier denoted γ is introduced. Therefore, the evolution equations of the internal variables are determined as: )

sign C C f () f f (74)
Finally, in order to determine the Lagrange multiplier γ ̇the consistency and Kuhn-Tucker conditions are used:

⩾⩽ = = γ ϕ γϕ γϕ ̇0, 0, ̇0, ̇̇0 (75)
If the softening mechanism is active, the multiplier γ ̇is strictly positive, which means that ϕ ̇must be zero to respect = γϕ ̇̇0. Using Eq. ( 62) and = ϕ ̇0 result to:

= γ S Cs i g nC ̇1 ̇( ) b f f (76)
The cohesive force is finally given by:

= ⎧ ⎨ ⎪ ⎩ ⎪ = +> < >= - C Du γ σS ξ s i g n C γ q σ γq σ pour ̇0 () ( ) p o u r ̇0a n d 0p o u r ̇0a n d f bu b f bu bu 1 ( 77
)
More information on the elaboration of these equations can be found in [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF][START_REF] Bitar | Modeling the failure of reinforced concrete structures using generalized and multifiber finite element beams[END_REF].

Elasto-plastic modelcohesive model

Continuous part

A classical 1D elasto-plastic model is briefly presented hereafter (for more details see [START_REF] Borja | Plasticity, modeling & computation[END_REF]). As usual, the partition of (regular) strains into an elastic and a plastic component is assumed:

=+ εxt ε xt ε xt (, ) (, ) (, ) ep (78)
The elastic and hardening behaviours are assumed to decoupled and therefore the free energy Ψ can be also decoupled into an elastic and a hardening term as follows: (79)

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ =+ εξ Eε Hξ Ψ, 1 2
with E s the Young modulus and H s the hardening modulus. The elastic threshold is expressed as:

=--ϕσ σ q () , sy (80) with σ sy the elastic limit and q the stress-variable hardening variable. ϕ is always negative when in elasticity. Once the elastic limit σ sy is reached, the material enters the plastic domain. As a result, the internal values associated with plasticity are activated and begin to evolve. The plastic dissipation is written as:

∫ =+ Dσ ε q ξ d [ ̇]Ω pT p Ω (81)
In order to determine the evolution equations of the internal variables, the classical principle of maximum dissipation using the Lagrange multiplier γ is applied. It is found that:

== == ∂ ∂ ∂ ∂ εγ s i g n σ γ ξγ γ ̇̇() ṗ ϕ σ ϕ q (82)
Using the consistency and the Kuhn-Tucker conditions, the rate of the Lagrange multiplier γ ̇is finally found:

= + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ γ Eε EH ̇φ σ s ϕ σ s ϕ σ ϕ q s ϕ q (83)

Cohesive part

The discontinuity represents a cohesive zone of zero thickness and appears at the integration points = xx dp g of the fiber when the stress exceeds a critical value σ u . The failure surface is written in terms of the cohesive force C as follows:

=-+ ⩽ ϕCU C σ S u (, ) ( ) 0 u (84)
The total dissipation energy is:

∫ =+= + + DD D σ ε q ξ d x C u [ ̇] ṫot T p T Ω (85)
As before, an optimization under constraint problem is solved resulting to:

= u γsign C ̇̇( ) (86)
In order to find the rate of Lagrange multiplier γ ̇, the Kuhn-Tucker and the consistency conditions are used:

⩾⩽ = = γ ϕ γϕ γϕ ̇0, 0, ̇0, ̇̇0 (87)
It is found that:

== γ S C S Csign C ̇1 ̇1 ̇() (88)

Shear behaviour

A multifiber Timoshenko beam can also account for the shear behavior. In the following and for simplicity reasons, shear is considered elastic and decoupled from the axial behavior. Therefore, 4. Computational procedure

⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ = ⎧ ⎨ ⎩ τxyt G yγxyt G y G G ,, ( ) ,

Linearization of the equilibrium equations

The first step toward the numerical implementation of the higher order enhanced multifiber beam is to linearize the equilibrium Eq. ( 45). This operation gives the following system:

⎜ ⎟ ⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ -- ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎫ ⎬ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪ = + + + yt yt yt KKK K KK KK KK d u u u FF A 00 00 •• • • •• • • 00 Δ Δ(, ) Δ(, ) • • Δ, () 0 0 • • 0 e n BB BG
with

∑ ⎜ ⎟ = ⎛ ⎝ ⎞ ⎠ = ty t KK () , BB f n BB 1 f f (91) ∫ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ∂ ∂ ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎞ ⎠ yt xy σ ε xyt xyd KB B ,( , ) , , , Ω BB fT x x ff Ω ff (92) ∫ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ∂ ∂ ⎛ ⎝ ⎞ ⎠ yt xy σ ε xyt xd KB G ,( , ) , , ( ) Ω BG fT x x r f Ω ff (93) ∫ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ∂ ∂ ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎞ ⎠ yt x σ ε xyt xyd KGB ,( ) , , , Ω GB v T x x ff Ω ff (94) ∫ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ∂ ∂ ⎛ ⎝ ⎞ ⎠ + ∂ ∂ yt x σ ε xyt xd A KG G C u ,( ) , , ( ) Ω GG v T x x r f f x f Ω ff f (95)
where A f is the fiber area and

∂ ∂ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ S S C u 0 0 f x C u C u C u C u (1) (2) f x f x f x f x (1), (1)
(1),

(2) (2), (1)
(2), [START_REF] Scordelis | Computer models for nonlinear analysis of reinforced and prestressed concrete structures[END_REF] (96) where S d () is the softening modulus for the cohesive models defined at integration points. The variables yt K (, )

GG ff and ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ yt u Δ ,
x f e take respectively a matrix and vector form when both discontinuities are simultaneously active on a fiber. The existence of two discontinuities per fiber necessitates to explain in details the numerical integration of the constitutive laws. The different integration paths are discussed hereafter.

Numerical integration of the elasto-plastic model and the cohesive model

An elastic prediction (trial) is first made. The cohesive trial stresses

+ C t ft r i a l 1 , (1),
and

+ C t ft r i a l 1 , (2),
are written as follows (the enhancement functions G x ()

v d ()
were determined in section (2.4) such as the cohesive stress C fd ,( ) at the discontinuity level is equal to the continuous model stress σ x () d , see Eq. ( 50)):

∫ =- = ++ + C G xσ xd σ x () () Ω ( ) t ft r i a l vt ft r i a l f t ft r i a l 1 ,(1), Ω (1) 1 , 1 , 1 f (97) ∫ =- = ++ + C G xσ xd σ x () () Ω ( ) t ft r i a l vt ft r i a l f t ft r i a l 1 ,(2), Ω (2) 1 , 1 , 2 f (98)
Assuming that each fiber is made of a single material and therefore = E xE () , the elastic predictor of the continuous stress

+ σ x ( ) t ft r i a l pg 1 , is: == + - ++ + σ xE εxE x y x ε Bd G u () () [(, ) ()
Developing expression (99) for the two integration points gives (calculation details are given in appendices (A.1) and (A.2)) provides: and therefore to check whether one or both discontinuities are activated:

= ++- = ++- + + + + C E x y Gx u x Gx u x ε x C E x y Gx u x Gx u x ε x Bd Bd ( (,) () () () () () ) ( (,) () () () () (
=- - + + + ϕσ Cq () t ft r i a l t ft r i a l u t ft r i a l 1 , 1 , 1 , ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ -- -- ⎤ ⎦ ⎥ ⎥ + + + + + + ϕ ϕ Cσ q Cσ q () () t ft r i a l t ft r i a l t ft r i a l u t ft r i a l t ft r i a l u t ft r i a l 1 , (1), 1 ,(1), 1 ,(1), 1 ,(2), 1 ,(2), 1 ,(2), (101)
Depending on the results, the following situations may occur, see Table 2: • Case 1: no active discontinuities in the fiber.

• Cases 2 and 3: one active discontinuity in the fiber. Calculation of the internal variables and stresses at the corresponding integration point.

• Case 4: two active discontinuities in the fiber. Calculation of the internal variables and stresses at two integration points.

Two active discontinuities per fiber

The Lagrange multipliers of the two failure surfaces are γx ( ) 1 and γx ( ) 2 . Considering that the cohesive models evolve independently at the two integration points x 1 and x 2 and using a backward Euler numerical integration scheme we get:

=+ =+ + + + ux u x γ x s i g n C ξx ξ x γ x () () Δ() ( ) () () Δ() tt t ft r i a l tt 11 1 1 1 , (1), 1 11 1 (102)
and

=+ =+ + + + ux u x γ x s i g n C ξx ξ x γ x () () Δ() ( ) () () Δ() tt t ft r i a l tt 12 2 2 1 ,(2), 1 22 2
(103)

In order to simplify the expressions, we set

= + sx s i g nC () ( ) t ft r i a l 1 1 , (1), and =
+ sx s i g nC () ( ) t ft r i a l 2 1 , (2),
. The new stresses at the two (104)

=+ + =+ + ++ ++ Cσ x E G x γ x s x E G x γ x s x Cσ x E G x γ x s x E G x γ x s x () () Δ()() () Δ()() () () Δ()() () Δ()() t f t
Eq. ( 104) can be written in a matrix form as follows:

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + + + + σx σx σx σx EG x s x EG x s x EG x s x EG x s x γx γx () () () () ()() ()() ()() ()() Δ() Δ()
⎡ ⎣ ⎢ ⎤ ⎦ ⎥ =- ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ - + + γx γx EG x S EG x s x s x E G xsxsx E G x S ϕx ϕx Δ() Δ() () ()()() ()()() () () () rr rr t trial t trial 1 2 (1) 1 (2) 121 (1) 212 (2) 2 1 1 1 1 2 (106)
According to Eq. ( 106), the calculation of Lagrange multipliers is interdependent. This means numerically that the calculation for both discontinuities will take place simultaneously and successively. The final step consists in updating the internal variables following Eqs. ( 102) and (103).

Complete failure at the two discontinuities

If both discontinuities are fully open, the two failure surfaces become zero ( -= + σ qx () 0

u t pg 1
) and no stress transfer occurs. Both cohesive stresses are therefore zero and this results to (see Appendix (A.6) for the details): (107)

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ =- ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ - + + γx γx EG x s x EG x s x EG x s x EG x s x σx σx Δ() Δ() ()() ()() ()() ()() () ()

Complete failure at one discontinuity

For the case of a complete failure at one discontinuity, the corresponding cohesive stress cancels out while the other discontinuity continues its opening process. Two equations have to be used, the first corresponding to the fully open discontinuity (

= + σ x () 0 t 11
) and the second resulting from the failure surface + ϕ t 1 of the still evolving discontinuity. This gives the following system:

⎧ ⎨ ⎩ = = + + σx ϕx () 0 () 0 t t 11 1 2 (108)
The first equation is similar to the first equation of the system in Eq. ( 107), while the second equation is the one in Eq. ( 106).We finally get (see Appendix (A.7) for the details): (109)

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ =- ⎡ ⎣ ⎢ - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ - + + γx γx EG x s x EG x s x E G xsxsx E G x S σx ϕx Δ() Δ() ()() ()() ()()() () () ()

Numerical integration of the damage mechanics model and the cohesive model

The numerical integration of the damage mechanics model associated with the cohesive model is presented hereafter. The evolution of the internal variables of each discontinuity is independent. Nevertheless, the calculation of strains and stresses as well as the corresponding Lagrange multipliers are coupled. Therefore, equations illustrating the evolution of the cohesive variables are hereafter adopted for each discontinuity.

The particularity of the damage model lies in the fact that the trial values of the displacement jumps ( + u x ()

t trial 1 1 and + u x () t trial 1
2 ) at the discontinuities at time step + t 1 are not the same as the ones (u x ( ) t 1 and u x ( ) t 2 ) of the previous time step t (this remark is important for a successful numerical implementation of the model). Indeed, (110)

== ≠ = == ≠ = ++ + + ++ + + u x D xC Dxσ x Dxσx ux u x D xC Dxσ x Dxσx ux () () () () () () () () () () () () () ()
t
The first step is the elastic prediction of stresses. We obtain the following system (see Appendix B for the details):

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ + + -- -- - - + - + σx σx Dx Gx D x Dx Gx D x Dx Gx D x Dx Gx D x Dx x y Dx x y Bd Bd () () () () () 1 () () () () () () () () () 1 () (,) () (,)
showing that elastic predictions of the stress states at the integration points are coupled.

The next step is to introduce the calculated stresses in the failure surfaces + ϕx ( )

t trial 1 1 and + ϕx ( ) t trial 1 2 : ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ + + + + ϕx ϕx σx σx σq x σq x () () () () () () . t trial t trial t trial t trial u t u t 1 1 1 2 1 1 1 2 1 2 (112)
and the stress vector of Eq. ( 111) in Eq. ( 112):

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ × ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + + -- -- - - + - + ϕx ϕx Dx Gx D x Dx Gx D x s x s x D x G x D x sx sx D x G x D x Dx x y s x Dx x y s x σx σx Bd Bd () () () () () 1 () () ()()() () () ()()() () () () 1 () (,) () () (,) () () ()
t trial t trial trt trt trt trt t f et t f et t t 1 1 1 2 1 1 (1) 11 1 1 (2) 12 1 2 1 2 (1) 21 1 2 1 2 (2) 22 1 1 11, 1 1 1 22 , 1 2 1 2 (113)
The different cases presented in Section 4.2 and Table 2 are found.

Two active discontinuities per fiber

The two failure surfaces ( + ϕx ( )

t trial 1 1 and + ϕx ( ) t trial 1
2 ) are positive. The internal variables associated with the continuous damage model are thus frozen while the internal variables associated with the cohesive model must be updated. The new stresses obtained at the integration points take the following form (see Appendices (B.5) and (B.6) for the calculation details):

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + + - + - + -- -- -- -- σx σx Dx x y Dx x y Dx Gx D x s xDx Gx D x s x Dx Gx D x s xDx Gx D x s x σx σx Dx Gx s xDx Gx s x Dx Gx s xDx Gx s x γx γx Bd Bd () () () (,) () (,) () () ()() () () ()() () () ()() () () ()() () () () ()() () ()() () ()() () ()()
In order to determine the Lagrange multipliers γx Δ ( ) 1 and γx Δ ( ) 2 , Eq. ( 114) are introduced to the failure surfaces:

= = ⇔ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ = ⎡ ⎣ ⎤ ⎦ + + + + + + ϕx ϕx σx σx σqx σqx ( ) 0 () 0 () () () () 0 0
t t t t u t u t 1 1 1 2 11 12 1 1 1 2 (115)
The Lagrange multipliers are calculated as follows (see Appendix B.1.1 for the details):

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ × ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ -- -- - -- -- + + γx γx D xG x S D xG xsxsx D x G x sx sx D x G x S Dx Gx D x Dx Gx D x s x s x D x G x D x sx sx D x G x D x ϕx ϕx Δ() Δ() () () () ()()() () ()()() () () () () () 1 () () ()()() () () ()()() () () () 1 () (
The internal variables associated with the cohesive models are finally updated.

Complete failure at the two discontinuities

When both discontinuities are fully open, we have:

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎤ ⎦ + + σx σx () () 0 0 t t 11 12 (117)
It can be found (see Appendix B.3 for the details):

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -= ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ × ⎛ ⎝ ⎜ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎞ ⎠ ⎟ -- -- - - + - + -- -- γx γx Dx Gx s xDx Gx s x Dx Gx s xDx Gx s x Dx x y Dx x y Dx Gx D x s xDx Gx D x s x Dx Gx D x s xDx Gx D x s x σxσx Bd Bd Δ() Δ() () ()() () ()() () ()() () ()() () (,) () (,) () () ()() () () ()() () () ()() () () ()() [()() ]
tr tr

tr tr t f et t f et trt trt trt trt tt 1 2 1 1 (1) 11 1 1 (2) 12 1 2 (1) 21 1 2 (2) 22 1 1 11, 1 1 22 , 1 1 1 (1) 11 1 1 1 (2) 12 2 1 2 (1) 21 1 1 2 (2) 22 2 12
(118)

Complete failure at a one discontinuity

As in Section 4.2.3, the following system is found:

⎧ ⎨ ⎩ = = + + σx ϕx () 0 () 0 t t 11 1 2 (119)
The corresponding Lagrange multipliers are:

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ =- ⎡ ⎣ ⎢ - ⎤ ⎦ ⎥ × ⎛ ⎝ ⎜ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎞ ⎠ ⎟ -- -- - - + - + -- -- γx γx Dx Gx s x Dx Gx s x D x G x sx sx D x G x S Dx x y Dx x y s x Dx Gx D x s x Dx Gx D x s x D x G x D x sx sx D x G x D x σx σx Bd Bd Δ() Δ() () ()() () ()() () ()()() () () () (,) () (,) () () () ()() () () ()() () () ()()() () () () 1 () ()
For the case when a constant stress/strain occurs in the fiber the choice made is to allow only one discontinuity.

Numerical applications

In order to validate the numerical implementation of the higher order enhanced Timoshenko multifiber beam and to study its performance, two numerical applications are presented: (1) a cantilever beam structure (2) a reinforced concrete framed structure tested by [START_REF] Cranston | Tests on reinforced concrete frames: Pinned portal frames[END_REF] and simulated in [START_REF] Saje | Large displacements and instability of beam-like structural systems, engineering structures under extreme conditions: multi-physics and multi-scale computer models in nonlinear analysis and optimal design[END_REF][START_REF] Jukić | Finite elements for modeling of localized failure in reinforced concrete[END_REF].

Cantilever beam structure

A cantilever beam structure of length = L 2.5 m submitted to a transverse displacement v y at its free end is studied hereafter, see Fig. 7. The aim of this example is to prove the numerical efficiency of the higher order multifiber beam (FCQ) and to compare its performance with first order multifiber beams (FLI) existing in the literature [START_REF] Guedes | A Fibre/Timoshenko beam element in Castem[END_REF].

A first case is considered where the fibers are not enhanced and follow a bilinear softening constitutive law (Fig. 8a). Symmetrical thresholds are considered in tension and compression. Then, a second case is studied with a linear elastic constitutive law coupled with a linear cohesive model (Fig. 8b). The material properties are summarized in Table 3.

The structure is discretized with different number (NE) of FLI or FCQ multifiber Timoshenko beams. Each section is discretized with 20 fibers. Results are presented in Fig. 9.

Fig. 9 illustrates the effect of mesh refinement on the global (bending momenttransversal imposed displacement) structural response considering or not fiber enhancement. Without enhancement, the results of the FLI and FCQ elements are similar and depend on the mesh size, see Fig. 9a andb. Fig. 9a andc show that the fiber enhancement improves the FLI performance. Nevertheless, mesh dependency is still significant in terms of softening response but also ultimate moment. Finally, Fig. 9b andd clearly show that the fiber enhancement of the FCQ Timoshenko beams makes the answer (almost) identical, independent of the number of elements.

Reinforced concrete frame

The non linear behavior of a reinforced concrete frame is studied hereafter. The frame has been tested by [START_REF] Cranston | Tests on reinforced concrete frames: Pinned portal frames[END_REF] and numerically modeled by several authors using classical finite elements [START_REF] Bažant | Softening in reinforced concrete beams and frames[END_REF][START_REF] Saje | Large displacements and instability of beam-like structural systems, engineering structures under extreme conditions: multi-physics and multi-scale computer models in nonlinear analysis and optimal design[END_REF] or FLI multifiber Timoshenko beam elements [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF]. The frame geometry is illustrated in Fig. 10. The geometry and the applied loading being symmetric, half of the frame is considered for the numerical model. All frame sections are of rectangular shape with a height = h 0.1524 m and width = b 0.1016 m. The top and bottom fiber concrete covering is 0.014 m. The steel reinforcement is shown in Fig. 10.

Concrete fibers are modeled with the damage model presented in Section 3.1.1. The steel fibers are modeled with the elasto-7. Cantilever beam structure submitted to a transverse displacement. 4 and5.

The frame is modeled using only 5 FCQ multifiber Timoshenko beam elements (NE = 5): two elements for the column and three Fig. 10. Reinforced concrete frame tested by Cranston [START_REF] Cranston | Tests on reinforced concrete frames: Pinned portal frames[END_REF].

elements for the beam. Each section is discretized with 20 fibers (18 for concrete and 2 for steel). The thickness of the steel fibers is calculated by respecting the steel ratio of the section, given in Fig. 10. Perfect bonding is considered between steel and concrete and 3D phenomenaas concrete confinement due to the stirrupsare not taken into account. In Fig. 11 (right), the notation S(eg , ) means the section S at the integration point g of element e. Also, f(i) denotes the fiber of number i. Due to the symmetric configuration of the frame structure (see Fig. 10, left), only half of it is simulated (see Fig. 10, right) and thus a vertical increasing displacement v y is imposed as shown in Fig. 11.

The global response (force versus applied vertical displacement) of the frame is illustrated in Fig. 12. Results of the NE = 5 FCQ multifiber Timoshenko beams model are compared with experimental results [START_REF] Cranston | Tests on reinforced concrete frames: Pinned portal frames[END_REF], classical finite element results [START_REF] Bažant | Softening in reinforced concrete beams and frames[END_REF] and the NE = 104 FLI multifiber Timoshenko beams model [START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF]. The FCQ model response is in good agreement with the experimental results Fig. 11. Reinforced concrete frame -Detail of the frame (left) and cross section (right) meshes.

[58] and the finite element numerical results [START_REF] Bažant | Softening in reinforced concrete beams and frames[END_REF]. The FLI model results are less accurate, despite the important number of beam elements adopted.

In order to analyze the local behaviour at the fiber level and to see how it affects the global response, the results of the FCQ model are isolated in Fig. 13. The force/displacement coordinates of the points A-G in Fig. 13 are summarized in Table 6. Each point corresponds to a local phenomenon that affects the global response of the structure. More specifically:

• Point A corresponds to the first change of the slope. It is related to the beginning of the hardening phase of the lower steel fiber in element 5 (see Fig. 14c).

• Point B corresponds to the same fiber when it reaches its ultimate stress (310 MPa), see Fig. 14c.

• Point C denotes the beginning of hardening of the upper steel fiber in Section 2 of the upper steel fiber (S(2,2)/f (2), Fig. 11)o f element 2 (Fig. 14a).

• Point D corresponds to the hardening of Section 1 (S(3,1), Fig. 11) of the upper steel fiber (f(2), Fig. 11) of element 3.

• Points E and F represent the beginning of hardening of the upper steel fibers f(2) of sections S(2,1) and S(3,2) of elements 2 and 3 respectively.

• Point G indicates the time step where the top-steel fibers of sections S(2,2) and S(3,1) of elements 2 and 3 reach their ultimate values ( = σ 310 su MPa).

++ + + ++ + C σxE εxε x Ex y G x uxG x uxε x E x y G xux G xux ε x Bd Bd () ( () () ) ( (,) () () () () () ) ( (,) () () () () ()
(A.1) (A.2)

= = - = ++- = ++- ++ + + ++ + C σxE εxε x Ex y G x uxG x uxε x E x y G xux G xux ε x Bd Bd () ( () () ) ( (,) () () () () () ) ( (,) () () () () (

A.1. Two active discontinuities per fiber

For the case of two active discontinuities per fiber, the new values of the cohesive stresses (104) are: (A.4)

= = - = ++- = ++= + + + ++ + + + ++ C σ xE ε xε x Ex y G x u xG x u xε x σx E G x u x E G x u x σx E G x γ x s x E G x γ x s x Bd () ( () () ) ( (,) () () () () () ) () () Δ() () Δ() () () Δ()() () Δ()()
(A.3) = = - = ++- = ++= + + + ++ + + + ++ C σ xE ε xε x Ex y G x u xG x u xε x σx E G x u x E G x u x σx E G x γ x s x E G x γ x s x Bd () ( () () ) ( (,) () () () () () ) () () Δ() () Δ() () () Δ()() () Δ()()
The Lagrange multipliers obtained in (106) are calculated as: (A.6)

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⇔ ⎡ ⎣ ⎤ ⎦ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⇔ ⎡ ⎣ ⎤ ⎦ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⇔ ⎡ ⎣ ⎤ ⎦ = ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ + ⎡ ⎣ ⎢ - - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + + + + + + + + + + + + ϕx ϕx σx σx σqx σqx σx s x σx s x E G x E G x sx sx E G xsxsx E G x γx γx σq x σq x Sγx Sγx σx s x σx s x σq x σq x EG x S EG x s x s x E G xsxsx E G x S γx γx ϕx ϕx EG x S EG x s x s x E G xsxsx E G x S γx γx () () () () () () 0 0 ()() ()() () ()()() ()()() () Δ() Δ() () () Δ() Δ() 0 0 ()() ()() () () () ()()() ()()() () Δ() Δ() 0 0 () () () ()()() ()()() () Δ() Δ()
⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎤ ⎦ ⇔ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎤ ⎦ + + + + σx σx σx σx EG x s x EG x s x EG x s x EG x s x γx γx () () 0 0 () () ()() ()() ()() ()() Δ() Δ() 0
From this system, one can easily deduce the two Lagrange multipliers (107).

A.3. Complete failure at a single discontinuity level

Developing the system of Eq. ( 108) gives: This equation can be reformulated as follows:

⎡ ⎣ ⎤ ⎦ = ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ + ⎡ ⎣ ⎢ - ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ + + σx ϕx EG x s x EG x s x E G xsxsx E G x S γx γx 0 0 () () ()() ()() ()()() () Δ() Δ()
-

+ = - - + - + - + Dx Gx D x σ x Dx Gx D x σ x Dx x y Bd ( () () () 1 ) () ( () () () ) () () (,) trt t trial trt t trial t f et 1 1 (1) 11 1 1 1 1 (2) 121 2 1 11, 1 (B.2)
The same procedure is applied for the other integration point. As in (B.1) the expression of

+ C t ft r i a l 1 , (2),
is calculated and we obtain: (B.4)

+- = - - + - + - + Dx Gx D x σ x Dx Gx D x σ x Dx x y Bd ( () () () ) () ( () () () 1 ) () () (,)
) () () 1 () () () () () () () () () 1 () () () (,) () (,)
The elastic stress trial values are deduced as in (111).

B.1. Two active discontinuities per fiber

For the case of both discontinuities active per fiber, the new values of the damage cohesive stresses (114) are: The latter equation is reformulated as follows: The Lagrange multipliers are obtained as in (116).

== + + =+ + + + + - + - ++ + - + σx Dx ε x Dx x y Gx u x Gx u x Dx x y Gx D x σ x γ xs x Gx D x σ x γ xs x Bd Bd () () () () ( (,) () () () () ) () ( (,) () [ () () Δ() ]() () [ () () Δ() ]() )
== + + =+ + + + + - + - ++ + - + σx Dx ε x Dx x y Gx u x Gx u x Dx x y Gx D x σ x γ xs x Gx D x σ x γ xs x Bd Bd () () () () ( (,) () () () () ) () ( (,) () [ () () Δ() ]() () [ () () Δ() ]() )

B.2. Complete failure at two discontinuities

The system of Eq. (117) gives: The Lagrange multipliers are finally deduced as in (120).

Fig. 2 .

 2 Fig. 2. Modified Timoshenko FCQ beam.

Fig. 3 .

 3 Fig. 3. Two multifiber beam formulations: (a) fiber level (Path a), (b) beam level axis (Path b).

Fig. 5 .

 5 Fig. 5. (a) Three-node fiber. (b) Three-node fiber: two integration points, one discontinuity by integration point.

  additional verification is required to validate the choice of the equilibrium operator G x ( )

  softening (negative) modulus, G b c and G b t the compression and tensile failure energies and ξ the strain-like internal variable whose role is to control softening.

  with D the compliance variable associated with cohesive model. This variable increases progressively with the discontinuity =

Fig. 6 .

 6 Fig. 6. Cohesive model.

Fig. 8 .

 8 Fig. 8. Constitutive laws (a) case without enhancement (b) case with enhancement.

Fig. 9 .

 9 Fig. 9. Cantilever beam structure -Bending moment versus transversal imposed displacement as a function of the number (NE) and type (FLI or FCQ) of multifiber Timoshenko beam elements.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Reinforced concrete frame -Force versus vertical imposed displacement: experimental, finite elements (FE), FLI and FCQ model results.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Reinforced concrete frame -Evolution of stress in the fibers, FCQ model.

(B. 3 )

 3 Eqs. (B.2) and (B.3) can be written as:

t trial 1 1

 1 ) derived from the elastic prediction multiplied by a factor of -( 1). This makes possible to get a simpler expression:

  Principle of virtual work 2.3.2.1. Fiber level. The internal work of each fiber is:

	γx y t , , xy	* γx t , xy	x	* Vx () Θ() * x yz	x	* xt Nd () () v ee * * t xt N d Bd () () () γ Θ e	(37)
	2.3.2.

  ]

	t	ft r i a l 1 ,	pg	t	ft r i a l 1 ,	pg	f pg	e t ,1	r pg	t	t	p

  ) )

	t	ft r i a l 1 ,(1),	f	et 1, 1	r (1)	t 11	r (2)	t 12	t	p	1
	t	ft r i a l 1 ,(2),	f	et 2, 1	r (1)	t 21	r (2)	t 22	t	p	2		(100)
	The above stresses (100) are used to calculate the trial failure surfaces + ϕ t ft r i a l 1 ,(1),	and + ϕ t ft r i a l 1 ,(2),	of the vector + ϕ t ft r i a l 1 ,

Table 2

 2 Different cases depending on the trial failure surfaces.

	Case 1	ϕ t	+ ft r i a l 1 ,(1),	<	0	and	ϕ t	+ ft r i a l 1 ,(2),	<	0
	Case 2	ϕ t	+ ft r i a l 1 ,(1),	>	0	and	ϕ t	+ ft r i a l 1 ,(2),	<	0
	Case 3	ϕ t	+ ft r i a l 1 ,(1),	<	0	and	ϕ t	+ ft r i a l 1 ,(2),	>	0
	Case 4	ϕ t	+ ft r i a l 1 ,(1),	>	0	and	ϕ t	+ ft r i a l 1 ,(2),	>	0

Table 3

 3 Material properties.

	Properties	Symbol	Value
	Young modulus	E	× 4 10 7 kN m -2
	Poisson coefficient	ν	0.2
	Softening modulus	S	-× 81 0 5 kN m -3
	Ultimate elastic stress	σ u	5000 kN m -2

Table 4

 4 Reinforced concrete frame -Concrete properties.

	Properties	Symbol	Value
	Young modulus	E b	× 31.5 10 6 kN m -2
	Poisson's coefficient	ν b	0.2
	Ultimate stress in traction	σ btu	950 kN m -2
	Elastic stress in compression	σ bcy	14,600 kN m -2
	Ultimate stress in compression	σ bcu	36,500 kN m -2
	Hardening modulus in compression	H bc	× 29 10 6 kN m -2
	Softening modulus in compression	S bc	-× 40 10 6 kN m -3
	Softening modulus in tension	S bt	-× 11 10 6 kN m -3

Table 5

 5 Reinforced concrete frame -Steel properties.

	Properties	Symbol	Value
	Young Modulus	E a	× 200 10 6 kNm -2
	Poisson's coefficient	ν a	0.3
	Elastic stress limit	σ ay	293,000 kN m -2
	Ultimate stress	σ au	310,000 kN m -2
	Hardening modulus	H a	× 5.02 10 5 kN m -2
	Softening modulus	S a	-× 42 10 6 kNm -3

Table 6

 6 Reinforced concrete frame -The coordinates of the green points on the Fig.13.

	Point	0	A	B	C	D	E	F	G
	v y (m)	0	0.0084	0.01504	0.01904	0.02208	0.0272	0.0432	0.0546
	F y (kN)	0	15.93	20.49	21.71	22.51	21.98	19.17	16.12

  )

	t	ft r i a l 1 ,	t trial	t trial	t	p	f	et		rt trial	rt trial	t	p
			f	et	r	t	r	t	t	p

  ) )

	t	ft r i a l 1 ,(2),	t trial 1	2	t trial 1	t 22 p	f	2	et , 1			rt trial (1) 2 1	1	rt trial (2) 2 1	t 22 p
			f		et 2, 1	r (1)	t 21	r (2)	t 22	t	p	2

  When the two discontinuities are fully open, there is no stress transfer and therefore the stresses are zero:

	A.2. Complete failure at two discontinuities
	t	1 1	t	11	u	t	1 1
	t	1 2	t	12	u	t	1 2
			t trial 1	11				rr (1) 1 (2)	121	1	u	t	1	1
			t trial 1	22		rr (1) 212 (2)	2	2	u	t	2	2
			t trial 1	11	u		t	1		rr (1) 1 (2)	121	1
			t trial 1	22	u		t	2	rr (1) 212 (2)	2	2
			t trial 1	1	rr (1) 1 (2)	121	1
			t trial 1	2	rr (1) 212 (2)	2	2	(A.5)
										22

  0

	t	11
	t	12
		t trial 1	1	rr (1) 11 (2)	12	1
		t trial 1	2	rr (1) 21 (2)	22	2

  The use of this notation in the failure surfaces equations results:

	-= σq x () u t 11 σ x () t
	-= σq x () u t 22 σ x () t												(B.7)
	⎡ ⎣ ⎢	-Dx () (,) x y Bd t f 1 11, 1 + et -+ Dx x y Bd () (,) t f et 1 22 , 1	s x () 1 s x () 2	⎤ ⎦ ⎥ +	⎡ ⎣ ⎢	---Dx Gx D x Dx Gx D x s x s x -() () () 1 () () ()()() trt trt 1 1 (1) 11 1 1 (2) 12 2 1 --D x G x D x sx sx D x G x D x () () ()()() () () () 1 trt trt 1 2 (1) 21 1 2 1 2 (2) 22	⎤ ⎦ ⎥	⎡ ⎣ ⎢	σx () t 1 σx () t 2	⎤ ⎦ ⎥
		+	⎡ ⎣ ⎢	---D xG x S D xG xsxsx -() () () ()()() tr tr 1 1 (1) 1 1 1 (2) 121 --D x G x sx sx D x G x S () ()()() () () tr tr 1 2 (1) 212 1 2 (2) 2	⎤ ⎦ ⎥	⎡ ⎣ ⎢	γx Δ() 1 γx Δ() 2	⎤ ⎦ ⎥	= ⎡ ⎣	0 0	⎤ ⎦
	t	12		t	t 1 21 2	t	f 1 22	et , 1	r (1)	t 21 1	r (2)	t 21 2
				t	f 1 22	et , 1	r (1)	tt 21	1			1	1	r (2)	tt 22	2	2	2	(B.6)
	B.1.1. Calculation of the Lagrange multipliers
	The calculation of the Lagrange multipliers γx Δ ( ) 1 and γx Δ ( ) 2 can be done using (114) and (115):
	⎡ ⎣ ⎢	-Dx () (,) x y Bd t f 1 11, 1 + et -+ Dx x y Bd () (,) t f et 1 22 , 1	s x () 1 s x () 2	⎤ ⎦ ⎥ +	⎡ ⎣ ⎢	--Dx Gx D x Dx Gx D x s x s x () () () () () ()()() trt trt 1 1 (1) 11 1 1 (2) 12 2 1 --D x G x D x sx sx D x G x D x () () ()()() () () () trt trt 1 2 (1) 21 1 2 1 2 (2) 22	⎤ ⎦ ⎥	⎡ ⎣ ⎢	σx () t 1 σx () t 2	⎤ ⎦ ⎥
		+	⎡ ⎣ ⎢	--Dx Gx Dx Gx s x s x () () () ()()() tr tr 1 1 (1) 1 1 1 (2) 121 --D x G x sx sx D x G x () ()()() () () tr tr 1 2 (1) 212 1 2 (2) 2	⎤ ⎦ ⎥	⎡ ⎣ ⎢	γx Δ() 1 γx Δ() 2	⎤ ⎦ ⎥	-⎡ ⎣ ⎢	--σq x () u t 1 σq x () u t 2	⎤ ⎦ ⎥	-⎡ ⎣ ⎢	Sγx Δ() 1 Sγx Δ() 2	⎤ ⎦ ⎥	= ⎡ ⎣	0 0	⎤ ⎦
	Furthermore, since we have:

  From Eq. (B.8) it can be noted that the first two terms are the expressions of the failure surface ( + ϕx ( )

	trt 1 1 (1) 11		trt 1 1 (2) 12	2	1	1	t	f 1 11, 1 et	1	t	1
	trt 1 2 (1) 21	1	2	trt 1 2 (2) 22				t	f 1 22	et , 1	2	t	2
		() trt () () 1 1 1 (1) 11		() trt () ()()() 1 1 (2) 12 2 1	1	() tr 1 1 (1)	() 1	() tr 1 1 (2)	()()() 121
	() trt () ()()() 1 2 (1) 21 1 2	() trt () () 1 1 2 (2) 22		() tr 1 2 (1)	()()() 212	() tr 1 2 (2)	() 2
	Δ() 1	0
	Δ() 2	0									(B.8)
											t trial 1	1 and + ϕx ( )

  Therefore, the system (118) is deduced.B.3. Complete failure at a single discontinuity levelThe development of the first equation of the system (119) gives:

	⎡ ⎣ ⎢	-Dx () (,) x y Bd t f 1 11, 1 + et -+ Dx x y s x Bd () (,) () t f et 1 22 , 1 2	⎤ ⎦ ⎥ +	⎡ ⎣ ⎢	-Dx Gx D x s x --Dx Gx D x s x () () ()() () () ()() trt trt 1 1 (1) 11 1 1 1 (2) 12 2 --D x G x D x sx sx D x G x D x () () ()()() () () () 1 trt trt 1 2 (1) 21 1 2 1 2 (2) 22	⎤ ⎦ ⎥	⎡ ⎣ ⎢	σx () t 1 σx () t 2	⎤ ⎦ ⎥
		+	⎡ ⎣ ⎢	-Dx Gx s x --Dx Gx s x () ()() () ()() tr tr 1 1 (1) 11 1 1 (2) 12 --D x G x sx sx D x G x S () ()()() () () tr tr 1 2 (1) 212 1 2 (2) 2	⎤ ⎦ ⎥	⎡ ⎣ ⎢	γx Δ() 1 γx Δ() 2	⎤ ⎦ ⎥	= ⎡ ⎣	0 0	⎤ ⎦	(B.11)
												1 1	(1)	11	1	1 1	(2)	12	2	1
		1 22		, 1				1 2	(1)	21	1	1 2	(2)	22	2	2
					1 1	(1)	11		1 1	(2)	12	1
					1 2	(1)	21		1 2	(2)	22	2
	⎜ ⎛ ⎝ x y Bd ⎟ ⎞ ⎠ , f 1 11, 1 -+ Dx () t et	+	[	--Dx Gx D x s xDx Gx D x s x () () ()() () () ()() ] trt trt 1 1 (1) 11 1 1 1 (2) 12 2	⎡ ⎣ ⎢	σx σx () () t t 1 2	⎤ ⎦ ⎥
		+	[	--Dx Gx s xDx Gx s x () ()() () ()() ] tr tr 1 1 (1) 11 1 1 (2) 12	⎡ ⎣ ⎢	γx γx Δ() Δ() 1 2	⎤ ⎦ ⎥	=	0	(B.9)
	The second equation of the system (119) becomes:
	⎜ ⎛ ⎝ x y Bd ⎟ ⎞ ⎠ , f 1 21, 1 -+ Dx () t et	+ () [ s x 2	-⎡ ⎣ ⎢ () () 1 ] D x G x D x sx sx D x G x D x --() () ()()() () trt trt 1 2 (1) 21 1 2 1 2 (2) 22	σx σx () () t t 1 2	⎤ ⎦ ⎥
		+	[	--D x G x sx sx D x G x () ()()() () () tr tr 1 2 (1) 212 1 2 (2) 2	-⎡ ⎣ ⎢ S ]	γx γx Δ() Δ() 1 2	⎤ ⎦ ⎥	=	0	(B.10)
	The two equations developed above are written as:

Acknowledgments

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).

Following the previous analysis, the behaviour of the steel fibers dominates the global reinforced concrete frame response. Nevertheless, discontinuities in concrete fibres start developing from the very beginning. Indeed, Fig. 15 illustrates the number of fibers with discontinuities developed in each element (ELE) as a function of the imposed vertical displacement. This numerical study proves the ability of the higher order multifiber FCQ Timoshenko beam to reproduce the non linear behavior of realistic reinforced concrete structures and to follow the evolution of non-linearities both at global and local levels, reason for considering it as a high valuable multi-scale numerical tool.

Conclusion

A novel high order enhanced multifiber Timoshenko beam is introduced in this article. The element is based on the Timoshenko finite element beam formulation of [START_REF] Caillerie | A Timoshenko finite element straight beam with internal degrees of freedom[END_REF]. The novelty of the enhanced beam lies in its ability to reproduce several discontinuities along the fibers ensuring compliance with the stress limit values. The variational formulation, the constitutive models and several computational issues are covered. Numerical examples are given to validate the model and to study its performance.

As a perspective of this work, one can mention the enhancement of the transverse component of the displacement field at the fiber level in order to study torsion or shear problems. Coupling between axial and transverse discontinuities can be also useful to simulate two-dimensional softening behaviour at the fiber level. Finally, extension of the model for cyclic and dynamic loadings requires to take into account crack opening and closure.