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

Androniki-Anna Doulgeroglou, Panagiotis Kotronis, Giulio Sciarra, Catherine Bouillon

► To cite this version:

Androniki-Anna Doulgeroglou, Panagiotis Kotronis, Giulio Sciarra, Catherine Bouillon. A generalized Timoshenko beam with embedded rotation discontinuity coupled with a 3D macroelement to assess the vulnerability of reinforced concrete frame structures. Finite Elements in Analysis and Design, 2024, 242, pp.104234. 10.1016/j.finel.2024.104234. hal-04671938

HAL Id: hal-04671938

https://cnrs.hal.science/hal-04671938v1

Submitted on 16 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Androniki-Anna Doulgeroglou<sup>a,c</sup>, Panagiotis Kotronis<sup>b,*</sup>, Giulio Sciarra<sup>b</sup>, Catherine Bouillon<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Macroelement

Reinforced concrete

Vulnerability

Timoshenko beam

Strong discontinuity

## ABSTRACT

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.

## 1. 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–8] for generalized finite element (FE) beams.

<footnote><sup>a</sup> Arts et Métiers Institute of Technology, Université de Bordeaux, CNRS, Institut de Mécanique et d'Ingénierie (I2M), UMR 5295, Talence, F-33405, France</footnote>

<footnote><sup>b</sup> Nantes Université, École Centrale Nantes, CNRS, Institut de Recherche en Génie Civil et Mécanique (GeM), UMR 6183, Nantes, F-44000, France</footnote>

