A generalized Timoshenko beam with embedded rotation discontinuity coupled with a 3D macroelement to assess the vulnerability of reinforced concrete frame structures

A generalized finite element beam with an embedded rotation discontinuity coupled with a 3D macroelement is proposed to assess, till complete failure (no stress transfer), the vulnerability of symmetrically reinforced concrete frame structures subjected to static (monotonic, cyclic) or dynamic loading. The beam follows the Timoshenko beam theory and its sectional behavior is described in terms of generalized forces and generalized strains. The beam response up to the peak is described by a macroelement, based on plasticity theory, that adopts a 3D failure criterion expressed in terms of axial force, shear force and bending moment. The Embedded Finite Element Method is then adopted to reproduce bending dominated failure, with a global cohesive model that links the cohesive moment to a rotational jump. The formulation allows for remedy of localization phenomena and significant reduction of the necessary computational time. The performance of the proposed simplified strategy is illustrated by comparison with experimental results.


Introduction
Numerical simulation is nowadays a powerful tool which allows conducting numerical experiments and predicting the behavior of reinforced concrete (RC) structures.This is not only essential for the design of a new structure, but also for evaluating the response of an existing one subjected to static (monotonic [1] or cyclic [2]) or dynamic excitations (e.g.earthquakes, see for instance [3]).Given the significant developments in the constitutive description of concrete and steel and the numerical methods over the last five decades, the Finite Element Method (FEM) and enriched formulations such as the Embedded Finite Element Method (E-FEM) [4] are capable of assessing the vulnerability of RC structures till failure.Numerical simulations may however exhibit prohibitive calculation times, and therefore they are not always suitable for industrial design applications.To overcome this drawback a simplified approach is proposed hereafter based on (1) the original macroelement concept introduced by Nova et Montrasio [5] and (2) the kinematically enhanced formulation introduced by Armero and Ehrlich [6][7][8] for generalized finite element (FE) beams.
Introduced in geomechanics by Nova and Montrasio [5] for soil structure interaction problems, the macroelement considers the global behavior of a foundation and of the soil volume interacting with it merged into a single, 3D stress-resultant constitutive equation relating the forces/moments on the foundation with the corresponding displacements/rotations. Subsequent developments of macroelements for soil structure interaction problems can for example be found in the work of [9][10][11] for shallow foundations, [12][13][14][15] for piles and [16] for caisson foundations respectively, following either the elasto-plasticity or the hypoplasticity computational framework.In structural engineering, the macroelement concept has been applied to study the nonlinear response of RC elements in terms of forces/moments -displacements/rotations [17,18] The word macroelement is sometimes used in the literature to describe uniaxial (uncoupled) nonlinear springs [19,20].In such cases, the behavior of each uniaxial spring is independent of the other springs.In the following however, the word macroelement is used in accordance with the original definition of Nova and Montrasio [5], i.e. a (coupled) 3D stress-resultant constitutive equation.
In order to simulate the behavior up to complete failure (no stress transfer), kinematically enhanced formulations coupled with cohesive models were first introduced by Armero and Ehrlich [6][7][8] for generalized Euler-Bernoulli and Timoshenko beams, initially for steel structures.Their work is based on E-FEM, initially developed by [4] and inspired by the Mixed Assumed Strain Approach [21].The main idea is that a crack can pass through the FE, while the enhanced variables are defined at the element and not at the nodal level.More specifically: Pham et al. [22,23] used the FE Timoshenko beam presented in Section 2 for RC structures.The authors adopted uncoupled (1D) stress-resultant constitutive models for the continuous part and a linear generalized cohesive model describing the relation between moment and rotation jump to represent the bending behavior until complete failure.The parameter identification of the generalized bending cohesive law was carried out via kinematically enhanced multi-fiber Timoshenko FE beams, enriched by axial displacement discontinuity at the fiber level.The longitudinal dimension of the beams was relatively small with respect to the transversal dimensions in order to represent a sectional analysis.In their formulation, different constitutive laws (continuum and cohesive, the latter given in terms of a stress-axial displacement jump) at a local level were assigned to each fiber to represent the 1D concrete and steel behaviors.
Jukić et al. [24] proposed a stress-resultant Euler-Bernoulli beam with linear shape functions for the axial component and high order shape functions for the transversal and rotational components for the continuous part of the finite element formulation.The enhanced formulation integrated the discontinuity of the rotational component.The constitutive description at the discontinuity was given by a linear cohesive model in terms of moment-rotation jump.Jukić [25] developed a multi-fiber Euler-Bernoulli beam and with co-authors [26] a multi-fiber Timoshenko beam for the analysis of RC beams and frames until failure.In both works, the axial kinematical field of the fibers was enriched by discontinuities.The continuous constitutive description of the concrete layer was given by a damage model while for the discontinuity a softening-damage traction-separation cohesive model was adopted.An elasto-plastic model and a softening-plasticity traction-separation cohesive law were used to describe the continuous and the discontinuous part of the steel.
Bui et al. [27] investigated the response of RC beam-columns by considering two modes of failure: bending and shear.A generalized Timoshenko beam element enriched by two discontinuity variables (transversal displacement and rotation jumps) was used.The results of the numerical simulations of a four-point bending test were presented by considering only shear failure, only bending failure and both shear and bending failures The authors showed that when considering only shear failure, with the bending response remaining elastic, the overall resistance was highly overestimated and the response was brittle.When considering only bending failure or both shear and bending failures, the obtained results were quite similar, more realistic and the response was more ductile.The authors explained that the type of failure depended on the peak resistances of the generalized shear and bending models.The dominant failure mechanism is therefore mainly driven by the geometrical and mechanical properties of the RC element.This conclusion is the basis of the assumptions of the enhanced formulation presented in Section 4.1.
Bitar and co-authors [28][29][30] used a high-order Timoshenko beam finite element as the basis of the enhanced formulation.This beam was initially developed by [31] and presents the advantage that it is free of shear locking and that only one element provides the exact displacements at the nodes for whatever loading in linear elasticity.In the first part of their work [28,29], a generalized beam was considered and uncoupled 1D stress-resultant models for the continuous part, while a rotational jump was introduced to account for bending failure.In [28,30], the same authors presented the enhanced formulation of a multi-fiber beam.The fibers were enriched at the Gauss points by axial discontinuities whose behavior was described by a linear cohesive law linking the axial stress with the displacement jump.The continuous part of the uniaxial behavior was expressed by a damage law for the concrete fibers and an elasto-plastic law for the steel fibers, while the shear behavior was assumed elastic.The authors finally compared the results of their enhanced high-order Timoshenko beam with the results of the FE Timoshenko beam presented in Section 2.
In the following, we propose to couple the macroelement approach with the kinematically enhanced formulation for a Timoshenko beam in order to assess, with a simplified modeling strategy, the vulnerability of symmetrically reinforced concrete frame structures subjected to static (monotonic, cyclic) or dynamic loading till complete failure (no stress transfer).The coupling of the two approaches presents the following advantages: on one hand, the use of the macroelement permits to take into account the various multidirectional couplings up to the peak.On the other hand, the enhanced formulation makes possible to simulate the evolution of the plastic hinges even for extreme loadings, while providing a way to regularize the finite element results.More specifically, the global behavior of a RC beam section is lumped into a single, integral, constitutive equation linking the evolution of the stress-resultant forces on the sectional gravity center to the corresponding strain histories.The proposed macroelement (stressresultant model) is then implemented in a Timoshenko FE beam where the softening behavior till complete failure (no stress transfer) is reproduced using a global cohesive model, within the framework of the E-FEM.The global cohesive model describes the response in terms of a cohesive moment -rotational jump, assuming that the RC structures under study are dominated by a bending failure mode.
The outline of the article is the following: Section 2 briefly presents the adopted Timoshenko FE beam formulation and Section 3 the novel macroelement for RC beams and columns sections, based on the 3D failure surface introduced in [32,33] and the plasticity theory.The procedure to identify its parameters is carefully detailed.Section 4 describes the way to simulate bending dominated failure using an embedded rotation discontinuity and a global cohesive model.Finally, the performance of the proposed simplified strategy is illustrated in Section 5, using the experimental results of a two-story RC frame subjected to monotonic and cyclic loads, (Section 5.1), and a RC simply supported beam subjected to an earthquake, (Section 5.2).
In the following, the symbol Ā defines the vector , Ā the matrix , Â the variable  in a standardized (dimensionless) space,