<footnote><sup>c</sup> Groupe Essor, 2 rue Pierre Gilles de Gennes, Lons, F-64140, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: androniki-anna.doulgeroglou@ensam.eu (A.-A. Doulgeroglou), Panagiotis.Kotronis@ec-nantes.fr (P. Kotronis), Giulio.Sciarra@ec-nantes.fr (G. Sciarra).</footnote>
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–11] for shallow foundations, [12–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–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–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 (stress resultant 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 A defines the vector A, A the matrix A, A the variable A in a standardized (dimensionless) space, A the derivative of A with respect to time, A' the derivative of A with respect to the coordinate x and T the symbol of the transpose. In Section 4, A is a variable assigned to a continuous field while A' a variable assigned to a discontinuous field.

## 2. 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 $u_p = (u_x, u_y)$ of a point P $\equiv (x, y)$ of the beam is expressed in terms of the displacement of the gravity center G of the section $U_G^T = (U_x, U_y)$ and the rotation $\Theta_z$ of the section as follows:

$$u_x(x, y) = U_x(x) - y\Theta_z(x), \quad u_y(x, y) = U_y(x) \qquad (1)$$

Using the small strain assumption, the strain components are given by:

$$\frac{\partial u_x}{\partial x} = U'_x(x) - y\Theta'_z(x), \quad \frac{\partial u_y}{\partial x} + \frac{\partial u_x}{\partial y} = U'_y(x) - \Theta_z(x) = \beta_y \neq 0 \qquad (2)$$

Consider now a 2D FE beam of length $L_e$, with two nodes i and j and three degrees of freedom per node: axial displacement ($U_{xi}$ and $U_{xj}$), transversal displacement ($U_{yi}$ and $U_{yj}$) and rotation ($\Theta_{zi}$ and $\Theta_{zj}$), see Fig. 1 (right). The generalized displacement vector $U_s$ is function of the nodal displacement vector $U$ as follows [34,35]:

$$U_s = \underline{N} U \qquad (3)$$

$$\begin{bmatrix} U_x \\ U_y \\ \Theta_z \end{bmatrix} = \begin{bmatrix} N_1^a & 0 & 0 & N_4^a & 0 & 0 \\ 0 & N_2^t & 0 & 0 & N_5^t & 0 \\ 0 & 0 & N_3^\theta & 0 & 0 & N_6^\theta \end{bmatrix} \begin{bmatrix} U_{xi} \\ U_{yi} \\ \Theta_{zi} \\ U_{xj} \\ U_{yj} \\ \Theta_{zj} \end{bmatrix}, \qquad (4)$$

$$N_1^a = N_2^t = N_3^\theta = 1 - \frac{x}{L_e}, \quad N_4^a = N_5^t = N_6^\theta = \frac{x}{L_e} \qquad (5)$$

with $\underline{N}$ the matrix of the linear shape functions, $x \in [0, L_e]$, while the upper indexes a, t, $\theta$ denote the axial, transversal and rotational components respectively. By derivation of the displacement field and according to Eq. (2) the generalized strain vector $\epsilon_s$ is obtained as follows:

$$\epsilon_s = \begin{bmatrix} \epsilon_x \\ \beta_y \\ \kappa_\theta \end{bmatrix} = \begin{bmatrix} U'_x \\ \beta_y \\ \Theta'_z \end{bmatrix} = \begin{bmatrix} U'_x \\ U'_y - \Theta_z \\ \Theta'_z \end{bmatrix} = \underline{B}_{exact} U \qquad (6)$$

$$\begin{bmatrix} \epsilon_x \\ \beta_y \\ \kappa_\theta \end{bmatrix} = \begin{bmatrix} -\frac{1}{L_e} & 0 & 0 & \frac{1}{L_e} & 0 & 0 \\ 0 & -\frac{1}{L_e} & -1 + \frac{x}{L_e} & 0 & \frac{1}{L_e} & -\frac{x}{L_e} \\ 0 & 0 & -\frac{1}{L_e} & 0 & 0 & \frac{1}{L_e} \end{bmatrix} \begin{bmatrix} U_{xi} \\ U_{yi} \\ \Theta_{zi} \\ U_{xj} \\ U_{yj} \\ \Theta_{zj} \end{bmatrix} \qquad (7)$$

with $\epsilon_x$, $\beta_y$ and $\kappa_\theta$ the axial strain, the shear strain and the curvature at the beam axis that passes through the point G respectively. Considering that $\beta_y$ can be written as

$$\beta_y = \left( \frac{U_{yj} - U_{yi}}{L_e} - \frac{\Theta_{zj} + \Theta_{zi}}{2} \right) + \frac{1}{2}(\Theta_{zi} - \Theta_{zj})\zeta \quad \text{with } \zeta = \frac{2}{L_e}x - 1 \qquad (8)$$
<figure><img src="image_4.png" /><figcaption>Fig. 1. 2D Timoshenko beam (left), 2-node FE Timoshenko beam (right).</figcaption></figure>

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 $\underline{B}$ below (see also [34,35]):

$$ \begin{bmatrix} \epsilon_x \\ \beta_y \\ \kappa_\theta \end{bmatrix} = \begin{bmatrix} -\frac{1}{L_e} & 0 & 0 & \frac{1}{L_e} & 0 & 0 \\ 0 & -\frac{1}{L_e} & -\frac{1}{2} & 0 & \frac{1}{L_e} & -\frac{1}{2} \\ 0 & 0 & -\frac{1}{L_e} & 0 & 0 & \frac{1}{L_e} \end{bmatrix} \begin{bmatrix} U_{xi} \\ U_{yi} \\ \Theta_{zi} \\ U_{xj} \\ U_{yj} \\ \Theta_{zj} \end{bmatrix} = \underline{B} \underline{U} \qquad (9) $$

Using the following form of the principle of virtual work [34,35]

$$ \int_{0}^{L_e} \left( F_x \frac{d}{dx} \delta U_x + F_y \frac{d}{dx} \delta \beta_y + M_{\theta} \frac{d}{dx} \delta \Theta_z \right) dx - w_{ext} = 0 \qquad (10) $$

together with the constitutive description of the section in terms of generalized forces (see Section 3.2), the following expression is obtained:

$$ \int_{0}^{L_e} \delta [U'_{x} \quad \beta_{y} \quad \Theta'_{z}] \underline{K}_{s} [U'_{x} \quad \beta_{y} \quad \Theta'_{z}]^{T} dx - w_{ext} = 0 \qquad (11) $$

Combining expressions (6) and (11), the element stiffness matrix is given by:

$$ \underline{K}_e = \int_{0}^{L_e} \underline{B}^T \underline{K}_s \underline{B} dx \qquad (12) $$

with $\underline{K}_s$ the stiffness matrix of the section (see Section 3.2).

As usual for FE beams, the element mass matrix is

$$ \underline{M}_e = \int_{0}^{L_e} \underline{N}^T \underline{M}_s \underline{N} dx \qquad (13) $$

with $\underline{M}_s$ the mass matrix of the section given as:

$$ \underline{M}_s = \begin{bmatrix} \rho S & 0 & 0 \\ 0 & \rho S & 0 \\ 0 & 0 & \rho I \end{bmatrix} \qquad (14) $$

where $\rho$ the material density, $S$ the area and $I$ 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

$$ \underline{F}_e = \int_{0}^{L_e} \underline{B}^T \underline{F}_s dx \qquad (15) $$

with $\underline{F}_s$ 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].
### 3. 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.

#### 3.1. 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 $\hat{F}_x$, shear force $\hat{F}_y$ and bending moment $\hat{M}_\theta$ are defined as:

$$ \hat{F}_x = \frac{F_x - F_{x0}}{F_x^*}, \quad \hat{F}_y = \frac{F_y}{F_y^*}, \quad \hat{M}_\theta = \frac{M_\theta}{M_\theta^*} \qquad (16) $$

with $F_{x0}$ the reference parameter for the axial component and $F_x^*$, $F_y^*$, $M_\theta^*$ the shift parameters for the axial, shear and bending components respectively.

The reference parameter $F_{x0}$ is defined as the midpoint between the maximum axial force $F_{x,max}^t$ (tensile force of positive sign) and the minimum axial force $F_{x,max}^c$ (compressive force of negative sign):

$$ F_{x0} = (F_{x,max}^t + F_{x,max}^c)/2 \qquad (17) $$

The shift parameter $F_x^*$ is defined as:

$$ F_x^* = (F_{x,max}^t - F_{x,max}^c)/2 \qquad (18) $$

The following four (4) parameters remaining to identify are:

$$ F_{x,max}^t, \quad F_{x,max}^c, \quad F_y^*, \quad M_\theta^* $$

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

#### 3.1.1. 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.

#### 3.1.2. Parameter identification

The analytical computation of the maximum axial tensile force $F_{x,max}^t$ 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. $F_{x,max}^c$ 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 $F_y^*$ and $M_\theta^*$ is not obvious as compound bending is prevailing for RC beam and column sections. The following simplified procedure is therefore proposed to calibrate $F_{x,max}^t$, $F_{x,max}^c$, $F_y^*$ and $M_\theta^*$:
<figure><img src="image_5.png" /><figcaption>Fig. 2. Normalized strength curves, adapted from [40], used to calculate $F_s^*$.</figcaption></figure>

1. Maximum axial tensile force $F_{x,max}^t$: it is calculated considering the strength and number of the steel rebars:

$$F_{x,max}^t = f_s^{ult} A_s \qquad (19)$$

where $f_s^{ult}$ is the steel strength and $A_s$ 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.

2. Maximum axial compressive force $F_{x,max}^c$: it is found according to the following formula

$$F_{x,max}^c = f_c A_c + f_s^{cur} A_s \qquad (20)$$

where $f_c$ is the concrete compressive strength, $A_c$ the area of the concrete section and $A_s$ the total area of the steel rebars. $f_s^{cur}$ is the current steel stress defined as follows:

$$f_s^{cur} = E_s \epsilon_{c,crush} \quad \text{if} \quad \epsilon_{c,crush} < \epsilon_{steel,yield} \qquad (21)$$

$$f_s^{cur} = f_s^{ult} \quad \text{if} \quad \epsilon_{c,crush} \ge \epsilon_{steel,yield} \qquad (22)$$

with $\epsilon_{steel,yield}$ the steel yield strain, dependent on the steel type and $\epsilon_{c,crush}$ the concrete crushing strain, taken either from experimental data or from the Eurocode.

3. Shift parameter $M_0^*$: 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 $F_s^*$: it is found using the method introduced by Rahal [40], based on a simplification of the Modified Compression Field Theory [41] for pure shear. $F_s^*$ is identified as the shear strength $v_u$ multiplied by the effective sectional area $(b_u d_v)$ and is provided by the graph of Fig. 2, given the compressive concrete strength $f_c$ and the dimensionless variables $\omega_L$ and $\omega_t$ calculated as follows:

$$\omega_L = \frac{A_L f_{yL}}{b_u d_v f_c}, \quad \omega_t = \frac{A_v f_{yt}}{b_u s f_c} \qquad (23)$$

with $f_{yL}$ and $f_{yt}$ the strengths of longitudinal and transverse steel bars respectively. $A_L$ is the total area of the symmetrical longitudinal steel, $A_v$ the area of the shear reinforcement bars within distance $s$ ($s$ being the spacing of the stirrups) and $b_u$ and $d_v$ the effective shear width and depth respectively.

## 3.2. 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:

### 3.2.1. Elastic behavior

The elastic behavior of the section is defined as:

$$F_s = \underline{K}_s^{el} \epsilon_s = \underline{K}_s^{el} \epsilon_s^{el} \qquad (24)$$
<table><caption>**Table 1**
Coefficients of monomials for the failure envelope expression of degree 6 [32] for a symmetrically RC section.</caption><thead><tr><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>a<sub>4</sub></th><th>a<sub>5</sub></th><th>a<sub>6</sub></th><th>a<sub>7</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>14.03</td><td>0.03</td><td>12.26</td><td>0.02</td><td>1</td></tr><tr><td>a<sub>8</sub></td><td>a<sub>9</sub></td><td>a<sub>10</sub></td><td>a<sub>11</sub></td><td>a<sub>12</sub></td><td>a<sub>13</sub></td><td>a<sub>14</sub></td></tr><tr><td>0.01</td><td>-12.73</td><td>0</td><td>-17.97</td><td>-0.06</td><td>-3.34</td><td>8.29</td></tr><tr><td>a<sub>15</sub></td><td>a<sub>16</sub></td><td>a<sub>17</sub></td><td>a<sub>18</sub></td><td>a<sub>19</sub></td><td>a<sub>20</sub></td><td>a<sub>21</sub></td></tr><tr><td>-0.05</td><td>35.83</td><td>0.13</td><td>11.09</td><td>0</td><td>-22.46</td><td>-0.18</td></tr><tr><td>a<sub>22</sub></td><td>a<sub>23</sub></td><td>a<sub>24</sub></td><td>a<sub>25</sub></td><td>a<sub>26</sub></td><td>a<sub>27</sub></td><td>a<sub>28</sub></td></tr><tr><td>15.42</td><td>5.56</td><td>0.1</td><td>12.69</td><td>-0.02</td><td>-5.51</td><td>1</td></tr></tbody></table>

where $\underline{F}_s$ is the generalized force vector, $\underline{\xi}_s$ the generalized strain vector (Eq. (6)), $\underline{\xi}_s^{el}$ the generalized elastic strain vector and $\underline{K}_s^{el}$ the elastic stiffness matrix of the section. They are given as follows:

$$ \underline{F}_s = \begin{bmatrix} F_x \\ F_y \\ M_\theta \end{bmatrix}, \quad \underline{\xi}_s = \begin{bmatrix} \underline{\xi}_x \\ \underline{\beta}_y \\ \underline{\kappa}_\theta \end{bmatrix}, \quad \underline{K}_s^{el} = \begin{bmatrix} K_x & 0 \\ 0 & K_y \\ 0 & 0 \end{bmatrix} \qquad (25) $$

with $F_x$, $F_y$ and $M_\theta$ the axial force, the shear force and the in-plane bending moment of the section respectively and $K_x$, $K_y$ and $K_\theta$ the axial, transversal and bending stiffness of the section.

### 3.2.2. Plastic behavior

The total generalized strain vector is split into an elastic and a plastic part:

$$ \underline{\xi}_s = \underline{\xi}_s^{el} + \underline{\xi}_s^{pl} \qquad (26) $$

Using Eqs. (24) and (26) the rate form of the generalized force vector is deduced:

$$ \dot{\underline{F}}_s = \underline{K}_s^{el} (\dot{\underline{\xi}}_s - \dot{\underline{\xi}}_s^{pl}) \qquad (27) $$

### 3.2.3. Failure surface

The failure surface $f_U$ is a homogeneous polynomial expression of degree 6 and it is written in the standardized (dimensionless) space defined in Section 3.1:

$$ f_U = \sum_{\substack{i,j,k=0 \\ |i+j+k|=6}}^{6} a_{I(i,j,k)} \hat{F}_x^i \hat{F}_y^j \hat{M}_\theta^k - 1 \qquad (28) $$

equivalent to:

$$ \begin{align} f_U ={}& a_1 \hat{F}_x^6 + a_2 \hat{F}_x^5 \hat{M}_\theta + a_3 \hat{F}_x^4 \hat{M}_\theta^2 + a_4 \hat{F}_x^3 \hat{M}_\theta^3 + a_5 \hat{F}_x^2 \hat{M}_\theta^4 + a_6 \hat{F}_x \hat{M}_\theta^5 + a_7 \hat{M}_\theta^6 \nonumber \\ & + a_8 \hat{F}_x^5 \hat{F}_y + a_9 \hat{F}_x^4 \hat{M}_\theta \hat{F}_y + a_{10} \hat{F}_x^3 \hat{M}_\theta^2 \hat{F}_y + a_{11} \hat{F}_x^2 \hat{M}_\theta^3 \hat{F}_y + a_{12} \hat{F}_x \hat{M}_\theta^4 \hat{F}_y + a_{13} \hat{M}_\theta^5 \hat{F}_y \nonumber \\ & + a_{14} \hat{F}_x^4 \hat{F}_y^2 + a_{15} \hat{F}_x^3 \hat{M}_\theta \hat{F}_y^2 + a_{16} \hat{F}_x^2 \hat{M}_\theta^2 \hat{F}_y^2 + a_{17} \hat{F}_x \hat{M}_\theta^3 \hat{F}_y^2 + a_{18} \hat{M}_\theta^4 \hat{F}_y^2 + a_{19} \hat{F}_x^3 \hat{F}_y^3 \nonumber \\ & + a_{20} \hat{F}_x^2 \hat{M}_\theta \hat{F}_y^3 + a_{21} \hat{F}_x \hat{M}_\theta^2 \hat{F}_y^3 + a_{22} \hat{M}_\theta^3 \hat{F}_y^3 + a_{23} \hat{F}_x^2 \hat{F}_y^4 + a_{24} \hat{F}_x \hat{M}_\theta \hat{F}_y^4 + a_{25} \hat{M}_\theta^2 \hat{F}_y^4 \nonumber \\ & + a_{26} \hat{F}_x \hat{F}_y^5 + a_{27} \hat{M}_\theta \hat{F}_y^5 + a_{28} \hat{F}_y^6 - 1 \tag{29} \end{align} $$

where $i, j, k$ are the exponents of every component of the monomials and $a_{I(i,j,k)}$ the coefficients of the exponents calculated in [32] via the resolution of a convex optimization problem. They are provided in Table 1:

### 3.2.4. Loading surfaces and hardening evolution laws

The loading surfaces have the same form with the failure surface:

$$ f_L = \sum_{\substack{i,j,k=0 \\ |i+j+k|=6}}^{6} a_{I(i,j,k)} \left( \left( \frac{\hat{F}_x}{r_x} \right)^i \left( \frac{\hat{F}_y}{r_y} \right)^j \left( \frac{\hat{M}_\theta}{r_\theta} \right)^k - 1 \right) \qquad (30) $$

with $r_x, r_y, r_\theta$ the different uncoupled hardening parameters for each loading direction $(x, y, \theta)$.

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:

$$ \underline{r} = \begin{bmatrix} r_x \\ r_y \\ r_\theta \end{bmatrix} = \begin{bmatrix} 1 + (r_{x0} - 1)e^{-a_x p_x} \\ 1 + (r_{y0} - 1)e^{-a_y p_y} \\ 1 + (r_{\theta0} - 1)e^{-a_\theta p_\theta} \end{bmatrix} \qquad (31) $$
with $r_{x0}$, $r_{y0}$, $r_{\theta0}$ (see Section 3.2.7) parameters that define the initial elastic domain in the standardized space while the $a_x$, $a_y$, $a_\theta$ (see Section 3.2.7) parameters that control the hardening rate. $p_x$, $p_y$, $p_\theta$ are the internal hardening variables, non-negative functions of the cumulative plastic flow of each loading direction, defined as:

$$ \underline{p} = \begin{bmatrix} p_x \\ p_y \\ p_\theta \end{bmatrix}, \quad \dot{\underline{p}} = \begin{bmatrix} \dot{p}_x \\ \dot{p}_y \\ \dot{p}_\theta \end{bmatrix} = \begin{bmatrix} |\dot{\epsilon}_x^{pl}| \\ |\dot{\rho}_y^{pl}| \\ |\dot{k}_\theta^{pl}| \end{bmatrix} \qquad (32) $$

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.

### 3.2.5. Plastic potential function and normality condition

An associative flow rule is assumed and therefore the plastic generalized strain evolution is computed as:

$$ \underline{\epsilon}_s^{pl} = \dot{\lambda} \frac{\partial f_L}{\partial F_s} \Leftrightarrow \begin{bmatrix} \dot{\epsilon}_x^{pl} \\ \dot{\rho}_y^{pl} \\ \dot{k}_\theta^{pl} \end{bmatrix} = \dot{\lambda} \begin{bmatrix} \frac{\partial f_L}{\partial F_s} \\ \frac{\partial f_L}{\partial F_y} \\ \frac{\partial f_L}{\partial M_\theta} \end{bmatrix} \qquad (33) $$

Combining Eqs. (31) and (33) the 1st order differential equations of the different hardening evolution laws are:

$$ \dot{\underline{r}} = \begin{bmatrix} \dot{r}_x \\ \dot{r}_y \\ \dot{r}_\theta \end{bmatrix} = \begin{bmatrix} a_x(1-r_x)\dot{p}_x \\ a_y(1-r_y)\dot{p}_y \\ a_\theta(1-r_\theta)\dot{p}_\theta \end{bmatrix} = \dot{\lambda} \underbrace{\begin{bmatrix} a_x(1-r_x) \frac{\partial f_L}{\partial F_s} \\ a_y(1-r_y) \frac{\partial f_L}{\partial F_y} \\ a_\theta(1-r_\theta) \frac{\partial f_L}{\partial M_\theta} \end{bmatrix}}_{h_L(F_s, \underline{r})} \qquad (34) $$

where $h_L(F_s, \underline{r})$ is the vectorial expression of the hardening evolution laws.

### 3.2.6. Persistency condition

The persistency condition ensures that for nonzero $\underline{\epsilon}_s^{pl}$, $\dot{\lambda}$ is positive and the actual stress point remains on the boundary of the loading surface, i.e. $f_L = 0$ [45]:

$$ \dot{\lambda} f_L = 0 \quad (if \quad f_L = 0) \qquad (35) $$

Eq. (35) implies:

$$ \dot{f}_L = 0 \Leftrightarrow \frac{\partial f_L}{\partial F_s} \dot{F}_s + \frac{\partial f_L}{\partial r} \dot{r} = 0 \qquad (36) $$

Combination of Eqs. (36) and (27) permits to calculate $\dot{\lambda}$ as follows:

$$ \begin{align} \frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\epsilon}_s (\dot{\epsilon}_s - \dot{\epsilon}_s^{pl}) + \frac{\partial f_L}{\partial r} \dot{r} &= 0 \Leftrightarrow \nonumber \\ \frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\epsilon}_s \dot{\epsilon}_s - \frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\epsilon}_s^{pl} + \frac{\partial f_L}{\partial r} \dot{r} &= 0 \Leftrightarrow \nonumber \\ \frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\epsilon}_s \dot{\epsilon}_s - \frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\lambda} \frac{\partial f_L}{\partial F_s} + \frac{\partial f_L}{\partial r} \dot{r} &= 0 \Leftrightarrow \nonumber \\ \dot{\lambda} &= \frac{\frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\epsilon}_s}{\frac{\partial f_L}{\partial F_s} \underset{\underline{F}_s}{\overset{K}{\equiv}} \dot{\lambda} \frac{\partial f_L}{\partial F_s} - \frac{\partial f_L}{\partial r} h_L} \tag{37} \end{align} $$