Ȧ
the derivative of  with respect to time,  ′ the derivative of  with respect to the coordinate  and  the symbol of the transpose.In Section 4, Ā is a variable assigned to a continuous field while Ā a variable assigned to a discontinuous field.

A Timoshenko beam finite element
The reference configuration of a 2D beam is shown in Fig. 1 (left).According to the Timoshenko theory [34,35] the displacement vector field ū  = (  ,   ) of a point  ≡ (, ) of the beam is expressed in terms of the displacement of the gravity center  of the section Ū   = (  ,   ) and the rotation   of the section as follows: Using the small strain assumption, the strain components are given by: Consider now a 2D FE beam of length   , with two nodes  and  and three degrees of freedom per node: axial displacement (  and   ), transversal displacement (  and   ) and rotation (  and   ), see Fig. 1 (right).The generalized displacement vector Ū is function of the nodal displacement vector Ū as follows [34,35]: (5) with N the matrix of the linear shape functions,  ∈ [0,   ], while the upper indexes , ,  denote the axial, transversal and rotational components respectively.By derivation of the displacement field and according to Eq. (2) the generalized strain vector ε is obtained as follows: with   ,   and   the axial strain, the shear strain and the curvature at the beam axis that passes through the point  respectively.Considering that   can be written as  and that Eq. ( 7) leads to shear locking, [36] proposed to solve the problem by neglecting the linear terms of Eq. ( 8).The strains are therefore approximately calculated with the simplified matrix B below (see also [34,35]): Using the following form of the principle of virtual work [34,35] ∫ together with the constitutive description of the section in terms of generalized forces (see Section 3.2), the following expression is obtained: Combining expressions (6) and (11), the element stiffness matrix is given by: with K the stiffness matrix of the section (see Section 3.2).As usual for FE beams, the element mass matrix is with M the mass matrix of the section given as: where  the material density,  the area and  the quadratic moment of the section.
The element nodal forces due to the stress resultants at the beam section (gauss point) are finally calculated as with F the vector of forces in the section (see Section 3.2).The numerical integration of the above integrals is done considering one (1) Gauss integration point in the middle of the beam [34].A study on the performance of the FE Timoshenko beam introduced above can be found in [37].
Remark.The displacement based Timoshenko beam finite element of [34,35] is adopted in this article due to its numerical performance (no shear locking, [37]) and its simplicity (linear shape functions, Eqs. ( 5)) that makes relatively easy the introduction of the enhanced formulation of Section 4.1 to model up to failure.The enhanced formulation can be also used with higher-order displacement based Timoshenko beam finite elements, as the one introduced by [31], where only one element in linear elasticity suffices to recover the exact nodal displacements for whatever loading (see [28] or [30] in a multifiber beam context).One can finally imagine the extension of the enhanced formulation to a flexibility based Timoshenko beam model, as the one introduced by [38].