### 3.2.7. Parameter identification

In order to illustrate the influence of $r_{x0}$, $r_{y0}$, $r_{\theta0}$ and $a_x$, $a_y$, $a_\theta$, 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 $r_{x0}$, $r_{y0}$, $r_{\theta0}$ and $a_x$, $a_y$, $a_\theta$ (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 $r_{x0}$, $r_{y0}$, $r_{\theta0}$ in the range of 0.35–0.8 and $a_x = 500$, $a_y = 250$, $a_\theta = 250$.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Influence of $r_{x0}$, $r_{y0}$, $r_{\theta0}$ and $a_s$, $a_y$, $a_\theta$ on the global response of a RC cantilever beam subjected to a transverse load.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

### 3.2.8. 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 $f_{L,n+1} = 0$, a non-linear equation that it is solved using the Newton-Raphson algorithm.

### 3.2.9. 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–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.e^{-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_s^{per}$ is given as:

$$K_s^{per}(i, i) = \frac{F_s^{per}(i) - F_s^{conv}(i)}{d\xi_s^{per}(i)} \qquad (38)$$

where $i = 1, 2, 3$ are the loading directions, $F_s^{per}$ the perturbed force vector, $F_s^{conv}$ the converged force vector and $d\xi_s^{per}$ the applied strain perturbation vector in the loading direction $i$.

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.
<figure><img src="image_8.png" /><figcaption>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).</figcaption></figure>