A macroelement for RC beam and columns sections
A novel macroelement for RC beams and RC columns is introduced hereafter.The proposed macroelement is a 3D stress-resultant constitutive equation linking the evolution of the resultant forces/moments to the corresponding generalized strains at the level of a Timoshenko generalized beam section.The constitutive description of the RC section is built within the classical elasto-plastic computational framework while the adopted failure criterion, developed for symmetrically reinforced concrete square sections, was introduced by Doulgeroglou et al. [32].The analytical convex polynomial expression of the failure surface is of degree six (6), links the axial force, shear force, and bending moment and was numerically reproduced via 3D non-linear finite element simulations [32].As usual in the literature for macroelements, equations are written in a standardized space, that is first presented hereafter (Section 3.1) .The elasto-plastic formulation of the macroelement at the FE beam section follows in Section 3.2.

Standardized space
A standardized (dimensionless) space is adopted such that under uniaxial loading conditions the stress-resultant forces are equal to 1 for the positive loading direction, and −1 for the negative loading direction [32].In the proposed standardized space, the macroelement's failure surface is symmetrical with respect to the origin while in the real (dimensional) space, the failure surface presents an asymmetry with respect to the center of the axes for uniaxial loading conditions, see Section 3.2.The standardization procedure requires therefore a reference parameter and a shift parameter for the axial component, while for the shear and bending components only shift parameters are needed.
The standardized axial force F , shear force F and bending moment M are defined as: with  0 the reference parameter for the axial component and  *  ,  *  ,  *  the shift parameters for the axial, shear and bending components respectively.
The reference parameter  0 is defined as the midpoint between the maximum axial force   ,  (tensile force of positive sign) and the minimum axial force   ,  (compressive force of negative sign): The shift parameter  *  is defined as: The following four (4) parameters remaining to identify are:

𝜃
For this, the following characteristic states are introduced [33]:

Characteristic states
• 1st characteristic state: corresponds to the concrete elastic domain.
• 2nd characteristic state: corresponds to the peak values of the generalized forces-generalized displacement curves.The driving idea is that steel response controls the behavior of the RC section when tension is predominant, while concrete crushing controls compression.For the combined bending to tension area, failure is chosen to correspond to a maximum longitudinal steel strain equal to 7.5%, as in the Eurocode [39].For the combined bending to compression region, failure is considered reached at the maximum (absolute value) axial force.
Remark.The use of a criterion in compression in terms of concrete crushing and not buckling is chosen in order to limit the buckling risk.Furthermore, the use of a criterion in terms of concrete maximum compression stress proved to be highly conservative.Instead, a failure criterion based on the maximum (absolute value) axial force is reached only after several material points attain their ultimate resistance.