<table><caption>**Table 2**
Macroelement parameters and proposed identification procedures.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Identification</th></tr></thead><tbody><tr><td colspan="3">Linear elastic behavior</td></tr><tr><td>Axial stiffness</td><td>K<sub>x</sub></td><td>Geometrical &amp; mechanical properties</td></tr><tr><td>Shear stiffness</td><td>K<sub>y</sub></td><td></td></tr><tr><td>Bending stiffness</td><td>K<sub>&theta;</sub></td><td></td></tr><tr><td>Axial elastic limit</td><td>r<sub>x0</sub></td><td>Selected in the range</td></tr><tr><td>Shear elastic limit</td><td>r<sub>y0</sub></td><td>0.35 to 0.8</td></tr><tr><td>Bending elastic limit</td><td>r<sub>&theta;0</sub></td><td></td></tr><tr><td colspan="3">Standardization</td></tr><tr><td>Maximum axial tensile force</td><td>F<sup>l</sup><sub>x,max</sub></td><td>From material strength formulas</td></tr><tr><td>Maximum axial compressive force</td><td>F<sup>c</sup><sub>x,max</sub></td><td></td></tr><tr><td>Shear shift parameter</td><td>F<sup>s</sup><sub>y</sub></td><td>MCFT [40]</td></tr><tr><td>Bending shift parameter</td><td>M<sup>s</sup><sub>&theta;</sub></td><td>Moment-axial force interaction diagram [39]</td></tr><tr><td colspan="3">Hardening evolution laws</td></tr><tr><td>Axial rate of hardening</td><td>a<sub>x</sub></td><td>a<sub>x</sub> = 500, a<sub>y</sub> = 250, a<sub>&theta;</sub> = 250</td></tr><tr><td>Shear rate of hardening</td><td>a<sub>y</sub></td><td>or calibration by numerical simulations</td></tr><tr><td>Bending rate of hardening</td><td>a<sub>&theta;</sub></td><td>or experiments.</td></tr><tr><td colspan="3">Stiffness degradation (cyclic loading of constant sign)</td></tr><tr><td>Axial stiffness of steel</td><td>K<sup>steel</sup><sub>x</sub></td><td>Geometrical &amp; mechanical properties</td></tr><tr><td>Shear stiffness of steel</td><td>K<sup>steel</sup><sub>y</sub></td><td></td></tr><tr><td>Bending stiffness of steel</td><td>K<sup>steel</sup><sub>&theta;</sub></td><td></td></tr><tr><td colspan="3">Stiffness degradation (cyclic loading of alternate sign)</td></tr><tr><td>Geometrical property of the rebars</td><td>c<sub>1</sub></td><td>Found in the technical report of [50]</td></tr><tr><td>Mechanical property of the steel</td><td>c<sub>2</sub></td><td></td></tr></tbody></table>

### 3.2.10. 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 $r_x^{lim}$, $r_y^{lim}$, $r_\theta^{lim}$, such as if one or more of the following conditions are satisfied:

$$r_x = r_x^{lim} \quad , \quad r_y = r_y^{lim} \quad , \quad r_\theta = r_\theta^{lim} \qquad (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:

$$K_x = K_x^{\text{steel}}, \quad K_y = K_y^{\text{steel}}, \quad K_{\theta} = K_{\theta}^{\text{steel}} \qquad (40)$$

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:

$$E = E_0(c_1 + (1 - c_1)e^{-c_2 p^2}) \qquad (41)$$

where $E_0$ is the steel Young modulus, $p$ the plastic strain and $c_1$ and $c_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:

$$\underline{K}_{s}^{el}(i, i) = \underline{K}_{s0}^{el}(i, i)(c_{1}(i) + (1 - c_{1}(i))e^{-c_{2}(i)p^{2}(i)}) \qquad (42)$$

where $\underline{K}_{s0}^{el}$ 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].

### 3.2.11. 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.

## 4. 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.

### 4.1. Enhanced formulation

The adopted enhanced formulation is the one introduced by [7]. In particular, among the four formulations introduced by the authors [7], the N0M0S0 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:

$$\Theta_z(x, t) = \bar{\Theta}_z(x, t) + \bar{\bar{\Theta}}_z(x, t) = \underline{N}^0(x)\underline{U}_e(t) + \underline{M}^0(x)\bar{\underline{U}}_e(t) \qquad (43)$$

where $\Theta_z$ is the total rotation of the section decomposed into a continuous part $\bar{\Theta}_z$ and a discontinuous part $\bar{\bar{\Theta}}_z$, $\underline{U}_e$ the element nodal displacements and $\bar{\underline{U}}_e$ the elementary displacement jumps defined as follows:

$$\underline{U}_e = \underbrace{\begin{bmatrix} U_{xi} \\ U_{yi} \\ \Theta_{zi} \\ U_{xj} \\ U_{yj} \\ \Theta_{zj} \end{bmatrix}}_{\text{nodal displacements}}, \quad \bar{\underline{U}}_e = \underbrace{\begin{bmatrix} \bar{\underline{U}}_x \\ \bar{\underline{U}}_y \\ \bar{\bar{\Theta}}_z \end{bmatrix}}_{\text{elementary displacement jumps}} \qquad (44)$$
<figure><img src="image_9.png" /><figcaption>Fig. 6. Sketch of the plastic hinge mode: initiation of the plastic hinge (up) and fully opened hinge (bottom).</figcaption></figure>

$N^θ = [0 0 N_3^θ 0 0 N_6^θ]$ 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:

$$ \underline{M}^θ = [M_1^θ \quad M_2^θ \quad M_3^θ] = \underbrace{[\bar{M}_1^θ \quad \bar{M}_2^θ \quad \bar{M}_3^θ]}_{\underline{M^θ(x)}} + \underbrace{[0 \quad 0 \quad H_{x_d}^θ]}_{\bar{M}^θ(x)} \quad (45) $$

with $H_{x_d}^θ$ the Heaviside function and $\bar{M}^θ(x)$ 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).

### 4.1.1. Identification procedure

The curvature field $κ_θ(x, t)$ of the Timoshenko beam's section is obtained by derivation of the rotational field of Eq. (43):

$$ κ_θ(x,t) = \bar{κ}_θ(x,t) + \bar{\bar{κ}}_θ(x,t) = \underline{B}^k(x)\underline{U}_e(t) + \underline{G}_r^k(x)\bar{\underline{U}}_e(t) \quad (46) $$

where $\underline{B}^k(x)$ contains the derivatives of the shape functions, with respect to $x$, for the curvature component of the Timoshenko beam (see Eqs. (6), (9)) and $\underline{G}_r^k(x)$ is the matrix of the derivatives of the enhancement interpolation functions, with respect to $x$, of the following form:

$$ \underline{G}_r^k = [\underline{G}_{r1}^k \quad \underline{G}_{r2}^k \quad \underline{G}_{r3}^k] = \underbrace{[\bar{\underline{G}}_{r1}^k \quad \bar{\underline{G}}_{r2}^k \quad \bar{\underline{G}}_{r3}^k]}_{\bar{G}_r(x)} + \underbrace{[0 \quad 0 \quad \delta_{x_d}^θ]}_{\bar{G}_r(x)} \quad (47) $$

with $δ_{x_d}^θ$ the Dirac function.

Using the above decomposition of $\underline{G}_r^k$ into continuous and discontinuous parts, the curvature field corresponding to Eq. (46) is rewritten as:

$$ κ_θ(x,t) = \underbrace{\underline{B}^k(x)\underline{U}_e(t) + \bar{\underline{G}}_r^k(x)\bar{\underline{U}}_e(t)}_{\text{regular part } \bar{κ}_θ(x,t)} + \underbrace{\bar{\underline{G}}_r^k(x)\bar{\underline{U}}_e(t)}_{\text{singular part}} \quad (48) $$

The $\bar{\underline{G}}_r^k$ 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 $L_e$ can be written as (see Fig. 6):

$$ Θ_{zj} = Θ_{zi} + \bar{Θ}_z \quad (49) $$

with $x_d$ the geometrical position of the discontinuity within the finite element (in our case $x_d = L_e/2$). Based on Eq. (49), the rotational jump component is thus given as:

$$ \bar{Θ}_z = Θ_{zj} - Θ_{zi} \quad (50) $$

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:

$$ \bar{κ}_θ(x,t) = \underline{B}^k(x)\underline{U}_e(t) + \bar{\underline{G}}_r^k(x)\bar{\underline{U}}_e(t) = 0 \quad (51) $$

Considering that $\underline{B}^k(x) = [0 \quad 0 \quad -1/L_e \quad 0 \quad 0 \quad 1/L_e]$, at the final state, the regular part of the curvature is expressed as:

$$ \bar{κ}_θ(x,t) = \frac{1}{L_e}(\Theta_{zj} - \Theta_{zi}) + \bar{\underline{G}}_r^k(x)\bar{\underline{U}}_e(t) = 0 \quad (52) $$

By introducing the rotational jump definition of Eq. (50) in Eq. (52), $\bar{\underline{G}}_r^k$ is deduced:

$$ \bar{\underline{G}}_r^k(x) = [\bar{\underline{G}}_{r1}^k \quad \bar{\underline{G}}_{r2}^k \quad \bar{\underline{G}}_{r3}^k] = [0 \quad 0 \quad -\frac{1}{L_e}] \quad (53) $$
<figure><img src="image_10.png" /><figcaption>Fig. 7. Global cohesive model: cohesive moment — rotational jump for cyclic loading of constant (left) and alternate (right) sign.</figcaption></figure>

## 4.2. 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.

### 4.2.1. 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 :

$$ \kappa_{\theta} = \kappa_{\theta}^{act} \qquad (54) $$

where the critical value assigned to the curvature is given as:

$$ \kappa_{\theta}^{act} = \frac{\Theta_{z}^{act}}{100L} \qquad (55) $$

where L is the length of the RC element and the denominator is multiplied by 100 in order to transform $\Theta_{z}^{act}$ given as a percentage in Eq. (56). The identification procedure of $\Theta_{z}^{act}$ is given in Section 4.2.2.

Further increase of the applied loading results in decrease of the cohesive moment $M_{\theta}^{coh}$ and in increase of the rotational jump $\bar{\Theta}_{z}$. 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).

### 4.2.2. Parameter identification

The ultimate rotational capacity $\Theta_{z}^{act}$ (in percent) of Eq. (55) is calculated from the empirical relation provided in [52], based on regression analysis on experimental on RC rectangular sections:

$$ \Theta_{z}^{act}(\%) = 0.52(L/d)^{0.93} \rho^{-0.27} \rho_{w}^{0.48} n_{o}^{-0.48} f_{c}^{-0.15} \qquad (56) $$

with:

* L/d the shear span ratio (L is the length of the RC element and d is the depth of the RC section).

* $\rho$ the normalized steel ratio, defined as $\frac{A_s}{A_c} \frac{f_y}{f_c}$ (where $A_s$ is the area of the longitudinal reinforcement in tension and $A_c$ is the cross-sectional area).

* $f_c$ the concrete compressive strength (in ksi).

* $f_y$ the steel yield limit (in ksi).

* $n_o$ the normalized axial force defined as $P/P_o$, with P the axial force (positive for compression) and $P_o = bdf_c$ (where b is the width and d the depth of the RC element).
* $\rho_w$ the confinement ratio (in percent), defined as $\frac{100A_{sx}}{bs_h}$ (where $A_{sx}$ is the area of the transverse steel parallel to the direction of loading and $s_h$ the spacing of the stirrups).

Finally, the softening modulus $S$ 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].

### 4.3. 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 $M_{\theta}^{coh}$ and the rotational jump $\bar{\Theta}_z$:

$$M_{\theta}^{coh} = M_{\theta}^{u} + S\bar{\Theta}_{z} \qquad (57)$$

where $S$ is the softening modulus (of negative value). The activation criterion is given by:

$$f_c = |M_{\theta}^{coh}| - (M_{\theta}^{u} + S\bar{\Theta}_{z}) \le 0 \qquad (58)$$

where $M_{\theta}^{u}$ corresponds to $\kappa_{\theta}^{act}$ 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 $f_c$. The evolution of the elementary variable of the rotational jump is given by:

$$\dot{\bar{\Theta}}_z = \dot{\bar{\lambda}}\mathrm{sign}(M_{\theta}^{coh}) \qquad (59)$$

where $\dot{\bar{\lambda}}$ 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 $\dot{\bar{\lambda}}$ is positive and thus the failure criterion is satisfied. This implies that:

$$\begin{gathered} \dot{\bar{\lambda}}f_c = 0 \Rightarrow f_c = 0 \Rightarrow \\ |\dot{M}_{\theta}^{coh}| - S\dot{\bar{\Theta}}_z = 0 \Rightarrow |\dot{M}_{\theta}^{coh}| - S\dot{\bar{\lambda}}\mathrm{sign}(M_{\theta}^{coh}) = 0 \Rightarrow \\ \dot{\bar{\lambda}} = \frac{1}{S}|\dot{M}_{\theta}^{coh}| \Rightarrow \dot{\bar{\lambda}} = \frac{1}{S}\dot{M}_{\theta}^{coh}\mathrm{sign}(M_{\theta}^{coh}) \end{gathered} \qquad (60)$$

and therefore the expression for the calculation of the plastic multiplier is known. The cohesive model internal variable $\bar{\Theta}_z$ is computed at every time step (once the activation criterion is satisfied) for every Gauss point at every iteration $k$. The elastic predictor of the current step $n+1$ is equal to the previous step $n$ value:

$$\bar{\Theta}_{z_{n+1}}^{trial} = \bar{\Theta}_{z_n} \qquad (61)$$

and the trial cohesive moment is given by the following expression satisfying the local equilibrium at the discontinuity level within the element:

$$M_{\theta_{n+1}}^{coh,trial} = - \int_{L_c} \bar{G}_{r3}^{\kappa}(x) M_{\theta_{n+1}}^{trial}(x) dx = - \frac{L_c}{2} \sum_{pg=1}^{npg} \bar{G}_{r3}^{\kappa}(x_{pg}) M_{\theta_{n+1}}^{trial}(x_{pg}) w(x_{pg}) \qquad (62)$$

where $npg$ is the number of integration points of the beam element (in our case $npg = 1$) and $w(x_{pg})$ is the integration weight of the Gauss point with coordinate $x_{pg}$ (in our case the Gauss point is in the middle of the FE).

The expression of $M_{\theta_{n+1}}^{trial}(x_{pg})$ is:

$$M_{\theta_{n+1}}^{trial}(x_{pg}) = K_{\theta}(x_{pg}) \left( B^{\kappa} U_{e_{n+1}}^{k-1} + \bar{G}_{r3}^{\kappa}(x_{pg}) \bar{\Theta}_{z_{n+1}}^{trial} - \kappa_{\theta_n}^{pl} \right) \qquad (63)$$

The trial failure criterion is evaluated:

$$f_{e_{n+1}}^{trial} = |M_{\theta_{n+1}}^{coh,trial}| - (M_{\theta}^{u} + S\bar{\Theta}_{z_{n+1}}^{trial}) \le 0 \qquad (64)$$

If this condition (64) holds the discontinuity is not activated and thus the variables of the cohesive model are updated as:

$$\begin{gathered} \bar{\Theta}_{z_{n+1}} = \bar{\Theta}_{z_{n+1}}^{trial} \\ M_{\theta_{n+1}}^{coh} = M_{\theta_{n+1}}^{coh,trial} \end{gathered} \qquad (65)$$
If discontinuity is activated the rotational jump at the current step is given by:

$$ \bar{\bar{\Theta}}_{n+1} = \bar{\bar{\Theta}}_n + \frac{\dot{\bar{\lambda}}\mathrm{sign}(M_{\theta_{n+1}}^{\mathrm{coh,trial}})}{\bar{\bar{\lambda}}} \qquad (66) $$

The updated cohesive force at the step n + 1 can be expressed as:

$$ \begin{aligned} M_{\theta_{n+1}}^{\mathrm{coh}} &= -\int_{L_e} \bar{\bar{G}}_{\nu3}^{\kappa}(x) M_{\theta_{n+1}}(x) dx \\ &= -\int_{L_e} \bar{\bar{G}}_{\nu3}^{\kappa}(x) K_{\theta} \left( \underline{\underline{B}}^{\kappa}(x) \underline{U}_{\epsilon_{n+1}} + \bar{\bar{G}}_{\nu3}^{\kappa}(x) \bar{\bar{\Theta}}_{z_{n+1}} - \kappa_{\theta_n}^{\mathrm{pl}} \right) dx \\ &= -\underbrace{\int_{L_e} \bar{\bar{G}}_{\nu3}^{\kappa}(x) K_{\theta} \left( \underline{\underline{B}}^{\kappa}(x) \underline{U}_{\epsilon_{n+1}} + \bar{\bar{G}}_{\nu3}^{\kappa}(x) \bar{\bar{\Theta}}_{z_n} - \kappa_{\theta_n}^{\mathrm{pl}} \right) dx}_{M_{\theta_{n+1}}^{\mathrm{coh,trial}}} - \underbrace{\int_{L_e} \bar{\bar{G}}_{\nu}^{\kappa}(x) K_{\theta} \bar{\bar{G}}_{\nu3}^{\kappa}(x) dx \Delta \bar{\bar{\Theta}}_z}_{K_m} \end{aligned} \qquad (67) $$

and thus:

$$ M_{\theta_{n+1}}^{\mathrm{coh}} = M_{\theta_{n+1}}^{\mathrm{coh,trial}} - K_m \frac{\dot{\bar{\lambda}}\mathrm{sign}(M_{\theta_{n+1}}^{\mathrm{coh,trial}})}{\bar{\bar{\lambda}}} \qquad (68) $$

The above expression allows for the calculation of the plastic multiplier:

$$ \frac{\dot{\bar{\lambda}}}{\bar{\bar{\lambda}}} = \begin{cases} \frac{f_{\epsilon_{n+1}}^{\mathrm{trial}}}{K_m + S} & \text{if } |S\bar{\bar{\Theta}}_{z_n}| < M_{\theta}^u \\ \frac{|M_{\theta_{n+1}}^{\mathrm{coh,trial}}|}{K_m} & \text{if } |S\bar{\bar{\Theta}}_{z_n}| = M_{\theta}^u \end{cases} \qquad (69) $$

and the cohesive tangent modulus is defined as:

$$ \frac{d M_{\theta_{n+1}}^{\mathrm{coh}}}{d \bar{\bar{\Theta}}_z} = \begin{cases} \text{not defined if } \dot{\bar{\lambda}} = 0 \\ S & \text{if } \dot{\bar{\lambda}} > 0 \text{ and } |S\bar{\bar{\Theta}}_{z_{n+1}}| < M_{\theta}^u \\ 0 & \text{if } \dot{\bar{\lambda}} > 0 \text{ and } |S\bar{\bar{\Theta}}_{z_{n+1}}| = M_{\theta}^u \end{cases} \qquad (70) $$

# 5. 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.

## 5.1. 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.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Two story RC frame: geometry and reinforcement details.
*Source: Adapted from [55].</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Two story RC frame: FE mesh and boundary conditions, (left) node numbers, (right) element numbers.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Two story RC frame: Numerical global response (left) and element's contribution (right) (monotonic loading).</figcaption></figure>
<table><caption>**Table 3**
Model parameters for the two story RC frame.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Value</th><th>SI unit</th></tr></thead><tbody><tr><td colspan="4">Linear elastic behavior</td></tr><tr><td>Axial stiffness</td><td>K<sub>x</sub></td><td>3.92e09</td><td>N</td></tr><tr><td>Shear stiffness</td><td>K<sub>y</sub></td><td>1.62e09</td><td>N</td></tr><tr><td>Bending stiffness</td><td>K<sub>θ</sub></td><td>5.67e07</td><td>N m<sup>2</sup></td></tr><tr><td>Axial elastic limit</td><td>r<sub>x0</sub></td><td>0.5</td><td>–</td></tr><tr><td>Shear elastic limit</td><td>r<sub>y0</sub></td><td>0.5</td><td>–</td></tr><tr><td>Bending elastic limit</td><td>r<sub>θ0</sub></td><td>0.5</td><td>–</td></tr><tr><td colspan="4">Standardization</td></tr><tr><td>Maximum axial tensile force</td><td>F<sup>t</sup><sub>x<sub>max</sub></sub></td><td>1.50e06</td><td>N</td></tr><tr><td>Maximum axial compressive force</td><td>F<sup>s</sup><sub>x<sub>max</sub></sub></td><td>-4.66e06</td><td>N</td></tr><tr><td>Shear shift parameter</td><td>F<sup>s</sup><sub>y</sub></td><td>4.14e05</td><td>N</td></tr><tr><td>Bending shift parameter</td><td>M<sup>*</sup><sub>θ</sub></td><td>2.91e05</td><td>N m</td></tr><tr><td colspan="4">Hardening evolution laws</td></tr><tr><td>Axial rate of hardening</td><td>a<sub>x</sub></td><td>500</td><td>–</td></tr><tr><td>Shear rate of hardening</td><td>a<sub>y</sub></td><td>250</td><td>–</td></tr><tr><td>Bending rate of hardening</td><td>a<sub>θ</sub></td><td>250</td><td>–</td></tr><tr><td colspan="4">Stiffness degradation<br>(cyclic loading of constant sign)</td></tr><tr><td>Axial stiffness of steel</td><td>K<sup>steel</sup><sub>x</sub></td><td>4.84e08</td><td>N</td></tr><tr><td>Shear stiffness of steel</td><td>K<sup>steel</sup><sub>y</sub></td><td>1.86e08</td><td>N</td></tr><tr><td>Bending stiffness of steel</td><td>K<sup>steel</sup><sub>θ</sub></td><td>1.09e07</td><td>N m<sup>2</sup></td></tr><tr><td colspan="4">Cohesive model</td></tr><tr><td>Concrete compressive strength</td><td>f<sub>c</sub></td><td>4.35</td><td>ksi</td></tr><tr><td>Steel yield limit</td><td>f<sub>y</sub></td><td>60.63</td><td>ksi</td></tr><tr><td>Normalized axial force</td><td>n<sub>θ</sub></td><td>0.194</td><td>–</td></tr><tr><td>Normalized steel ratio</td><td>ρ</td><td>0.146</td><td>–</td></tr><tr><td>Depth of the section</td><td>d</td><td>0.40</td><td>m</td></tr><tr><td>Shear span ratio (beam & column)</td><td>L/d</td><td>8.75 & 5</td><td>–</td></tr><tr><td>Confinement ratio in percent</td><td>ρ<sub>u</sub></td><td>2.51%</td><td>–</td></tr><tr><td>Softening modulus</td><td>S</td><td>-3.97e06</td><td>N m</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Two story RC frame: Numerical vs. experimental global response (cyclic loading with constant sign).</figcaption></figure>

## 5.2. RC simply supported beam