Parameter identification
The analytical computation of the maximum axial tensile force   ,  of the RC section corresponding to the 2nd characteristic state is not straightforward.Actually, the 3D numerical simulations [33] showed that the maximum axial tensile force in the RC section corresponds to a point beyond the concrete elastic limit, where concrete undergoes softening and steel hardening.  ,  can be analytically computed given the elastic modulus, the elastic limit, the ultimate strength and strain both for concrete and steel.Finally the calculation of  *  and  *  is not obvious as compound bending is prevailing for RC beam and column sections.

Maximum axial tensile force 𝐹 𝑡 𝑥, 𝑚𝑎𝑥
: it is calculated considering the strength and number of the steel rebars: where    is the steel strength and   the total area of the steel rebars.Numerical calculations [32] showed that this assumption can underestimate the maximum axial tensile force.Nevertheless, tensile forces, especially of large magnitude, are rarely expected in typical RC beams and columns.It can therefore be deduced that the adopted simplified approach does not significantly affect the performance of the macroelement.

Maximum axial compressive force 𝐹 𝑐 𝑥, 𝑚𝑎𝑥
: it is found according to the following formula where   is the concrete compressive strength,   the area of the concrete section and   the total area of the steel rebars.   is the current steel stress defined as follows: with  , the steel yield strain, dependent on the steel type and  ,ℎ the concrete crushing strain, taken either from experimental data or from the Eurocode.

Shift parameter 𝑀 *
: it is found from the moment-axial force interaction diagrams of Eurocode [39], without the security coefficients, as this permits to better reproduce the 3D FE numerical simulations [32,33].4. Shift parameter  *  : it is found using the method introduced by Rahal [40], based on a simplification of the Modified Compression Field Theory [41] for pure shear. *  is identified as the shear strength   multiplied by the effective sectional area (    ) and is provided by the graph of Fig. 2, given the compressive concrete strength   and the dimensionless variables   and   calculated as follows: with   and   the strengths of longitudinal and transverse steel bars respectively.  is the total area of the symmetrical longitudinal steel,   the area of the shear reinforcement bars within distance  ( being the spacing of the stirrups) and   and   the effective shear width and depth respectively.

Formulation
The stress-resultant plasticity model of the section (macroelement) is expressed in terms of generalized forces -generalized strains and is built within the plasticity framework.The failure surface is taken from [32], the loading surfaces are deduced from the failure surface and their convexity is guaranteed.Details are given hereafter:

Elastic behavior
The elastic behavior of the section is defined as: where F is the generalized force vector, ε the generalized strain vector (Eq.( 6)), ε  the generalized elastic strain vector and K  the elastic stiffness matrix of the section.They are given as follows: with   ,   and   the axial force, the shear force and the in-plane bending moment of the section respectively and   ,   and   the axial, transversal and bending stiffness of the section.

Plastic behavior
The total generalized strain vector is split into an elastic and a plastic part: Using Eqs. ( 24) and ( 26) the rate form of the generalized force vector is deduced:

Failure surface
The failure surface   is a homogeneous polynomial expression of degree 6 and it is written in the standardized (dimensionless) space defined in Section 3.1: equivalent to: where , ,  are the exponents of every component of the monomials and  (,,) the coefficients of the monomials calculated in [32] via the resolution of a convex optimization problem.They are provided in Table 1:

Loading surfaces and hardening evolution laws
The loading surfaces have the same form with the failure surface: with   ,   ,   the different uncoupled hardening parameters for each loading direction (, , ).
The proposed hardening evolution laws are exponential functions which tend asymptotically to unity.This choice presents two advantages: it correctly represents the global response of a typical RC section up to the peak, and it guarantees the non-interpenetration of the failure surface by the loading surface.The proposed hardening evolution laws are given by: with  0 ,  0 ,  0 (see Section 3.2.7)parameters that define the initial elastic domain in the standardized space while the   ,   ,   (see Section 3.2.7)parameters that control the hardening rate.  ,   ,   are the internal hardening variables, non-negative functions of the cumulative plastic flow of each loading direction, defined as: The proposed formulation does not adopt a kinematic hardening law, as this could lead to interpenetration of the failure surface by the loading surface without the appropriate tangent rule.A tangent rule [42] has been used in the previous works of [43,44], where the failure surface was of 2nd degree and by a change of variables the surfaces reduce to circles.In this work however, given the complexity of the failure surface, the tangent rule is not obvious and therefore an isotropic formulation is adopted in the following.

Plastic potential function and normality condition
An associative flow rule is assumed and therefore the plastic generalized strain evolution is computed as: Combining Eqs. ( 31) and ( 33) the 1st order differential equations of the different hardening evolution laws are: where h ( F , r) is the vectorial expression of the hardening evolution laws.

Persistency condition
The persistency condition ensures that for nonzero ̄ε   , λ is positive and the actual stress point remains on the boundary of the loading surface, i.e.   = 0 [45]: Eq. ( 35) implies: Combination of Eqs. ( 36) and ( 27) permits to calculate λ as follows:

Parameter identification
In order to illustrate the influence of  0 ,  0 ,  0 and   ,   ,   , see Eq. ( 31), a parametrical study on a RC cantilever beam subjected to a shear (transverse) load has been presented in [33].In particular, the force-displacement curves are obtained for different values of  0 ,  0 ,  0 and   ,   ,   (see Fig. 3).Following the results of the 3D FE numerical simulations for different RC sections [33] and comparison with experimental data lead us [33] to select the values of  0 ,  0 ,  0 in the range of 0.35-0.8and   = 500,   = 250,   = 250.

Numerical integration
A classical backward (implicit) Euler scheme of the Return Mapping Algorithm is used for the numerical implementation of the macroelement [45,46].The generalized forces are calculated for a given history of generalized strains at the section level.During each load increment, a trial elastic prediction of the sectional response is considered and a plastic correction is carried out if the predicted stress point is situated beyond the loading surface, to bring it back to the current loading surface according to  +1 = 0, a non-linear equation that it is solved using the Newton-Raphson algorithm.

Convergence rate
In order to improve the numerical performance of the proposed macroelement, a numerical tangent stiffness matrix is adopted, inspired by the approach implemented in the finite element code LAGAMINE (University of Liège) [47][48][49], where the terms of the compliance matrix are numerically calculated by finite differences between a perturbed state and a non-perturbed state.More specifically, the sectional stiffness matrix is updated applying strain perturbations successively for each loading direction.The magnitude of each directional strain perturbation is fixed at 1. −8 .A new state is calculated for each perturbation by applying the resolution algorithm of Section 3.2.8.The updated (at each iteration) sectional stiffness matrix K  is given as: where  = 1, 2, 3 are the loading directions, F   the pertubated force vector, F   the converged force vector and  ε  the applied strain perturbation vector in the loading direction .
Numerical simulations of a cantilever beam subjected to a transversal monotonic displacement at its free edge are conducted hereafter to illustrate the performance of the adopted strategy.One Timoshenko beam FE is used for the spatial discretization and the load is divided in 1000 steps.Fig. 4 presents the number of necessary iterations per step using the elastic matrix and the numerical stiffness matrix obtained through strain perturbations.It can be observed that the numerical performance is significantly improved.39) and ( 40) and of alternate sign, using the elastic stiffness and modified stiffness, Eq. ( 42).

Cyclic loading conditions
The proposed model is built within the elasto-plasticity framework and therefore, the unloading response is considered elastic.This does not provide a realistic response for earthquake loading as the model cannot capture different phenomena appearing during cyclic loads, such as stiffness degradation and pinching.Two simplified strategies are proposed hereafter to improve the performance of the model: the former for loading cycles of constant sign and the latter for loading cycles of alternate sign; the driving idea being that the response of the reinforced concrete element beyond a certain level of concrete damage is mainly controlled by steel.
• Cyclic loads of constant sign A limit value (equal to 0.8) is assigned to each hardening variable    ,    ,    , such as if one or more of the following conditions are satisfied: (39) then, accordingly, the axial, shear, and bending components of the sectional stiffness matrix equal the steel stiffness for loading and unloading conditions, introducing thus a stiffness degradation during unloading: The global response of a cantilever beam subjected to a cyclic transversal load of constant sign is presented in Fig. 5 (left) for two cases; using the elastic stiffness and the modified stiffness (Eq.( 40)).It can be observed that when the modified sectional stiffness is adopted for important applied displacements, the unloading slope decreases.

• Cyclic loads of alternate sign
For the case of cyclic loads of alternate sign, the method to introduce stiffness degradation is based on the well-known steel constitutive model developed by Menegotto and Pinto [51].This model provides Eq. ( 41) to describe the evolution of the steel elastic modulus at the semi-cycle under tension: where  0 is the steel Young modulus,  the plastic strain and  1 and  2 two material constants, the former depending on the geometric properties of the rebars and the latter on the steel mechanical properties.The values assigned to these constants can be found in the technical report of [50].
Eq. ( 41) is also adopted for the proposed macroelement, for each component independently, as: where K 0 is the initial sectional stiffness matrix.The global response of a cantilever beam subjected to a cyclic transversal displacement of alternate sign is presented in Fig. 5 (right), for two cases: using the elastic stiffness and the modified stiffness according to Eq. ( 42).It can be observed that Eq. ( 42) provides a simple way to introduce unloading stiffness degradation and pinching effects for the case of cyclic loads of alternate sign.This simplified approach was also validated through comparison with experimental results in [33].