### 5.2.1. 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)
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** RC simply supported beam: geometry and reinforcement details (left), vertical displacement time history (right).
Source: Adapted from [57].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** RC simply supported beam: FE mesh and boundary conditions, node (left) and element numbers (right).</figcaption></figure>

<table><caption>**Table 4**
Model parameters for the RC simply supported beam.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Value</th><th>SI unit</th></tr></thead><tbody><tr><td colspan="4">Linear elastic behavior</td></tr><tr><td>Axial stiffness</td><td>K<sub>x</sub></td><td>1.02e09</td><td>N</td></tr><tr><td>Shear stiffness</td><td>K<sub>y</sub></td><td>4.21e08</td><td>N</td></tr><tr><td>Bending stiffness</td><td>K<sub>θ</sub></td><td>3.66e06</td><td>N m<sup>2</sup></td></tr><tr><td>Axial elastic limit</td><td>r<sub>x0</sub></td><td>0.7</td><td>–</td></tr><tr><td>Shear elastic limit</td><td>r<sub>y0</sub></td><td>0.7</td><td>–</td></tr><tr><td>Bending elastic limit</td><td>r<sub>θ0</sub></td><td>0.7</td><td>–</td></tr><tr><td colspan="4">Standardization</td></tr><tr><td>Maximum axial tensile force</td><td>F<sup>t</sup><sub>xmax</sub></td><td>2.80e05</td><td>N</td></tr><tr><td>Maximum axial compressive force</td><td>F<sup>c</sup><sub>xmax</sub></td><td>-1.13e06</td><td>N</td></tr><tr><td>Shear shift parameter</td><td>F<sup>s</sup><sub>y</sub></td><td>8.67e04</td><td>N</td></tr><tr><td>Bending shift parameter</td><td>M<sup>*</sup><sub>θ</sub></td><td>3.86e04</td><td>N m</td></tr><tr><td colspan="4">Hardening evolution laws</td></tr><tr><td>Axial rate of hardening</td><td>a<sub>x</sub></td><td>500</td><td>–</td></tr><tr><td>Shear rate of hardening</td><td>a<sub>y</sub></td><td>250</td><td>–</td></tr><tr><td>Bending rate of hardening</td><td>a<sub>θ</sub></td><td>250</td><td>–</td></tr><tr><td colspan="4">Stiffness degradation<br>(cyclic loading of alternate sign)</td></tr><tr><td>Geometrical property of the rebars</td><td>c<sub>1</sub></td><td>0.3</td><td>–</td></tr><tr><td>Mechanical property of the steel</td><td>c<sub>2</sub></td><td>620</td><td>–</td></tr><tr><td colspan="4">Cohesive model</td></tr><tr><td>Concrete compressive strength</td><td>f<sub>c</sub></td><td>4.49</td><td>ksi</td></tr><tr><td>Steel yield limit</td><td>f<sub>y</sub></td><td>65.27</td><td>ksi</td></tr><tr><td>Normalized axial force</td><td>n<sub>θ</sub></td><td>0.001</td><td>–</td></tr><tr><td>Normalized steel ratio</td><td>ρ</td><td>0.110</td><td>–</td></tr><tr><td>Depth of the section</td><td>d</td><td>0.20</td><td>m</td></tr><tr><td>Shear span ratio</td><td>L/d</td><td>4.25</td><td>–</td></tr><tr><td>Confinement ratio in percent</td><td>ρ<sub>w</sub></td><td>1.064%</td><td>–</td></tr><tr><td>Softening modulus</td><td>S</td><td>-2.56e05</td><td>N m</td></tr></tbody></table>

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
<figure><img src="image_17.png" /><figcaption>Fig. 14. RC simply supported beam: Numerical vs. experimental (left) and numerical (calibrated) vs. experimental load/time response (right).</figcaption></figure>

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 $c1$ 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.

## 6. 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.

## CRediT authorship contribution statement

**Androniki-Anna Doulgeroglou:** Writing – original draft, Validation, Software, Investigation. **Panagiotis Kotronis:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Giulio Sciarra:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Catherine Bouillon:** Funding acquisition.

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

## Data availability

Data may be available upon request.

## Acknowledgments

The authors would like to acknowledge the financial support of Groupe Essor, France for funding the Cifre PhD of the first author.

## References

[1] H. Hippola, K. Wijesundara, R. Nascimbene, Response of shear critical reinforced concrete frames and walls under monotonic loading, Eng. Struct. 251 (2022) 113483.

[2] E. Brunesi, S. Peloso, R. Pinho, R. Nascimbene, Cyclic tensile testing of a three-way panel connection for precast wall-slab-wall structures, Struct. Concr. 20 (4) (2019) 1307-1315.

[3] R. Nascimbene, Numerical model of a reinforced concrete building: Earthquake analysis and experimental validation, Period. Polytech. Civ. Eng. 59 (4) (2015) 521-530.
[4] E.N. Dvorkin, A.M. Cuitino, G. Gioia, Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions, Int. J. Numer. Methods Eng. 30 (3) (1990) 541–564.

[5] R. Nova, L. Montrasio, Settlements of shallow foundations on sand, Géotechnique 41 (2) (1991) 243–256.

[6] F. Armero, D. Ehrlich, An analysis of strain localization and wave propagation in plastic models of beams at failure, Comput. Methods Appl. Mech. Engrg. 193 (30–32) (2004) 3129–3171.

[7] D. Ehrlich, F. Armero, Finite element methods for the analysis of softening plastic hinges in beams and frames, Comput. Mech. 35 (4) (2005) 237–264.

[8] F. Armero, D. Ehrlich, Numerical modeling of softening hinges in thin Euler–Bernoulli beams, Comput. Struct. 84 (10–11) (2006) 641–656.

[9] C. Cremer, A. Pecker, L. Davenne, Elaboration of a SSI macro-element with uplift of shallow foundation, in: Implications of Recent Earthquakes on Seismic Risk, World Scientific, 2000, pp. 127–138.

[10] C. Chatzigogos, A. Pecker, J. Salencon, Macroelement modeling of shallow foundations, Soil Dyn. Earthq. Eng. 29 (5) (2009) 765–781.

[11] S. Grange, P. Kotronis, J. Mazars, A macro-element to simulate 3D soil–structure interaction considering plasticity and uplift, Int. J. Solids Struct. 46 (20) (2009) 3651–3663.

[12] Z. Li, P. Kotronis, S. Escoffier, Numerical study of the 3D failure envelope of a single pile in sand, Comput. Geotech. 62 (2014) 11–26.

[13] Z. Li, P. Kotronis, S. Escoffier, C. Tamagnini, A hypoplastic macroelement for single vertical piles in sand subject to three-dimensional loading conditions, Acta Geotech. 11 (2) (2016) 373–390.

[14] Z. Li, P. Kotronis, S. Escoffier, C. Tamagnini, A hypoplastic macroelement formulation for single batter piles in sand, Int. J. Numer. Anal. Methods Geomech. 42 (12) (2018) 1346–1365.

[15] J. Perez-Herreros, *Dynamic Soil-Structure Interaction of Pile Foundations: Experimental and Numerical Study (Ph.D. thesis)*, École centrale de Nantes, France, 2020, https://theses.hal.science/tel-02972312.

[16] Z. Jin, Z.-Y. Yin, P. Kotronis, Z. Li, C. Tamagnini, A hypoplastic macroelement model for a caisson foundation in sand under monotonic and cyclic loadings, Mar. Struct. 66 (2019) 16–26.

[17] S.M. Elachachi, Sur l'élaboration d'une Méthode Simplifiée d'Analyse des Structures de Génie Civil par Macro-Elements Adaptés Aux Constructions Composites et Endommageables (Ph.D. thesis), Université Paris 6, Paris, France, 1992.

[18] L. Davenne, C. Brenet, Macro-Éléments de Poutres en Béton Armé, Tech. Rep., Rapport interne LMT - Cachan, num. 210, juin, France, 1998.

[19] S. Tajiri, H. Shiohara, F. Kusuhara, A new macroelement of reinforced concrete beam column joint for elasto-plastic plane frame analysis, in: Proc. 8th National Conf. Earthquake Engineering, Citeseer, 2006, paper no. 674.