Synopsis of the macroelement parameters
A synopsis of the eighteen ( 18) macroelement parameters and the proposed identification procedure that can be used by practitioners is presented in Table 2.

Modeling up to failure
The purpose of this article being the simplified evaluation of the response of RC beams and columns up to failure, the axial failure mode is excluded as it is the least probable to happen.Furthermore, as already stated, the dominant failure mechanism depends on the geometrical and mechanical properties of the RC element under study [27].Considering that the design of typical RC members is carried out by preventing shear failure (using an adequate number and spacing of stirrups), it is assumed in the following that bending is the dominant failure mode.The displacement field is therefore hereafter enhanced only by a strong discontinuity at the rotational component.The standard interpolation functions of the FE Timoshenko beam of Section 2, Eqs. ( 4) and (5), are used for the other components of the displacement field.Upon activation of the cohesive part, the behavior is no longer described by the macroelement; the bending component is described by the relation between cohesive moment and rotational jump, while the other two components are defined by linear elastic uniaxial relations.

Enhanced formulation
The adopted enhanced formulation is the one introduced by [7].In particular, among the four formulations introduced by the authors [7], the 000 formulation is chosen as it corresponds to a constant axial, bending and shear finite element with one-point integration.The plastic hinge is therefore considered activated at the center of the Timoshenko beam FE element of Section 2 (one Gauss integration point in the middle of the beam).
The enhanced generalized rotational field of the beam section is given by: where   is the total rotation of the section decomposed into a continuous part Θ and a discontinuous part Θ , Ū the element nodal displacements and Ū  the elementary displacement jumps defined as follows: are the Timoshenko beam shape functions of the rotational field (see Eqs. ( 3), ( 4), ( 5)) and M the enhancement shape functions of the rotational field.
Given the singularity of the corresponding strain field M is decomposed into the two following terms: with     the Heaviside function and ̄M  () containing additional functions to ensure the continuity of the rotational field between the finite elements.These functions are not a priori known.Their identification procedure is detailed in Section 4.1.1,based on the kinematic assumptions describing the discontinuity for a zero hinge mode (fully opened plastic hinge, i.e. fully opened crack).

Identification procedure
The curvature field   (, ) of the Timoshenko beam's section is obtained by derivation of the rotational field of Eq. ( 43): where B () contains the derivatives of the shape functions, with respect to , for the curvature component of the Timoshenko beam (see Eqs. ( 6), ( 9)) and Ḡ  () is the matrix of the derivatives of the enhancement interpolation functions, with respect to , of the following form: with     the Dirac function.Using the above decomposition of Ḡ  into continuous and discontinuous parts, the curvature field corresponding to Eq. ( 46) is rewritten as: The ̄Ḡ   functions are found considering the kinematics of the discontinuity at the final state.According to [7], at the zero hinge mode (fully opened plastic hinge), the end nodal rotations of the Timoshenko beam finite element of length   can be written as (see Fig. 6): with   the geometrical position of the discontinuity within the finite element (in our case   =   ∕2).Based on Eq. ( 49), the rotational jump component is thus given as: At the final state, the complete opening of the discontinuity results in a rigid body mode and the regular part of the enhanced curvature (see Eq. ( 48)) vanishes: Considering that , at the final state, the regular part of the curvature is expressed as: By introducing the rotational jump definition of Eq. ( 50) in Eq. ( 52), ̄Ḡ   is deduced:

Global cohesive model
The constitutive behavior of the discontinuity in the rotational field is described by means of a global cohesive model.This model reproduces a softening behavior that associates the moment at the plastic hinge to the rotation jump, and it is activated when the moment of the section reaches an ultimate value.

Activation and behavior
The curvature capacity threshold of the RC element, that will be used in Section 4.3 in the activation criterion (Eq.( 58)), is defined hereafter as : (54) where the critical value assigned to the curvature is given as: where  is the length of the RC element and the denominator is multiplied by 100 in order to transform    given as a percentage in Eq. (56).The identification procedure of    is given in Section 4.2.2.Further increase of the applied loading results in decrease of the cohesive moment  ℎ  and in increase of the rotational jump Θ .In particular, it is assumed that when the cohesive model is activated, the RC section has reached a high level of non-linearity and its behavior is controlled by the behavior of the reinforcing steel bars only (and not of concrete).Stiffness degradation and pinching effect are thus not taken into account by the cohesive model but by the modification of the unloading modulus of the continuous model (Section 3.2.10).As the RC sections are symmetric, the cohesive model is supposed identical for positive and negative bending.The behavior for a cyclic loading of constant sign is presented in the left part of Fig. 7, while for an alternate sign in the right.A vertical unloading/reloading slope is adopted to take into account permanent jumps (Fig. 7).