[20] M. Hemsas, S.M. Elachachi, D. Breysse, Modélisation par macroéléments du comportement non linéaire des voiles en béton armé, Eur. J. Environ. Civ. Eng. 13 (5) (2009) 615–640.

[21] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. (ISSN: 0029-5981) 29 (8) (1990) 1595–1638, http://dx.doi.org/10.1002/nme.1620290802, URL http://doi.wiley.com/10.1002/nme.1620290802,

[22] B.-H. Pham, L. Davenne, D. Brancherie, A. Ibrahimbegovic, Stress resultant model for ultimate load design of reinforced-concrete frames: combined axial force and bending moment, Comput. Concr. 7 (4) (2010) 303–315.

[23] B. Pham, D. Brancherie, L. Davenne, A. Ibrahimbegovic, Stress resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates, Comput. Mech. 51 (3) (2013) 347–360.

[24] M. Jukić, B. Brank, A. Ibrahimbegović, Embedded discontinuity finite element formulation for failure analysis of planar reinforced concrete beams and frames, Eng. Struct. 50 (2013) 115–125.

[25] M. Jukić, Finite Elements for Modeling of Localized Failure in Reinforced Concrete (Ph.D. thesis), École normale supérieure de Cachan-ENS Cachan; Univerza v Ljubljani. Fakulteta za gradbeništvo in geodezijo, France, Slovenia, 2013.

[26] M. Jukić, B. Brank, A. Ibrahimbegovic, Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity, Eng. Struct. 75 (2014) 507–527.

[27] N. Bui, M. Ngo, M. Nikolic, D. Brancherie, A. Ibrahimbegovic, Enriched timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames, Comput. Struct. 143 (2014) 9–18.

[28] I. Bitar, Modélisation de la Rupture Dans les Structures en Béton Armé par des Éléments Finis Poutres Généralisées et Multifibres (Ph.D. thesis), École centrale de Nantes, France, 2017, https://theses.hal.science/tel-01755118.

[29] I. Bitar, P. Kotronis, N. Benkemoun, S. Grange, A generalized timoshenko beam with embedded rotation discontinuity, Finite Elem. Anal. Des. 150 (2018) 34–50.

[30] I. Bitar, N. Benkemoun, P. Kotronis, S. Grange, A multifiber timoshenko beam with embedded discontinuities, Eng. Fract. Mech. 214 (2019) 339–364.

[31] D. Caillerie, P. Kotronis, R. Cybulski, A Timoshenko finite element straight beam with internal degrees of freedom, Int. J. Numer. Anal. Methods Geomech. 39 (16) (2015) 1753–1773.

[32] A.-A. Doulgeroglou, P. Kotronis, G. Sciarra, C. Bouillon, 3D interaction diagrams for symmetrically reinforced concrete square sections with various reinforcement ratios, Eng. Struct. 262 (2022) 114272.

[33] A.-A. Doulgeroglou, A Novel Macroelement to Assess the Vulnerability of Reinforced Concrete Frame Stuctures Under Severe Dynamic Loadings (Ph.D. thesis), École centrale de Nantes, France, 2022, https://theses.hal.science/tel-03952840.

[34] P. Pegon, A Timoshenko Simple Beam Element in Castem 2000, Special Publication Nr. I. 94.04. Italy: Applied Mechanics Unit, Institute for Safety Technology, Joint Research Centre, Commission of the European Communities, I-21020 ISPRA (VA), 1994.

[35] J. Guedes, P. Pegon, A.V. Pinto, A Fibre/Timoshenko Beam Element in Castem 2000, Special publication Nr. I.94.31, Applied Mechanics Unit, Institute for Safety Technology, Joint Research Centre, Commission of the European Communities, I-21020 ISPRA (VA), Italy, Ispra Italy, 1994.

[36] J. Donea, L.G. Lamain, A modified representation of transverse shear in C0 quadrilateral plate elements, Comput. Methods Appl. Mech. Engrg. 63 (2) (1987) 183–207.

[37] I. Bitar, S. Grange, P. Kotronis, N. Benkemoun, A comparison of displacement-based timoshenko multi-fiber beams finite element formulations and elasto-plastic applications, Eur. J. Environ. Civ. Eng. 22 (4) (2018) 464–490, http://dx.doi.org/10.1080/19648189.2016.1210031, arXiv:https://doi.org/10.1080/19648189.2016.1210031.

[38] E. Spacone, F.C. Filippou, F.F. Taucer, Fibre beam–column model for non-linear analysis of R/C frames: Part I. Formulation, Earthq. Eng. Struct. Dyn. 25 (7) (1996) 711–725.

[39] CEN, Eurocode 2: Design of Concrete Structures: Part 1-1: General Rules and Rules for Buildings, EN 1992-1-1, European Committee for Standardization, Brussels, 2004.

[40] K.N. Rahal, Shear strength of reinforced concrete: Part II—Beams subjected to shear, bending moment, and axial load, Struct. J. 97 (2) (2000) 219–224.

[41] F.J. Vecchio, M.P. Collins, The modified compression-field theory for reinforced concrete elements subjected to shear, ACI J. 83 (2) (1986) 219–231.

[42] J.-H. Prevost, Anisotropic undrained stress-strain behaviour of clays-closure, J. Geotech. Geoenviron. Eng. 106 (ASCE 13942) (1980).

[43] C. Cremer, A. Pecker, L. Davenne, Cyclic macro-element for soil–structure interaction: material and geometrical non-linearities, Int. J. Numer. Anal. Methods Geomech. 25 (13) (2001) 1257–1284.

[44] S. Grange, P. Kotronis, J. Mazars, A macro-element for a circular foundation to simulate 3D soil–structure interaction, Int. J. Numer. Anal. Methods Geomech. 32 (10) (2008) 1205–1227.

[45] J.C. Simo, T.J. Hughes, Computational Inelasticity, vol. 7, Springer Science & Business Media, 2006.
[46] S. Grange, Modélisation Simplifiée 3D de l'Interaction Sol-Structure: Application au Génie Parasismique (Ph.D. thesis), Institut National Polytechnique de Grenoble - INPG, France, 2008, https://theses.hal.science/tel-00306842.

[47] R. Charlier, Approche Unifiée de Quelques Problèmes Non Linéaires de Mécanique des Milieux Continus par la Méthode des Éléments Finis (Ph.D. thesis), Université de Liège, Liège, Belgique, 1987.

[48] J.-D. Barnichon, Finite Element Modelling in Structural and Petroleum Geology (Ph.D. thesis), Université de Liège, Liège, Belgium, 1998.

[49] F. Collin, Couplages Thermo-Hydro-Mécaniques Dans les Sols et les Roches Tendres Partiellement Saturés (Ph.D. thesis), Université de Liège, Liège, Belgique, 2003.

[50] J.-L. Fléjou, Relations de Comportement Non Linéaires 1D, Technical Report, Code Aster, EDF R & D, 2018.

[51] M. Menegotto, P. Pinto, Method of analysis for cyclically loaded reinforced concrete plane frames including changes in geometry and non-elastic of elements under combined normal force and bending, in: Proceedings, IABSE symposium, 1973, pp. 20-102.

[52] Y.-J. Park, A.H. Ang, Y.K. Wen, Damage-limiting aseismic design of buildings, Earthq. Spectra 3 (1) (1987) 1-26.

[53] C.B. Haselton, Assessing Seismic Collapse Safety of Modern Reinforced Concrete Moment Frame Buildings (Ph.D. thesis), Stanford University, USA, 2006.

[54] R.I. Borja, Plasticity: Modeling & Computation, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN: 978-3-642-38547-6, 2013, http://dx.doi.org/10.1007/978-3-642-38547-6_3.

[55] F.J. Vecchio, M.B. Emara, Shear deformations in reinforced concrete frames, ACI Struct. J. 89 (1) (1992) 46-56.

[56] MATLAB, Version: 9.13.0 (r2022b), The MathWorks Inc., Natick, Massachusetts, United States, 2022, https://www.mathworks.com.

[57] G. Lebon, Analyse de l'Endommagement des Structures de Génie Civil: Techniques de Sous-Structuration Hybride Couplées à un Modèle D'endommagement Anisotrope (Ph.D. thesis), École normale supérieure de Cachan-ENS Cachan, France, 2011.