Parameter identification
The ultimate rotational capacity    (in percent) of Eq. ( 55) is calculated from the empirical relation provided in [52], based on regression analysis on experimental on RC rectangular sections: with: • ∕ the shear span ratio ( is the length of the RC element and  is the depth of the RC section).
•  the normalized steel ratio, defined as (where   is the area of the longitudinal reinforcement in tension and   is the cross-sectional area).
•   the steel yield limit (in ).
•   the normalized axial force defined as  ∕  , with  the axial force (positive for compression) and   =   (where  is the width and  the depth of the RC element).
•   the confinement ratio (in percent), defined as (where   is the area of the transverse steel parallel to the direction of loading and  ℎ the spacing of the stirrups).
Finally, the softening modulus  of the cohesive moment-rotational jump is chosen equal to −7% of the initial bending modulus in Section 5, as proposed by [53].It should be mentioned however that a value equal to −0.7% of the initial bending modulus seems more suitable for certain cyclic loading tests, see for example [33].

Numerical integration
The 1D global cohesive model of Fig. 7 is numerically integrated using the classical computational plasticity framework [45,54].As stated in [7], the enhanced by strong discontinuities approach is a tool for modeling localized dissipative mechanisms and allows for capturing energy dissipation in large-scale problems.The displacement jump plays therefore a similar role to an internal variable and a backward (implicit) Euler scheme of the Return Mapping Algorithm families can be used for the numerical integration.The numerical implementation and the computational procedure (including static condensation at an elementary level) adopted are inspired by [28,29].The integration algorithm of the cohesive model is briefly presented hereafter (more detailed information can be found in [29]).
The adopted cohesive model is a linear type plastic softening law (see Fig. 7) that describes the relationship between the cohesive moment  ℎ  and the rotational jump Θ : where  is the softening modulus (of negative value).The activation criterion is given by: where    corresponds to    and it is numerically calculated by the FE code.It is assumed that when the rotational discontinuity is activated, it remains activated for the rest of the loading history and the failure criterion is the same as the activation criterion   .The evolution of the elementary variable of the rotational jump is given by: where ̇λ is the plastic multiplier defining the magnitude of the rotational jump rate.The persistency condition expresses that once the discontinuity is activated the plastic multiplier ̇λ is positive and thus the failure criterion is satisfied.This implies that: and therefore the expression for the calculation of the plastic multiplier is known.The cohesive model internal variable Θ is computed at every time step (once the activation criterion is satisfied) for every Gauss point at every iteration .The elastic predictor of the current step  + 1 is equal to the previous step  value: and the trial cohesive moment is given by the following expression satisfying the local equilibrium at the discontinuity level within the element: where  is the number of integration points of the beam element (in our case  = 1) and (  ) is the integration weight of the Gauss point with coordinate   (in our case the Gauss point is in the middle of the FE).

The expression of 𝑀 𝑡𝑟𝑖𝑎𝑙 𝜃 𝑛+1
(  ) is: The trial failure criterion is evaluated: If this condition (64) holds the discontinuity is not activated and thus the variables of the cohesive model are updated as: If discontinuity is activated the rotational jump at the current step is given by: The updated cohesive force at the step  + 1 can be expressed as: The above expression allows for the calculation of the plastic multiplier: and the cohesive tangent modulus is defined as:

Numerical applications
In order to validate the performance of the new generalized Timoshenko beam with embedded rotation discontinuity coupled with the 3D macroelement, two case studies are presented in this section.A two story RC frame subjected to a monotonic and a cyclic loading (Section 5.1), and a RC beam subjected to an earthquake (Section 5.2).In both cases, comparisons are made with the experimental results.
A mesh convergence study can be found in [33] (not reproduced hereafter), where it can be seen that the proposed finite element model provides the same results for a reinforced concrete cantilever column discretized with 4-8-16 and 32 beam elements.

A two story RC frame [55]
A two story RC frame was experimentally tested by [55].Its geometry and reinforcement details are depicted in Fig. 8.All beam and columns sections have the same geometrical and material properties, the only difference being the concrete cover (30 mm for the beams and 20 mm for the columns).
This structure has often been used for validation purposes, see for example [23,26,27,29].In these works, the researchers assigned different material parameters to the beams and the columns, justified by the fact that the columns are subjected to significant axial compressive loads resulting to higher bending resistances.This however is not necessary for the proposed new beam as the coupling between the different stress-resultant components is explicitly taken into account via the 3D failure surface.
Thirty (30) finite elements of 0.5 m length are used for the spatial discretization (seven (7) for each beam and eight (8) for each column).The columns are considered fixed at their base.Fig. 9 presents the FE mesh and the boundary conditions.Numerical simulations are conducted for monotonic and cyclic (constant sign) lateral loads, the latter to reproduce the experimental conditions.More specifically, equal constant axial loads are applied on top of the columns, while lateral displacements are incrementally imposed at the top of the north column (displacement increments of 0.06 mm) (see Figs. 8 and 9).The model parameters are given in Table 3.
The numerical response is compared to the experimental one for monotonic and cyclic applied lateral loads in Figs. 10 and 11 respectively.Numerical simulations were done with Matlab [56].Fig. 10 shows that the macroelement successfully captures the peak of the response.Fig. 10 (right) provides the individual and total (sum) contributions of elements 8 and 16 (up to 0.3 m, as no experimental data are available after this point).It illustrates that the cohesive model is activated for a displacement equal to 0.13 m for the element number 8, while for 0.19 m for the element number 16.Even though softening is not observed at the global response in Fig. 10 (left), there are elements that exhibit softening.Fig. 11 shows that for the case of a cyclic load without sign change, the model is also able to accurately reproduce loading and unloading.The global response of the RC structure for severe loading is mainly controlled by the reinforcing steel.The modification of the stiffness components such that they become equal to the reinforcement values near the peak, see Section 3.2.10,results therefore to a realistic prediction of the cyclic response.

A simply supported RC beam [57]
[57] experimentally studied an intermediate level beam, part of a two story RC frame, subjected to a vertical earthquake loading.During the experiment, pivot connections of the beam with the adjacent columns were assumed, and thus the studied beam case is similar to a three-point bending test.The geometry and reinforcing details of the beam are shown in Fig. 12 (left).The concrete cover is 0.01 m.
In the following, the seismic loading is numerically reproduced under quasi-static conditions by applying the experimental vertical displacement time history provided in Fig. 12 (right).The vertical displacements increments are of 0.06 mm.Four (4)   Timoshenko beam finite elements of 0.425 m length are used for the spatial discretization.In Fig. 13, the FE mesh and the boundary conditions are presented and the model parameters are presented in Table 4.
Fig. 14 (left) compares the numerical response with the experimental data.Numerical simulations were again made using Matlab [56].The cohesive model was not found activated in this example.Although the elastic phase is well captured, the model at certain moments over/underestimates the response and a calibration procedure is necessary to achieve a more satisfactory response.In order to improve the response, the stiffness degradation parameter 1 is calibrated and taken equal to 0.1.Fig. 14 (right) shows that by improving the unloading/reloading response, a better fit of the numerical response to the experimental curve is achieved at late stages of loading.

Conclusions
A generalized finite element beam with an embedded rotation discontinuity coupled with a 3D macroelement is proposed to assess, till complete failure the vulnerability of symmetrically RC frame structures subjected to quasi-static and dynamic loading.Parameter identification procedures are discussed, based on simplified methods to ensure every day practice.The numerical performance of this novel simulation tool is checked considering two case studies, a two story RC frame subjected to a monotonic and a cyclic loading, and a RC simply supported beam subjected to an earthquake.The comparison with the experimental results shows that the generalized finite element beam with an embedded rotation discontinuity coupled with the 3D macroelement captures accurately the peak of the response, the initiation (or not) of the softening phase and it is able to reproduce the behavior till complete failure.The new element can be used in design offices for fast calculations that take into account multidimensional complicated actions in a simplified way.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Influence of  0 ,  0 ,  0 and   ,   ,   on the global response of a RC cantilever beam subjected to a transverse load.

Fig. 4 .
Fig. 4. Number of necessary iterations per step for the convergence of the global Newton algorithm using the elastic matrix and the numerical stiffness matrix with strain perturbations.

Fig. 5 .
Fig. 5.Comparison of the global responses of a cantilever beam subjected to a cyclic transversal load of constant sign using the elastic stiffness and the modified stiffness, Eqs.(39) and (40) and of alternate sign, using the elastic stiffness and modified stiffness, Eq. (42).

Fig. 6 .
Fig. 6.Sketch of the plastic hinge mode: initiation of the plastic hinge (up) and fully opened hinge (bottom).

Table 1
[32]ficients of monomials for the failure envelope expression of degree 6[32]for a symmetrically RC section.

Table 2
Macroelement parameters and proposed identification procedures.

Table 3
Model parameters for the two story RC frame.Two story RC frame: Numerical vs. experimental global response (cyclic loading with constant sign).

Table 4
Model parameters for the RC simply supported beam.