# Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications

Maryam Trad, Ibrahim Bitar, Stéphane Grange, Benjamin Richard

## ► To cite this version:

Maryam Trad, Ibrahim Bitar, Stéphane Grange, Benjamin Richard. Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications. Finite Elements in Analysis and Design, 2024, 237, pp.104170. 10.1016/j.finel.2024.104170. hal-04574181

HAL Id: hal-04574181

https://hal.science/hal-04574181

Submitted on 14 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications

Maryam Trad, Ibrahim Bitar, Stéphane Grange, Benjamin Richard

## ► To cite this version:

Maryam Trad, Ibrahim Bitar, Stéphane Grange, Benjamin Richard. Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications. Finite Elements in Analysis and Design, 2024, 237, pp.104170. 10.1016/j.finel.2024.104170. hal-04574181

HAL Id: hal-04574181

https://hal.science/hal-04574181v1

Submitted on 14 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

# Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications

Maryam Trad<sup>a,b,*</sup>, Ibrahim Bitar<sup>a</sup>, Stéphane Grange<sup>b</sup>, Benjamin Richard<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Steel-concrete interface model

Kinematic enhancement

Bond model

Reinforced concrete structures

Beam finite element

Plate finite element

## ABSTRACT

Considering the interaction between concrete and steel reinforcement in numerical simulations of reinforced concrete structures is crucial for accurately capturing the concrete cracking process. This is particularly interesting when studying structures fulfilling functions that go beyond their simple mechanical resistance, such as waterproofing functions. While three-dimensional (3D) volumetric finite element modeling offers detailed insights into structural behavior, its computational intensity becomes prohibitive for large-scale structures. In such contexts, adopting beam and plate elements formulations proves computationally more efficient, due to their reduced number of degrees of freedom. This paper presents a kinematic enhancement technique designed to integrate steel-concrete interface behavior into beam and plate finite element formulations. The approach combines classical beam or plate elements representing concrete behavior, conventional beam or truss elements modeling steel reinforcement, and the incorporation of bond stresses at the interface. The paper provides comprehensive explanations of this enhancement technique along with a curated selection of numerical validation and application examples. These examples are supplemented by a comparison with experimental data, illustrating the efficiency of the proposed enhancement approach.

## 1. Introduction

The mechanical behavior of reinforced concrete relies on a stress transfer between steel and concrete, through their interface. As soon as the first cracks appear in concrete, the steel tensile properties come into action, provided that the steel-concrete interface transmits the corresponding internal forces. Thus, cracking in reinforced concrete structures is significantly influenced by the stress distribution at the steel-concrete interface. Taking this interface into account in numerical modeling has significant importance on the cracking process of reinforced concrete structures and the spatial distribution of cracks. This is particularly important in case of civil engineering structures exhibiting waterproofing requirements [1].

Various experimental setups are used in the literature to understand and characterize the interface behavior. The most common test types lie in pulling a steel bar out of a concrete specimen [2–7]. Nonetheless, push-in type tests can also be used to represent a compressed concrete state, similarly to its state in pre-stressed concrete structures [8]. In addition, some other tests aim to study the interface in bending elements [9,10].

<footnote><sup>a</sup> Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-EXP/SES/LMAPS, B.P. 17 - 92262 Fontenay-aux-Roses Cedex, France</footnote>

<footnote><sup>b</sup> INSA Lyon, GEOMAS, UR7495, 69621 Villeurbanne, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: maryamtrad100@hotmail.com (M. Trad).</footnote>
<figure><img src="image_2.png" /></figure>

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

# Enhanced beams and plates models incorporating the steel-concrete interface behavior for large-scale reinforced concrete structural applications

Maryam Trad<sup>a,b,*</sup>, Ibrahim Bitar<sup>a</sup>, Stéphane Grange<sup>b</sup>, Benjamin Richard<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Steel-concrete interface model

Kinematic enhancement

Bond model

Reinforced concrete structures

Beam finite element

Plate finite element

## ABSTRACT

Considering the interaction between concrete and steel reinforcement in numerical simulations of reinforced concrete structures is crucial for accurately capturing the concrete cracking process. This is particularly interesting when studying structures fulfilling functions that go beyond their simple mechanical resistance, such as waterproofing functions. While three-dimensional (3D) volumetric finite element modeling offers detailed insights into structural behavior, its computational intensity becomes prohibitive for large-scale structures. In such contexts, adopting beam and plate elements formulations proves computationally more efficient, due to their reduced number of degrees of freedom. This paper presents a kinematic enhancement technique designed to integrate steel-concrete interface behavior into beam and plate finite element formulations. The approach combines classical beam or plate elements representing concrete behavior, conventional beam or truss elements modeling steel reinforcement, and the incorporation of bond stresses at the interface. The paper provides comprehensive explanations of this enhancement technique along with a curated selection of numerical validation and application examples. These examples are supplemented by a comparison with experimental data, illustrating the efficiency of the proposed enhancement approach.

## 1. Introduction

The mechanical behavior of reinforced concrete relies on a stress transfer between steel and concrete, through their interface. As soon as the first cracks appear in concrete, the steel tensile properties come into action, provided that the steel-concrete interface transmits the corresponding internal forces. Thus, cracking in reinforced concrete structures is significantly influenced by the stress distribution at the steel-concrete interface. Taking this interface into account in numerical modeling has significant importance on the cracking process of reinforced concrete structures and the spatial distribution of cracks. This is particularly important in case of civil engineering structures exhibiting waterproofing requirements [1].

Various experimental setups are used in the literature to understand and characterize the interface behavior. The most common test types lie in pulling a steel bar out of a concrete specimen [2–7]. Nonetheless, push-in type tests can also be used to represent a compressed concrete state, similarly to its state in pre-stressed concrete structures [8]. In addition, some other tests aim to study the interface in bending elements [9,10].

<footnote><sup>a</sup> Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-EXP/SES/LMAPS, B.P. 17 - 92262 Fontenay-aux-Roses Cedex, France</footnote>

<footnote><sup>b</sup> INSA Lyon, GEOMAS, UR7495, 69621 Villeurbanne, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: maryamtrad100@hotmail.com (M. Trad).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Beam configuration.</figcaption></figure>

In parallel, several numerical models are proposed within different frameworks to describe the interface between steel and concrete (spring elements [2,11], 2D and 3D elements [12], joint elements [13–16], interface modeling within 2D and 3D homogenized reinforced concrete elements [17], etc.). A large variety of these models is incorporated in 2D and 3D finite elements analysis. A description of the steel-concrete interface at the level of a large-scale structure is highly time-consuming when considering these approaches. In this case, employing semi-global beam and plate finite elements can be computationally more efficient. Several studies have proposed different development strategies to include the steel-concrete interface behavior within beam and plate elements. Some models consider the effect of bonding implicitly by enhancing the steel fiber strain component within fiber-based beam elements [18–20] or incorporating the interface behavior within homogenized reinforced concrete plate elements [21,22]. Using these techniques may be challenging to provide a detailed localized description of the concrete cracking around the steel reinforcement [21,23]. On the other hand, explicit bond-slip behavior modeling is possible with other modeling strategies, where steel displacement degrees of freedom are defined independently of the concrete degrees of freedom, as conceived in [24–26] within a fiber beams framework. Moreover, multiscale modeling strategies have proven their numerical efficiency for a large variety of applications such as concrete damage [27] and steel-concrete composite components [28–31].

This paper proposes an extension of the displacement-based fiber element of [25] by creating a multiscale enhanced element formulation composed of an assembly of different finite elements representing: concrete (beam of plate concrete finite elements), steel (bar or beam steel elements), and steel-concrete bond stresses. The number of assembled concrete and steel elements is a user choice. Similarly to [25], independent degrees of freedom for steel and concrete are integrated within a single finite element. In contrast, the proposed enhancement technique in this work adopts a sub-structured approach, combining existing bar, beam, and plate elements from the literature into a unified enhanced element. This multiscale method facilitates the efficient definition of different reinforced concrete beam or plate elements configurations that consider a nonlinear steel-concrete behavior, while one single specific configuration of a reinforced concrete beam element is defined in [25].

This paper is organized as follows. Section 2 recalls the principles of classical beam and plate theories and details the formulation of the proposed enhancement technique. Sections 3 and 4 present model verification and application examples. The presented application of a four-point flexural beam test presents a comparison with experimental results to demonstrate the representativity of the proposed enhancement approach. The paper ends with conclusions and perspectives in Section 5.

## 2. Theoretical background

This section outlines the theoretical framework of classical beam and plate finite elements and introduces a novel enhancement approach. This approach involves an enhanced element that integrates contributions from standard finite elements representing concrete and steel, along with bond stresses between them.

### 2.1. Beam elements

The literature offers diverse beam models, with varying levels of complexity [32]. Kinematic assumptions relate displacements of point M of the beam to point G, the projection of M on the beam’s neutral axis (see Fig. 1):

$$ u_M(x) = u_G(x) + \theta(x) \wedge \overline{GM} = \begin{bmatrix} u - y\theta \\ v \end{bmatrix} \qquad (1) $$

where u and v are the displacements of point G in the beam longitudinal and normal directions x and y, respectively. $u_M(x)$ and $u_G(x)$ are the displacement vectors of points M and G. $\theta$ is the rotation of the beam section.

#### 2.1.1. Weak formulation

Let B be a beam of a length l and a cross-sectional area S. Solving the beam boundary value problem lies in finding the generalized displacement vector $u = [u \ v \ \theta]^T$ at each longitudinal position $0 \le x \le l$ of the beam, where T stands for the transpose operator.

A generalized strain vector is defined at each beam section as follows:

$$ \epsilon = \begin{bmatrix} \frac{\partial u(x)}{\partial x} & \beta_y & \frac{\partial \theta(x)}{\partial x} \end{bmatrix}^T \qquad (2) $$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Beam configuration.</figcaption></figure>

In parallel, several numerical models are proposed within different frameworks to describe the interface between steel and concrete (spring elements [2,11], 2D and 3D elements [12], joint elements [13–16], interface modeling within 2D and 3D homogenized reinforced concrete elements [17], etc.). A large variety of these models is incorporated in 2D and 3D finite elements analysis. A description of the steel-concrete interface at the level of a large-scale structure is highly time-consuming when considering these approaches. In this case, employing semi-global beam and plate finite elements can be computationally more efficient. Several studies have proposed different development strategies to include the steel-concrete interface behavior within beam and plate elements. Some models consider the effect of bonding implicitly by enhancing the steel fiber strain component within fiber-based beam elements [18–20] or incorporating the interface behavior within homogenized reinforced concrete plate elements [21,22]. Using these techniques may be challenging to provide a detailed localized description of the concrete cracking around the steel reinforcement [21,23]. On the other hand, explicit bond-slip behavior modeling is possible with other modeling strategies, where steel displacement degrees of freedom are defined independently of the concrete degrees of freedom, as conceived in [24–26] within a fiber beams framework. Moreover, multiscale modeling strategies have proven their numerical efficiency for a large variety of applications such as concrete damage [27] and steel-concrete composite components [28–31].

This paper proposes an extension of the displacement-based fiber element of [25] by creating a multiscale enhanced element formulation composed of an assembly of different finite elements representing: concrete (beam of plate concrete finite elements), steel (bar or beam steel elements), and steel-concrete bond stresses. The number of assembled concrete and steel elements is a user choice. Similarly to [25], independent degrees of freedom for steel and concrete are integrated within a single finite element. In contrast, the proposed enhancement technique in this work adopts a sub-structured approach, combining existing bar, beam, and plate elements from the literature into a unified enhanced element. This multiscale method facilitates the efficient definition of different reinforced concrete beam or plate elements configurations that consider a nonlinear steel-concrete behavior, while one single specific configuration of a reinforced concrete beam element is defined in [25].

This paper is organized as follows. Section 2 recalls the principles of classical beam and plate theories and details the formulation of the proposed enhancement technique. Sections 3 and 4 present model verification and application examples. The presented application of a four-point flexural beam test presents a comparison with experimental results to demonstrate the representativity of the proposed enhancement approach. The paper ends with conclusions and perspectives in Section 5.

## 2. Theoretical background

This section outlines the theoretical framework of classical beam and plate finite elements and introduces a novel enhancement approach. This approach involves an enhanced element that integrates contributions from standard finite elements representing concrete and steel, along with bond stresses between them.

### 2.1. Beam elements

The literature offers diverse beam models, with varying levels of complexity [32]. Kinematic assumptions relate displacements of point M of the beam to point G, the projection of M on the beam’s neutral axis (see Fig. 1):

$$ u_M(x) = u_G(x) + \theta(x) \wedge \overline{GM} = \begin{bmatrix} u - y\theta \\ v \end{bmatrix} \qquad (1) $$

where u and v are the displacements of point G in the beam longitudinal and normal directions x and y, respectively. $u_M(x)$ and $u_G(x)$ are the displacement vectors of points M and G. $\theta$ is the rotation of the beam section.

#### 2.1.1. Weak formulation

Let B be a beam of a length l and a cross-sectional area S. Solving the beam boundary value problem lies in finding the generalized displacement vector $u = [u \ v \ \theta]^T$ at each longitudinal position $0 \le x \le l$ of the beam, where T stands for the transpose operator.

A generalized strain vector is defined at each beam section as follows:

$$ \epsilon = \begin{bmatrix} \frac{\partial u(x)}{\partial x} & \beta_y & \frac{\partial \theta(x)}{\partial x} \end{bmatrix}^T \qquad (2) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Plate configuration.</figcaption></figure>

knowing that the shear strain $\beta_y$ is considered within the Timoshenko beam theory and neglected when following the Euler-Bernoulli theory. A generalized forces and moments vector $\mathbf{F} = [F_x \ F_y \ M_z]^T$ is defined for each beam section. This section is subjected to three types of loading: a normal force $F_x$, a shear force $F_y$, and a bending moment $M_z$:

$$F_x = \int_S \sigma_x dS, \quad F_y = \int_S \tau_{xy} dS, \quad M_z = -\int_S y \sigma_x dS \qquad (3)$$

$\sigma_x$ and $\tau_{xy}$ are the normal and shear stresses of the beam. The shear stresses $\tau_{xy}$ are not considered in the Euler-Bernoulli theory. The virtual power principle is presented as follows:

$$\int_x \varepsilon^{*T} \mathbf{F} dx = \mathbf{u}^{*T} \mathbf{F}_{\text{ext}} \qquad (4)$$

$\varepsilon^{*T}$ et $\mathbf{u}^{*T}$ are two strain and displacement virtual fields, respectively. $\mathbf{F}_{\text{ext}}$ is the assembly of the external beam volumetric and boundary forces. The strain field $\varepsilon$ derives from the displacement field, so $\varepsilon(x) = \mathbf{B}u$. Eq. (12) is then simplified and expressed as follows:

$$\int_x \mathbf{B}^T \mathbf{F} dx = \mathbf{F}_{\text{ext}} \qquad (5)$$

where $\mathbf{F}_{\text{int}} = \int_x \mathbf{B}^{*T} \mathbf{F} dx$ is the internal forces vector of a beam element. The tangent operator of a beam element is deduced by deriving its internal forces vector $\mathbf{F}_{\text{int}}$ with respect to its degrees of freedom vector $\mathbf{u}$:

$$\frac{\partial \mathbf{F}_{\text{int}}}{\partial \mathbf{u}} = \int_x \mathbf{B}^T \mathbf{k} \mathbf{B} \, dx \qquad (6)$$

where $\mathbf{k}$ is the derivative of $\mathbf{F}_{\text{int}}$ with respect to $\varepsilon$. The literature offers more developed beam theories such as fiber beams where a discretization is applied to the beam section [33–39].

## 2.2. Plate elements

Kinematic assumptions relate displacements of point M of the plate to point G, the projection of M on the plate’s reference plane $\Omega$ (see Fig. 2):

$$\mathbf{u}(\mathbf{x}, \mathbf{z}_M) = \begin{bmatrix} \mathbf{u}_x(\mathbf{x}, \mathbf{z}_M) \\ \mathbf{v}_y(\mathbf{x}, \mathbf{z}_M) \\ \mathbf{w}_z(\mathbf{x}, \mathbf{z}_M) \end{bmatrix} = \begin{bmatrix} \mathbf{u}_x(\mathbf{x}) \\ \mathbf{v}_y(\mathbf{x}) \\ \mathbf{w}_z(\mathbf{x}) \end{bmatrix} + \mathbf{z}_M \begin{bmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \theta_x(\mathbf{x}) \\ \theta_y(\mathbf{x}) \\ \theta_z(\mathbf{x}) \end{bmatrix} \qquad (7)$$

$\theta_x$, $\theta_y$, and $\theta_z$ represent the rotation of the plate section with respect to x, y, and z axes. x is the position vector of G in the plate mid-plane.

### 2.2.1. Weak formulation

Solving the plate boundary value problem lies in finding the generalized displacement vector $\mathbf{u} = [u_x(\mathbf{x}) \ v_y(\mathbf{x}) \ w_z(\mathbf{x}) \ \theta_x(\mathbf{x}) \ \theta_y(\mathbf{x}) \ \theta_z(\mathbf{x})]^T$ at each position of the plate mid-plane ($z = 0$).

A generalized plate strain vector is defined as $\varepsilon = [e \ k \ \gamma]^T$, where:

$$\mathbf{e} = \begin{bmatrix} e_x = \frac{\partial u_x}{\partial x} \\ e_y = \frac{\partial v_y}{\partial y} \\ 2 e_{xy} = \frac{\partial u_x}{\partial y} + \frac{\partial v_y}{\partial x} \end{bmatrix}, \quad \mathbf{k} = \begin{bmatrix} \kappa_x = \frac{\partial \theta_y}{\partial x} \\ \kappa_y = - \frac{\partial \theta_x}{\partial y} \\ 2 \kappa_{xy} = \frac{\partial \theta_y}{\partial y} - \frac{\partial \theta_x}{\partial x} \end{bmatrix}, \quad \mathbf{\gamma} = \begin{bmatrix} \gamma_x = \frac{\partial w_z}{\partial x} + \theta_y \\ \gamma_y = \frac{\partial w_z}{\partial y} - \theta_x \end{bmatrix} \qquad (8)$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Plate configuration.</figcaption></figure>

knowing that the shear strain $\beta_y$ is considered within the Timoshenko beam theory and neglected when following the Euler-Bernoulli theory. A generalized forces and moments vector $\mathbf{F} = [F_x \ F_y \ M_z]^T$ is defined for each beam section. This section is subjected to three types of loading: a normal force $F_x$, a shear force $F_y$, and a bending moment $M_z$:

$$F_x = \int_S \sigma_x dS, \quad F_y = \int_S \tau_{xy} dS, \quad M_z = - \int_S y \sigma_x dS \qquad (3)$$

$\sigma_x$ and $\tau_{xy}$ are the normal and shear stresses of the beam. The shear stresses $\tau_{xy}$ are not considered in the Euler-Bernoulli theory. The virtual power principle is presented as follows:

$$\int_x \varepsilon^{*T} \mathbf{F} dx = \mathbf{u}^{*T} \mathbf{F}_{\text{ext}} \qquad (4)$$

$\varepsilon^{*T}$ et $\mathbf{u}^{*T}$ are two strain and displacement virtual fields, respectively. $\mathbf{F}_{\text{ext}}$ is the assembly of the external beam volumetric and boundary forces. The strain field $\varepsilon$ derives from the displacement field, so $\varepsilon(x) = \mathbf{B}u$. Eq. (12) is then simplified and expressed as follows:

$$\int_x \mathbf{B}^T \mathbf{F} dx = \mathbf{F}_{\text{ext}} \qquad (5)$$

where $\mathbf{F}_{\text{int}} = \int_x \mathbf{B}^{*T} \mathbf{F} dx$ is the internal forces vector of a beam element. The tangent operator of a beam element is deduced by deriving its internal forces vector $\mathbf{F}_{\text{int}}$ with respect to its degrees of freedom vector $\mathbf{u}$:

$$\frac{\partial \mathbf{F}_{\text{int}}}{\partial \mathbf{u}} = \int_x \mathbf{B}^T \mathbf{k} \mathbf{B} \, dx \qquad (6)$$

where $\mathbf{k}$ is the derivative of $\mathbf{F}_{\text{int}}$ with respect to $\varepsilon$. The literature offers more developed beam theories such as fiber beams where a discretization is applied to the beam section [33–39].

## 2.2. Plate elements

Kinematic assumptions relate displacements of point M of the plate to point G, the projection of M on the plate’s reference plane $\Omega$ (see Fig. 2):

$$\mathbf{u}(\mathbf{x}, \mathbf{z}_M) = \begin{bmatrix} u_x(\mathbf{x}, \mathbf{z}_M) \\ v_y(\mathbf{x}, \mathbf{z}_M) \\ w_z(\mathbf{x}, \mathbf{z}_M) \end{bmatrix} = \begin{bmatrix} u_x(\mathbf{x}) \\ v_y(\mathbf{x}) \\ w_z(\mathbf{x}) \end{bmatrix} + \mathbf{z}_M \begin{bmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \theta_x(\mathbf{x}) \\ \theta_y(\mathbf{x}) \\ \theta_z(\mathbf{x}) \end{bmatrix} \qquad (7)$$

$\theta_x$, $\theta_y$, and $\theta_z$ represent the rotation of the plate section with respect to x, y, and z axes. x is the position vector of G in the plate mid-plane.

### 2.2.1. Weak formulation

Solving the plate boundary value problem lies in finding the generalized displacement vector $\mathbf{u} = [u_x(\mathbf{x}) \ v_y(\mathbf{x}) \ w_z(\mathbf{x}) \ \theta_x(\mathbf{x}) \ \theta_y(\mathbf{x}) \ \theta_z(\mathbf{x})]^T$ at each position of the plate mid-plane ($z = 0$).

A generalized plate strain vector is defined as $\varepsilon = [e \ k \ \gamma]^T$, where:

$$\mathbf{e} = \begin{bmatrix} e_x = \frac{\partial u_x}{\partial x} \\ e_y = \frac{\partial v_y}{\partial y} \\ 2 e_{xy} = \frac{\partial u_x}{\partial y} + \frac{\partial v_y}{\partial x} \end{bmatrix}, \quad \mathbf{k} = \begin{bmatrix} \kappa_x = \frac{\partial \theta_y}{\partial x} \\ \kappa_y = - \frac{\partial \theta_x}{\partial y} \\ 2 \kappa_{xy} = \frac{\partial \theta_y}{\partial y} - \frac{\partial \theta_x}{\partial x} \end{bmatrix}, \quad \mathbf{\gamma} = \begin{bmatrix} \gamma_x = \frac{\partial w_z}{\partial x} + \theta_y \\ \gamma_y = \frac{\partial w_z}{\partial y} - \theta_x \end{bmatrix} \qquad (8)$$
$\epsilon$ represents the membrane strains vector, and $k$ is the vector of bending strains. $\gamma$ components are shear strains, neglected by the Kirchhoff-Love plate theory and taken into account in the Mindlin-Reissner theory. The stress vector $\sigma$ is composed of two components: $\sigma_m$ representing membrane stresses, and $\sigma_f$ representing bending stresses:

$$ \sigma = \left[ \begin{matrix} \sigma_{xx} & \sigma_{yy} & \sigma_{xy} \\ \hline \sigma_{m}^{T} & \sigma_{xz} & \sigma_{yz} \end{matrix} \right]^{T} \qquad (9) $$

The generalized forces vector $\mathbf{F}$ of the plate is defined as $\mathbf{F} = |\mathbf{N} \ \mathbf{M} \ \mathbf{T}|^T$, where:

$$ \mathbf{N} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \sigma_m dz, \quad \mathbf{M} = \int_{-\frac{h}{2}}^{\frac{h}{2}} z \sigma_m dz, \quad \mathbf{T} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \sigma_f dz \qquad (10) $$

The material behavior law of the plate links the stress $\sigma$ to the strain $\epsilon$. Subsequently, the relationship between $\mathbf{F}$ and $\epsilon$ is deduced. The virtual power principle is presented as follows:

$$ \int_{\Omega} \epsilon^{*T} \mathbf{F} d\Omega = \mathbf{u}^{*T} \mathbf{F}_{\text{ext}} \qquad (11) $$

$\epsilon^{*T}$ et $\mathbf{u}^{*T}$ are two strain and displacement virtual fields, respectively. $\mathbf{F}_{\text{ext}}$ is the assembly of the external plate mid-plane surfacic and boundary forces. The strain field $\epsilon$ derives from the displacement field, so $\epsilon(x) = \mathbf{B}u$. Eq. (11) is then simplified and expressed as follows:

$$ \int_{\Omega} \mathbf{B}^{T} \mathbf{F} d\Omega = \mathbf{F}_{\text{ext}} \qquad (12) $$

where $\mathbf{F}_{\text{int}} = \int_{\Omega} \mathbf{B}^{T} \mathbf{F} d\Omega$ is the internal forces vector of the plate element. The tangent operator of a beam element is deduced by deriving its internal forces vector $\mathbf{F}_{\text{int}}$ with respect to its degrees of freedom vector $\mathbf{u}$:

$$ \frac{\partial \mathbf{F}_{\text{int}}}{\partial \mathbf{u}} = \int_{\Omega} \mathbf{B}^{T} \mathbf{k} \mathbf{B} \, d\Omega \qquad (13) $$

where $\mathbf{k}$ is the derivative of $\mathbf{F}_{\text{int}}$ with respect to $\epsilon$. The literature offers more developed multilayered plate models including perfect or imperfect interfaces between the different layers [40–42].

## 2.3. Kinematic enhancement technique

Beam and plate elements presented in Sections 2.1 and 2.2 can be enhanced in order to model reinforced concrete structural elements, where the steel-concrete interface behavior is taken into consideration.

### 2.3.1. Principle

The proposed enhancement technique defines an enhanced element that assembles the following contributions:

* one or several beam or plate concrete finite elements;

* one or several truss or beam steel elements;

* steel-concrete bond stresses.

The assembly of the steel and the concrete elements creates one enhanced finite element, where inner bond stresses are considered between steel and concrete. These stresses contribute to the internal forces vector of the enhanced element.

The main advantage of the methodology is that existing truss, beam, and plate elements are called within an assembly of one enhanced element. The used beam or plate theory is a user choice. In addition, the composition of the enhanced element (the number of assembled finite elements and the type of each element) is a user choice. For the compositions of the enhanced elements defining inner nodes, a sub-structured resolution involving a static condensation technique is applied since these inner nodes are not linked to other finite elements of the studied structure. This technique aims to condensate the enhanced element rigidity matrix at the level of the outer nodes of this element [43]. Fig. 3 shows several enhanced concrete beams configurations with and without the need for static condensation. Several concrete plate enhanced elements configurations are shown in Fig. 4.

### 2.3.2. Weak formulation

For an enhanced element with $N_s$ steel elements, and $N_c$ concrete elements, the virtual power principle is expressed as follows:

$$ \mathbf{A}_{e=1}^{N_c} \mathbf{P}_{\text{intc}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{ints}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{intg}}^{e*} = \mathbf{A}_{e=1}^{N_c} \mathbf{P}_{\text{extc}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{exts}}^{e*} \qquad (14) $$

where:

* $\mathbf{P}_{\text{intc}}^{e*}$: is the virtual power of the internal forces of concrete elements (beam or plate elements);

* $\mathbf{P}_{\text{ints}}^{e*}$: is the virtual power of the internal forces of steel elements (truss or beam elements);

* $\mathbf{P}_{\text{intg}}^{e*}$: is the virtual power of the internal forces due to the interface bond stresses around the steel elements;
$\epsilon$ represents the membrane strains vector, and $k$ is the vector of bending strains. $\gamma$ components are shear strains, neglected by the Kirchhoff-Love plate theory and taken into account in the Mindlin-Reissner theory. The stress vector $\sigma$ is composed of two components: $\sigma_m$ representing membrane stresses, and $\sigma_f$ representing bending stresses:

$$ \sigma = \left[ \begin{matrix} \sigma_{xx} & \sigma_{yy} & \sigma_{xy} \\ \hline \sigma_{m}^{T} & \sigma_{xz} & \sigma_{yz} \end{matrix} \right]^{T} \qquad (9) $$

The generalized forces vector $\mathbf{F}$ of the plate is defined as $\mathbf{F} = |\mathbf{N} \ \mathbf{M} \ \mathbf{T}|^T$, where:

$$ \mathbf{N} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \sigma_m dz, \quad \mathbf{M} = \int_{-\frac{h}{2}}^{\frac{h}{2}} z \sigma_m dz, \quad \mathbf{T} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \sigma_f dz \qquad (10) $$

The material behavior law of the plate links the stress $\sigma$ to the strain $\epsilon$. Subsequently, the relationship between $\mathbf{F}$ and $\epsilon$ is deduced. The virtual power principle is presented as follows:

$$ \int_{\Omega} \epsilon^{*T} \mathbf{F} d\Omega = \mathbf{u}^{*T} \mathbf{F}_{\text{ext}} \qquad (11) $$

$\epsilon^{*T}$ et $\mathbf{u}^{*T}$ are two strain and displacement virtual fields, respectively. $\mathbf{F}_{\text{ext}}$ is the assembly of the external plate mid-plane surfacic and boundary forces. The strain field $\epsilon$ derives from the displacement field, so $\epsilon(\mathbf{x}) = \mathbf{B}\mathbf{u}$. Eq. (11) is then simplified and expressed as follows:

$$ \int_{\Omega} \mathbf{B}^{T} \mathbf{F} d\Omega = \mathbf{F}_{\text{ext}} \qquad (12) $$

where $\mathbf{F}_{\text{int}} = \int_{\Omega} \mathbf{B}^{T} \mathbf{F} d\Omega$ is the internal forces vector of the plate element. The tangent operator of a beam element is deduced by deriving its internal forces vector $\mathbf{F}_{\text{int}}$ with respect to its degrees of freedom vector $\mathbf{u}$:

$$ \frac{\partial \mathbf{F}_{\text{int}}}{\partial \mathbf{u}} = \int_{\Omega} \mathbf{B}^{T} \mathbf{k} \mathbf{B} \, d\Omega \qquad (13) $$

where $\mathbf{k}$ is the derivative of $\mathbf{F}_{\text{int}}$ with respect to $\epsilon$. The literature offers more developed multilayered plate models including perfect or imperfect interfaces between the different layers [40–42].

## 2.3. Kinematic enhancement technique

Beam and plate elements presented in Sections 2.1 and 2.2 can be enhanced in order to model reinforced concrete structural elements, where the steel-concrete interface behavior is taken into consideration.

### 2.3.1. Principle

The proposed enhancement technique defines an enhanced element that assembles the following contributions:

* one or several beam or plate concrete finite elements;

* one or several truss or beam steel elements;

* steel-concrete bond stresses.

The assembly of the steel and the concrete elements creates one enhanced finite element, where inner bond stresses are considered between steel and concrete. These stresses contribute to the internal forces vector of the enhanced element.

The main advantage of the methodology is that existing truss, beam, and plate elements are called within an assembly of one enhanced element. The used beam or plate theory is a user choice. In addition, the composition of the enhanced element (the number of assembled finite elements and the type of each element) is a user choice. For the compositions of the enhanced elements defining inner nodes, a sub-structured resolution involving a static condensation technique is applied since these inner nodes are not linked to other finite elements of the studied structure. This technique aims to condensate the enhanced element rigidity matrix at the level of the outer nodes of this element [43]. Fig. 3 shows several enhanced concrete beams configurations with and without the need for static condensation. Several concrete plate enhanced elements configurations are shown in Fig. 4.

### 2.3.2. Weak formulation

For an enhanced element with $N_s$ steel elements, and $N_c$ concrete elements, the virtual power principle is expressed as follows:

$$ \mathbf{A}_{e=1}^{N_c} \mathbf{P}_{\text{intc}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{ints}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{intg}}^{e*} = \mathbf{A}_{e=1}^{N_c} \mathbf{P}_{\text{extc}}^{e*} + \mathbf{A}_{e=1}^{N_s} \mathbf{P}_{\text{exts}}^{e*} \qquad (14) $$

where:

* $\mathbf{P}_{\text{intc}}^{e*}$: is the virtual power of the internal forces of concrete elements (beam or plate elements);

* $\mathbf{P}_{\text{ints}}^{e*}$: is the virtual power of the internal forces of steel elements (truss or beam elements);

* $\mathbf{P}_{\text{intg}}^{e*}$: is the virtual power of the internal forces due to the interface bond stresses around the steel elements;
<figure><img src="image_6.png" /><figcaption>Fig. 3. Enhanced concrete beams configurations with (a) and without (b) static condensation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Enhanced concrete plates configurations with (a) and without (b) static condensation.</figcaption></figure>

* $P_{extc}^e$: is the virtual power of the external forces applied to the concrete elements;

* $P_{exts}^e$: is the virtual power of the external forces applied to the steel elements.

It is possible to define an elementary degrees of freedom vector $u_{el}$ for the enhanced element. $u_{el}$ is the vector of the steel and the concrete degrees of freedom. Eq. (14) is then written in a more condensed way as follows:

$$u_{el}^* T F_i^e(u_{el}) = u_{el}^* T F_e \qquad (15)$$

$u_{el}^*$ is a virtual elementary degrees of freedom vector. $F_i^e(u_{el})$ is the vector that combines all internal forces (steel, concrete, and internal forces due to steel-concrete interface behavior) of the enhanced element. $F_e$ represents the vector of external forces on the enhanced element.

For enhanced elements assembly choices involving internal nodes and thus requiring a sub-structured resolution and a static condensation, $F_i^e$ is condensed at the level of the enhanced element external nodes. These nodes are linked to the other finite elements of the adopted structural mesh. Similarly, the stiffness matrix of the enhanced element is condensed at the level of its outer nodes. Appendix A details the condensation technique formulation and the associated resolution.

The following subsections detail the formulation of three enhanced elements structures: two configurations of enhanced concrete beams and one configuration of an enhanced concrete plate.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Enhanced concrete beams configurations with (a) and without (b) static condensation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Enhanced concrete plates configurations with (a) and without (b) static condensation.</figcaption></figure>

* $P_{extc}^e$: is the virtual power of the external forces applied to the concrete elements;

* $P_{exts}^e$: is the virtual power of the external forces applied to the steel elements.

It is possible to define an elementary degrees of freedom vector $u_{el}$ for the enhanced element. $u_{el}$ is the vector of the steel and the concrete degrees of freedom. Eq. (14) is then written in a more condensed way as follows:

$$u_{el}^* T F_i^e(u_{el}) = u_{el}^* T F_e \qquad (15)$$

$u_{el}^*$ is a virtual elementary degrees of freedom vector. $F_i^e(u_{el})$ is the vector that combines all internal forces (steel, concrete, and internal forces due to steel-concrete interface behavior) of the enhanced element. $F_e$ represents the vector of external forces on the enhanced element.

For enhanced elements assembly choices involving internal nodes and thus requiring a sub-structured resolution and a static condensation, $F_i^e$ is condensed at the level of the enhanced element external nodes. These nodes are linked to the other finite elements of the adopted structural mesh. Similarly, the stiffness matrix of the enhanced element is condensed at the level of its outer nodes. Appendix A details the condensation technique formulation and the associated resolution.

The following subsections detail the formulation of three enhanced elements structures: two configurations of enhanced concrete beams and one configuration of an enhanced concrete plate.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Enhanced concrete beam configuration (a) and degrees of freedom (b).</figcaption></figure>

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

## 2.4. Enhanced element of a concrete beam and a steel truss

Fig. 5(a) describes the enhanced beam element composition with a concrete beam element, a truss steel element, and bond stresses. Two nodes are defined, with four degrees of freedom each (three concrete degrees of freedom and one steel degree of freedom), as shown in Fig. 5(b).

A steel truss element is selected, assigning one longitudinal degree of freedom per node ($u_{sl}$ for node i) for steel's longitudinal displacement. Alternatively, a steel beam element could be employed, which would involve defining three steel degrees of freedom per node: two displacements and one rotation.

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as $u_{el} = [u_1 \ v_1 \ \theta_1 \ u_{s1} \ u_2 \ v_2 \ \theta_2 \ u_{s2}]^T$. It is possible to identify a concrete degrees of freedom vector $u_c$ and a steel degrees of freedom vector $u_s$ such that:

$$u_c = [u_1 \ v_1 \ \theta_1 \ u_2 \ v_2 \ \theta_2]^T, \ u_s = [u_{s1} \ u_{s2}]^T \qquad (16)$$

where:

$$u_c = \begin{bmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 \end{bmatrix} u_{el} = I_c u_{el} \qquad (17)$$

and:

$$u_s = [0 \ 0 \ 0 \ 1 \ 0 \ 0 \ 0 \ 1] u_{el} = I_s u_{el} \qquad (18)$$

The degrees of freedom at a longitudinal position x are interpolated using the matrix of shape functions N as follows:

$$\begin{bmatrix} u(x) \\ v(x) \\ \theta(x) \\ u_s(x) \end{bmatrix} = N u_{el} = \begin{bmatrix} N_b u_c \\ N_t u_s \end{bmatrix} \qquad (19)$$

so:

$$\begin{bmatrix} u(x) \\ v(x) \\ \theta(x) \\ u_s(x) \end{bmatrix} = \begin{bmatrix} N_b I_c \\ N_t I_s \end{bmatrix} u_{el} = \begin{bmatrix} N_c \\ N_s \end{bmatrix} u_{el} \qquad (20)$$

$N_b$ is the shape functions matrix of the used concrete beam element which can be a user choice (for example generalized or fiber Timoshenko or Euler-Bernoulli beam). $N_t$ is the shape functions matrix of the steel truss element. Similarly:

$$\epsilon_c = \begin{bmatrix} \frac{du}{dx} \\ \beta_y \\ \frac{\partial\theta}{\partial x} \end{bmatrix} = B_b u_c \qquad (21)$$

and:

$$\epsilon_s = \frac{\partial u_s(x)}{\partial x} = B_t u_s \qquad (22)$$

$\epsilon_c$ is the generalized strain vector of the concrete beam (see Eq. (2)). $\epsilon_s$ is the steel longitudinal strain. b and t indices of $B_b$ and $B_t$ matrices stand for beam and truss. So:

$$\epsilon_c = B_b I_c u_{el} = B_c u_{el} \qquad (23)$$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Enhanced concrete beam configuration (a) and degrees of freedom (b).</figcaption></figure>

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

## 2.4. Enhanced element of a concrete beam and a steel truss

Fig. 5(a) describes the enhanced beam element composition with a concrete beam element, a truss steel element, and bond stresses. Two nodes are defined, with four degrees of freedom each (three concrete degrees of freedom and one steel degree of freedom), as shown in Fig. 5(b).

A steel truss element is selected, assigning one longitudinal degree of freedom per node ($u_{sl}$ for node i) for steel's longitudinal displacement. Alternatively, a steel beam element could be employed, which would involve defining three steel degrees of freedom per node: two displacements and one rotation.

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as $u_{el} = [u_1 \ v_1 \ \theta_1 \ u_{s1} \ u_2 \ v_2 \ \theta_2 \ u_{s2}]^T$. It is possible to identify a concrete degrees of freedom vector $u_c$ and a steel degrees of freedom vector $u_s$ such that:

$$u_c = [u_1 \ v_1 \ \theta_1 \ u_2 \ v_2 \ \theta_2]^T, \ u_s = [u_{s1} \ u_{s2}]^T \qquad (16)$$

where:

$$u_c = \begin{bmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 \end{bmatrix} u_{el} = I_c u_{el} \qquad (17)$$

and:

$$u_s = [0 \ 0 \ 0 \ 1 \ 0 \ 0 \ 0 \ 1] u_{el} = I_s u_{el} \qquad (18)$$

The degrees of freedom at a longitudinal position x are interpolated using the matrix of shape functions N as follows:

$$\begin{bmatrix} u(x) \\ v(x) \\ \theta(x) \\ u_s(x) \end{bmatrix} = N u_{el} = \begin{bmatrix} N_b u_c \\ N_t u_s \end{bmatrix} \qquad (19)$$

so:

$$\begin{bmatrix} u(x) \\ v(x) \\ \theta(x) \\ u_s(x) \end{bmatrix} = \begin{bmatrix} N_b I_c \\ N_t I_s \end{bmatrix} u_{el} = \begin{bmatrix} N_c \\ N_s \end{bmatrix} u_{el} \qquad (20)$$

$N_b$ is the shape functions matrix of the used concrete beam element which can be a user choice (for example generalized or fiber Timoshenko or Euler-Bernoulli beam). $N_t$ is the shape functions matrix of the steel truss element. Similarly:

$$\epsilon_c = \begin{bmatrix} \frac{du}{dx} \\ \beta_y \\ \frac{\partial\theta}{\partial x} \end{bmatrix} = B_b u_c \qquad (21)$$

and:

$$\epsilon_s = \frac{\partial u_s(x)}{\partial x} = B_t u_s \qquad (22)$$

$\epsilon_c$ is the generalized strain vector of the concrete beam (see Eq. (2)). $\epsilon_s$ is the steel longitudinal strain. b and t indices of $B_b$ and $B_t$ matrices stand for beam and truss. So:

$$\epsilon_c = B_b I_c u_{el} = B_c u_{el} \qquad (23)$$
and:

$$ \epsilon_s = \mathbf{B}_t \mathbf{I}_s \mathbf{u}_{el} = \mathbf{B}_s \mathbf{u}_{el} \qquad (24) $$

### 2.4.1. Steel-concrete slip evaluation

The slip $u_g(x)$ is computed as follows:

$$ u_g(x) = u_s(x) - u_c(x) \qquad (25) $$

Considering that $u_c(x)$ represents the concrete displacement near the steel bar, calculated according to the beam kinematics (see Eq. (1) and Fig. 1). Therefore:

$$ u_c(x) = \underbrace{[1 \quad -y]}_{a_s} \mathbf{N}_{13} \mathbf{u}_{el} \qquad (26) $$

where y is the eccentricity of the steel bar with respect to the concrete beam neutral axis. $\mathbf{N}_{13}$ is the matrix that holds the first and third rows of the shape functions matrix $\mathbf{N}$ (see Eq. (19)). So:

$$ u_g(x) = \mathbf{N}_s \mathbf{u}_{el} - \mathbf{a}_s \mathbf{N}_{13} \mathbf{u}_{el} = (\mathbf{N}_s - \mathbf{a}_s \mathbf{N}_{13}) \mathbf{u}_{el} \qquad (27) $$

### 2.4.2. Weak formulation

The virtual power principle takes into account the behavior of concrete, steel, and bond stresses within the enhanced element framework:

$$ \mathbf{u}_{el}^{*T} \mathbf{F}_{intc}^{e} + \int_{0}^{L_s} \varepsilon_s^* \sigma_s(\varepsilon_s) S_s dx + \int_{0}^{L_s} u_g^* \tau(u_g) P dx = \mathbf{u}_{el}^{*T} \mathbf{F}_e \qquad (28) $$

$\mathbf{F}_{intc}^e$ represents the internal forces vector of the concrete beam (see Eq. (13)). The length $L_s$ of the steel bar is equal to the length $l$ of the beam. The steel constitutive law links the steel stress $\sigma_s$ to its strain $\varepsilon_s$. $u_g^*$ is a virtual steel-concrete slip, and $P$ is the steel bar perimeter. The virtual power principle of the enhanced element as follows:

$$ \mathbf{u}_{el}^{*T} \mathbf{F}_{intc}^{e} + \mathbf{u}_{el}^{*T} \int_{0}^{L_s} \mathbf{B}_{s}^{T} \sigma_{s}(\varepsilon_{s}) S_{s} dx \mathbf{u}_{el} + \mathbf{u}_{el}^{*T} \int_{0}^{L_s} (\mathbf{N}_{s}^{T} - \mathbf{N}_{13}^{T} \mathbf{a}_{s}^{T}) \tau(u_{g}) P dx = \mathbf{u}_{el}^{*T} \mathbf{F}_{e} \qquad (29) $$

After simplifying Eq. (29) by $\mathbf{u}_{el}^{*T}$, the balance between internal and external forces of the enhanced element is described as follows:

$$ \mathbf{F}_{l}^{e} = \mathbf{F}_{e} \qquad (30) $$

The computation of the tangent operator of the enhanced element consists in calculating the derivative of its internal forces vector $\mathbf{F}_{l}^{e}$ with respect to the degrees of freedom vector $\mathbf{u}_{el}$:

$$ \frac{\partial \mathbf{F}_{l}^{e}}{\partial \mathbf{u}_{el}} = k_{sc}^{e} = k_{c}^{e} + k_{s}^{e} + k_{g}^{e} \qquad (31) $$

$k_c^e$ and $k_s^e$ represent the stiffness matrices of concrete and steel elements, respectively (classical beam and truss finite elements) and $k_g^e$ is the contribution of the interface bond behavior in the computation of the tangent operator $k_{sc}^e$ for the enhanced element:

$$ k_c^e = \int_0^l \mathbf{B}_c^T \mathbf{k}_c \mathbf{B}_c \, dx \qquad (32) $$

$$ k_s^e = \int_{0}^{L_s} \mathbf{B}_s^T \mathbf{k}_s \mathbf{B}_s \, dx \qquad (33) $$

$$ k_g^e = \int_{0}^{L_s} (\mathbf{N}_s^T - \mathbf{N}_{13}^T \mathbf{a}_s^T) \frac{\partial \tau(u_g)}{\partial \mathbf{u}_{el}} P dx \qquad (34) $$

with:

$$ \frac{\partial \tau(u_g)}{\partial \mathbf{u}_{el}} = \underbrace{\frac{\partial \tau(u_g)}{\partial u_g}}_{k_{int}} \frac{\partial u_g}{\partial \mathbf{u}_{el}} = k_{int}(\mathbf{N}_s - \mathbf{a}_s \mathbf{N}_{13}) \qquad (35) $$

$\mathbf{k}_c$ and $\mathbf{k}_s$ are concrete and steel constitutive laws matrices, respectively. $k_{int}$ is calculated due to the expression of a steel-interface constitutive bond law which links the bond stress $\tau$ to the steel-concrete slip $u_g$.

Pull-out tests are usually adopted to identify and characterize the bond laws. Several analytical expressions for constitutive bond laws are proposed in the literature [44–54]. The bond law of [54] is chosen for the numerical applications of this paper for its simplicity and its representativity for monotonic and cyclic bond behaviors, as shown in Fig. 6. Three parameters define this bond law: the maximum bond strength $\tau_1$, the slip $g_1$ for which $\tau_1$ is reached, and the slip $g_3$. For slip values bigger than $g_3$ the total stress of the monotonic version of the law remains constant.
and:

$$ \epsilon_s = \mathbf{B}_t \mathbf{I}_s \mathbf{u}_{el} = \mathbf{B}_s \mathbf{u}_{el} \qquad (24) $$

### 2.4.1. Steel-concrete slip evaluation

The slip $u_g(x)$ is computed as follows:

$$ u_g(x) = u_s(x) - u_c(x) \qquad (25) $$

Considering that $u_c(x)$ represents the concrete displacement near the steel bar, calculated according to the beam kinematics (see Eq. (1) and Fig. 1). Therefore:

$$ u_c(x) = \underbrace{[1 \quad -y]}_{a_s} \mathbf{N}_{13} \mathbf{u}_{el} \qquad (26) $$

where y is the eccentricity of the steel bar with respect to the concrete beam neutral axis. $\mathbf{N}_{13}$ is the matrix that holds the first and third rows of the shape functions matrix $\mathbf{N}$ (see Eq. (19)). So:

$$ u_g(x) = \mathbf{N}_s \mathbf{u}_{el} - \mathbf{a}_s \mathbf{N}_{13} \mathbf{u}_{el} = (\mathbf{N}_s - \mathbf{a}_s \mathbf{N}_{13}) \mathbf{u}_{el} \qquad (27) $$

### 2.4.2. Weak formulation

The virtual power principle takes into account the behavior of concrete, steel, and bond stresses within the enhanced element framework:

$$ \mathbf{u}_{el}^{*T} \mathbf{F}_{intc}^{e} + \int_{0}^{L_s} \varepsilon_s^* \sigma_s(\varepsilon_s) S_s dx + \int_{0}^{L_s} u_g^* \tau(u_g) P dx = \mathbf{u}_{el}^{*T} \mathbf{F}_e \qquad (28) $$

$\mathbf{F}_{intc}^e$ represents the internal forces vector of the concrete beam (see Eq. (13)). The length $L_s$ of the steel bar is equal to the length $l$ of the beam. The steel constitutive law links the steel stress $\sigma_s$ to its strain $\varepsilon_s$. $u_g^*$ is a virtual steel-concrete slip, and $P$ is the steel bar perimeter. The virtual power principle of the enhanced element as follows:

$$ \mathbf{u}_{el}^{*T} \mathbf{F}_{intc}^{e} + \mathbf{u}_{el}^{*T} \int_{0}^{L_s} \mathbf{B}_{s}^{T} \sigma_{s}(\varepsilon_{s}) S_{s} dx \mathbf{u}_{el} + \mathbf{u}_{el}^{*T} \int_{0}^{L_s} (\mathbf{N}_{s}^{T} - \mathbf{N}_{13}^{T} \mathbf{a}_{s}^{T}) \tau(u_{g}) P dx = \mathbf{u}_{el}^{*T} \mathbf{F}_{e} \qquad (29) $$

After simplifying Eq. (29) by $\mathbf{u}_{el}^{*T}$, the balance between internal and external forces of the enhanced element is described as follows:

$$ \mathbf{F}_{l}^{e} = \mathbf{F}_{e} \qquad (30) $$

The computation of the tangent operator of the enhanced element consists in calculating the derivative of its internal forces vector $\mathbf{F}_{l}^{e}$ with respect to the degrees of freedom vector $\mathbf{u}_{el}$:

$$ \frac{\partial \mathbf{F}_{l}^{e}}{\partial \mathbf{u}_{el}} = k_{sc}^{e} = k_{c}^{e} + k_{s}^{e} + k_{g}^{e} \qquad (31) $$

$k_c^e$ and $k_s^e$ represent the stiffness matrices of concrete and steel elements, respectively (classical beam and truss finite elements) and $k_g^e$ is the contribution of the interface bond behavior in the computation of the tangent operator $k_{sc}^e$ for the enhanced element:

$$ k_c^e = \int_0^l \mathbf{B}_c^T \mathbf{k}_c \mathbf{B}_c \, dx \qquad (32) $$

$$ k_s^e = \int_0^{L_s} \mathbf{B}_s^T \mathbf{k}_s \mathbf{B}_s \, dx \qquad (33) $$

$$ k_g^e = \int_0^{L_s} (\mathbf{N}_s^T - \mathbf{N}_{13}^T \mathbf{a}_s^T) \frac{\partial \tau(u_g)}{\partial \mathbf{u}_{el}} P dx \qquad (34) $$

with:

$$ \frac{\partial \tau(u_g)}{\partial \mathbf{u}_{el}} = \underbrace{\frac{\partial \tau(u_g)}{\partial u_g}}_{k_{int}} \frac{\partial u_g}{\partial \mathbf{u}_{el}} = k_{int}(\mathbf{N}_s - \mathbf{a}_s \mathbf{N}_{13}) \qquad (35) $$

$\mathbf{k}_c$ and $\mathbf{k}_s$ are concrete and steel constitutive laws matrices, respectively. $k_{int}$ is calculated due to the expression of a steel-interface constitutive bond law which links the bond stress $\tau$ to the steel-concrete slip $u_g$.

Pull-out tests are usually adopted to identify and characterize the bond laws. Several analytical expressions for constitutive bond laws are proposed in the literature [44–54]. The bond law of [54] is chosen for the numerical applications of this paper for its simplicity and its representativity for monotonic and cyclic bond behaviors, as shown in Fig. 6. Three parameters define this bond law: the maximum bond strength $\tau_1$, the slip $g_1$ for which $\tau_1$ is reached, and the slip $g_3$. For slip values bigger than $g_3$ the total stress of the monotonic version of the law remains constant.
<figure><img src="image_10.png" /><figcaption>Fig. 6. Bond law proposed in [54]: monotonic envelope curve (a); cyclic response (b).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Enhanced concrete beams configuration (a) and degrees of freedom (b).</figcaption></figure>

It is important to note here that the presented enhanced element is equivalent to the displacement beam element of [24] developed in a large displacement framework in [25]. However, the enhanced methodology of this work proposes a new way to construct the enhanced element, by calling existing concrete and steel elements from a finite elements code library. This approach facilitates the construction of the enhanced element and its implementation. In addition, it allows to assemble the concrete and steel elements finite elements of the user's choice. Moreover, considering large displacements can be done, in the framework of the proposed approach, by calling co-rotational beam elements. These elements carry on the large displacement aspect without the necessity of complicating the enhanced element formulation.

A sub-structured resolution is possible using the enhanced element method of this paper, which allows the transition from an enhanced element construction composed of a single concrete beam and a single steel truss to another one with two or multiple concrete beams or plates and steel bar or beam elements.

For applications where a higher interpolation degree is chosen for steel bars, the enhanced concrete beam element of the following subsection can be employed

## 2.5. Enhanced element of two concrete beams and a three-node steel bar

Fig. 7(a) describes the enhanced beam element composition with two concrete beam elements, a three-node steel bar element, and bond stresses. Three nodes are defined, with four degrees of freedom each (arranged as follows: three beam concrete degrees of freedom and one steel bar degree of freedom), as shown in Fig. 7(b). This enhanced element presents two external nodes and an internal one (see Fig. 7(b)). The formulation of this enhanced element is detailed in Appendix B.

## 2.6. Enhanced element of a concrete plate and two steel trusses

Fig. 8 describes the enhanced plate element composition with one concrete plate element, two truss steel elements, and steel bars-concrete bond stresses. Twelve nodes are defined: eight concrete nodes of six degrees of freedom each (three displacements and three rotations), and four steel nodes of one degree of freedom each (longitudinal displacement).

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as

$$ u_{el} = [u_1 \ v_1 \ w_1 \ \theta_{x1} \ \theta_{y1} \ \theta_{z1} \ u_2 \ v_2 \ w_2 \ \theta_{x2} \ \theta_{y2} \ \theta_{z2} \ \dots \ u_9 \ u_{10} \ u_{11} \ u_{12}]^T, \text{ such that:} $$

$$ u_{el} = [u_c \ u_{sx} \ u_{sy}]^T \qquad (36) $$

where:

$$ u_{sx} = [u_9 \ u_{10}]^T, \ u_{sy} = [u_{11} \ u_{12}]^T \qquad (37) $$
<figure><img src="image_10.png" /><figcaption>Fig. 6. Bond law proposed in [54]: monotonic envelope curve (a); cyclic response (b).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Enhanced concrete beams configuration (a) and degrees of freedom (b).</figcaption></figure>

It is important to note here that the presented enhanced element is equivalent to the displacement beam element of [24] developed in a large displacement framework in [25]. However, the enhanced methodology of this work proposes a new way to construct the enhanced element, by calling existing concrete and steel elements from a finite elements code library. This approach facilitates the construction of the enhanced element and its implementation. In addition, it allows to assemble the concrete and steel elements finite elements of the user's choice. Moreover, considering large displacements can be done, in the framework of the proposed approach, by calling co-rotational beam elements. These elements carry on the large displacement aspect without the necessity of complicating the enhanced element formulation.

A sub-structured resolution is possible using the enhanced element method of this paper, which allows the transition from an enhanced element construction composed of a single concrete beam and a single steel truss to another one with two or multiple concrete beams or plates and steel bar or beam elements.

For applications where a higher interpolation degree is chosen for steel bars, the enhanced concrete beam element of the following subsection can be employed

## 2.5. Enhanced element of two concrete beams and a three-node steel bar

Fig. 7(a) describes the enhanced beam element composition with two concrete beam elements, a three-node steel bar element, and bond stresses. Three nodes are defined, with four degrees of freedom each (arranged as follows: three beam concrete degrees of freedom and one steel bar degree of freedom), as shown in Fig. 7(b). This enhanced element presents two external nodes and an internal one (see Fig. 7(b)). The formulation of this enhanced element is detailed in Appendix B.

## 2.6. Enhanced element of a concrete plate and two steel trusses

Fig. 8 describes the enhanced plate element composition with one concrete plate element, two truss steel elements, and steel bars-concrete bond stresses. Twelve nodes are defined: eight concrete nodes of six degrees of freedom each (three displacements and three rotations), and four steel nodes of one degree of freedom each (longitudinal displacement).

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as

$$u_{el} = [u_1 \ v_1 \ w_1 \ \theta_{x1} \ \theta_{y1} \ \theta_{z1} \ u_2 \ v_2 \ w_2 \ \theta_{x2} \ \theta_{y2} \ \theta_{z2} \ \dots \ u_9 \ u_{10} \ u_{11} \ u_{12}]^T, \text{ such that:}$$

$$u_{el} = [u_c \ u_{sx} \ u_{sy}]^T \qquad (36)$$

where:

$$u_{sx} = [u_9 \ u_{10}]^T, u_{sy} = [u_{11} \ u_{12}]^T \qquad (37)$$
<figure><img src="image_12.png" /><figcaption>Fig. 8. Enhanced concrete plate configuration: 3d configuration (a) and top view (b).</figcaption></figure>

$u_c$ represent the concrete elementary degrees of freedom vector. $u_{sx}$ and $u_{sy}$ are the two steel bars degrees of freedom vectors.

The concrete degrees of freedom at a position x of the plate mid-plane are interpolated using the plate shape functions matrix $N_p$ as follows:

$$[u(x) \ v(x) \ w(x) \ \theta_x(x) \ \theta_y(x) \ \theta_z(x)]^T = N_p u_c \quad (38)$$

where:

$$u_c = \underbrace{[I_{d_{48 \times 48}} \quad 0_{48 \times 4}]}_{I_c} u_{el} \quad (39)$$

knowing that $I_{d_{i \times j}}$ is an $i \times i$ eye matrix and $0_{i \times j}$ is a $i \times j$ zero matrix. Hence:

$$[u(x) \ v(x) \ w(x) \ \theta_x(x) \ \theta_y(x) \ \theta_z(x)]^T = N_p I_c u_{el} = N_c u_{el} \quad (40)$$

The steel longitudinal displacement along the x steel bar is interpolated as follows:

$$u_{sx}(x) = N_{tx} u_{sx} \quad (41)$$

where:

$$u_{sx} = \underbrace{[0_{2 \times 48} \quad I_{d_{2 \times 2}} \quad 0_{d_{2 \times 2}}]}_{I_{sx}} u_{el} \quad (42)$$

Hence:

$$u_{sx}(x) = N_{tx} I_{sx} u_{el} = N_{sx} u_{el} \quad (43)$$

Similarly, the steel longitudinal displacement along the y steel bar is deduced as follows:

$$u_{sy}(y) = N_{ty} I_{sy} u_{el} = N_{sy} u_{el} \quad (44)$$

where:

$$I_{sy} = [0_{2 \times 50} \quad I_{d_{2 \times 2}}] \quad (45)$$

$N_{tx}$ and $N_{ty}$ are the shape functions of the x and y steel truss elements, respectively. The concrete and steel strain fields are derived from the corresponding generalized displacements fields, so:

$$\epsilon_c = B_p u_c = B_c u_{el} \quad \text{with} \quad B_c = B_p I_c \quad (46)$$

$$\epsilon_{sx} = B_{tx} u_{sx} = B_{sx} u_{el} \quad \text{with} \quad B_{sx} = B_{tx} I_{sx} \quad (47)$$

$$\epsilon_{sy} = B_{ty} u_{sy} = B_{sy} u_{el} \quad \text{with} \quad B_{sy} = B_{ty} I_{sy} \quad (48)$$

$\epsilon_{sx}$ and $\epsilon_{sy}$ are the longitudinal strains of the steel bars. $\epsilon_c$ is the generalized concrete plate strain vector (see Eq. (8) that defines the terms of $\epsilon_c$).
<figure><img src="image_12.png" /><figcaption>Fig. 8. Enhanced concrete plate configuration: 3d configuration (a) and top view (b).</figcaption></figure>

$u_c$ represent the concrete elementary degrees of freedom vector. $u_{sx}$ and $u_{sy}$ are the two steel bars degrees of freedom vectors.

The concrete degrees of freedom at a position x of the plate mid-plane are interpolated using the plate shape functions matrix $N_p$ as follows:

$$[u(x) \ v(x) \ w(x) \ \theta_x(x) \ \theta_y(x) \ \theta_z(x)]^T = N_p u_c \quad (38)$$

where:

$$u_c = \underbrace{\begin{bmatrix} I_{d_{48 \times 48}} & 0_{48 \times 4} \end{bmatrix}}_{I_c} u_{el} \quad (39)$$

knowing that $I_{d_{i \times j}}$ is an $i \times i$ eye matrix and $0_{i \times j}$ is a $i \times j$ zero matrix. Hence:

$$[u(x) \ v(x) \ w(x) \ \theta_x(x) \ \theta_y(x) \ \theta_z(x)]^T = N_p I_c u_{el} = N_c u_{el} \quad (40)$$

The steel longitudinal displacement along the x steel bar is interpolated as follows:

$$u_{sx}(x) = N_{tx} u_{sx} \quad (41)$$

where:

$$u_{sx} = \underbrace{\begin{bmatrix} 0_{2 \times 48} & I_{d_{2 \times 2}} & 0_{d_{2 \times 2}} \end{bmatrix}}_{I_{sx}} u_{el} \quad (42)$$

Hence:

$$u_{sx}(x) = N_{tx} I_{sx} u_{el} = N_{sx} u_{el} \quad (43)$$

Similarly, the steel longitudinal displacement along the y steel bar is deduced as follows:

$$u_{sy}(y) = N_{ty} I_{sy} u_{el} = N_{sy} u_{el} \quad (44)$$

where:

$$I_{sy} = \begin{bmatrix} 0_{2 \times 50} & I_{d_{2 \times 2}} \end{bmatrix} \quad (45)$$

$N_{tx}$ and $N_{ty}$ are the shape functions of the x and y steel truss elements, respectively. The concrete and steel strain fields are derived from the corresponding generalized displacements fields, so:

$$\epsilon_c = B_p u_c = B_c u_{el} \quad \text{with} \quad B_c = B_p I_c \quad (46)$$

$$\epsilon_{sx} = B_{tx} u_{sx} = B_{sx} u_{el} \quad \text{with} \quad B_{sx} = B_{tx} I_{sx} \quad (47)$$

$$\epsilon_{sy} = B_{ty} u_{sy} = B_{sy} u_{el} \quad \text{with} \quad B_{sy} = B_{ty} I_{sy} \quad (48)$$

$\epsilon_{sx}$ and $\epsilon_{sy}$ are the longitudinal strains of the steel bars. $\epsilon_c$ is the generalized concrete plate strain vector (see Eq. (8) that defines the terms of $\epsilon_c$).
## 2.6.1. Steel-concrete slip evaluation

The longitudinal slips of the two steel bars are independently calculated with respect to concrete. For the x bar:

$$u_{gx}(x) = u_{sx}(x) - u_{cx}(x) \qquad (49)$$

where:

$$u_{cx}(x) = \underbrace{[1 \ 0 \ 0 \ 0 \ z_{ax} \ 0]}_{a_x} N_p u_c = a_x N_c u_{el} \qquad (50)$$

which gives:

$$u_{gx}(x) = N_{sx}u_{el} - a_x N_c u_{el} = (N_{sx} - a_x N_c) u_{el} \qquad (51)$$

Similarly, the y steel bar longitudinal slip with respect to concrete is evaluated as follows:

$$u_{gy}(y) = u_{sy}(y) - u_{cy}(y) = (N_{sy} - a_y N_c) u_{el} \qquad (52)$$

with:

$$a_y = [0 \ 1 \ 0 \ -z_{ay} \ 0 \ 0] \qquad (53)$$

## 2.6.2. Weak formulation

The virtual power principle of the enhanced plate element is expressed as follows:

$$\begin{aligned} & \mathbf{u}_{el}^* \mathbf{F}_{intc}^e + \int_0^{L_x} \varepsilon_{sx}^* \sigma_{sx} (\varepsilon_{sx}) \varepsilon_{sx} dx + \int_0^{L_y} \varepsilon_{sy}^* \sigma_{sy} (\varepsilon_{sy}) \varepsilon_{sy} dy \\ & + \int_0^{L_x} u_{gx}^* \tau_x (u_{gx}) P_x dx + \int_0^{L_y} u_{gy}^* \tau_y (u_{gy}) P_y dy = \mathbf{u}_{el}^* \mathbf{F}_e \end{aligned} \qquad (54)$$

which gives after a simplification with respect to $\mathbf{u}_{el}^* \mathbf{F}_e$, the following expression:

$$\begin{aligned} & \mathbf{F}_{intc}^e + \int_0^{L_x} \mathbf{B}_{sx}^T \sigma_{sx} (\varepsilon_{sx}) \varepsilon_{sx} dx + \int_0^{L_y} \mathbf{B}_{sy}^T \sigma_{sy} (\varepsilon_{sy}) \varepsilon_{sy} dy \\ & + \int_0^{L_x} (N_{sx}^T - N_c^T a_x^T) \tau_x (u_{gx}) \mathbf{u}_{el} P_x dx + \int_0^{L_y} (N_{sy}^T - N_c^T a_y^T) \tau_y (u_{gy}) \mathbf{u}_{el} P_y dy = \mathbf{F}_e \end{aligned} \qquad (55)$$

where $\tau_x$ and $\tau_y$ are the bond stresses around the x and y bars of perimeters $P_x$ and $P_y$, respectively. The entirety of the expression on the left-hand side of Eq. (55) denotes for the internal forces vector $\mathbf{F}_e^i$ of the enhanced plate element. Its derivation with respect to $\mathbf{u}_{el}$ gives the tangent operator of the enhanced element:

$$\frac{\partial \mathbf{F}_e^e}{\partial \mathbf{u}_{el}} = k_{sc}^e = k_c^e + k_{sx}^e + k_{sy}^e + k_{gx}^e + k_{gy}^e \qquad (56)$$

$k_c^e$ is the tangent operator of the concrete plate element. $k_{sx}^e$ and $k_{sy}^e$ are the tangent operators of the two steel bars. $k_{gx}^e$ and $k_{gy}^e$ are the tangent operators due to the steel-concrete bond of the two steel bars, calculated as follows:

$$k_{gx}^e = \int_0^{L_x} (N_{sx}^T - N_c^T a_x^T) \frac{\partial \tau_x (u_{gx})}{\partial \mathbf{u}_{el}} P_x dx \qquad (57)$$

$$k_{gy}^e = \int_0^{L_y} (N_{sy}^T - N_c^T a_y^T) \frac{\partial \tau_y (u_{gy})}{\partial \mathbf{u}_{el}} P_y dy \qquad (58)$$

where:

$$\frac{\partial \tau_x (u_{gx})}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_x (u_{gx})}{\partial u_{gx}} \frac{\partial u_{gx}}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_x (u_{gx})}{\partial u_{gx}} (N_{sx} - a_x N_c) \qquad (59)$$

$$\frac{\partial \tau_y (u_{gy})}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_y (u_{gy})}{\partial u_{gy}} \frac{\partial u_{gy}}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_y (u_{gy})}{\partial u_{gy}} (N_{sy} - a_y N_c) \qquad (60)$$

knowing that $\frac{\partial \tau_x (u_{gx})}{\partial u_{gx}}$ and $\frac{\partial \tau_y (u_{gy})}{\partial u_{gy}}$ are the slopes of the bond laws considered for bond stresses around two steel bars. The formulation allows the utilization of two distinct laws.

# 3. Verification tests

## 3.1. Pull-out test

In this subsection, the pull-out experimental test of [55] is modeled.
## 2.6.1. Steel-concrete slip evaluation

The longitudinal slips of the two steel bars are independently calculated with respect to concrete. For the x bar:

$$u_{gx}(x) = u_{sx}(x) - u_{cx}(x) \qquad (49)$$

where:

$$u_{cx}(x) = \underbrace{[1 \ 0 \ 0 \ 0 \ z_{ax} \ 0]}_{a_x} N_p u_c = a_x N_c u_{el} \qquad (50)$$

which gives:

$$u_{gx}(x) = N_{sx}u_{el} - a_x N_c u_{el} = (N_{sx} - a_x N_c) u_{el} \qquad (51)$$

Similarly, the y steel bar longitudinal slip with respect to concrete is evaluated as follows:

$$u_{gy}(y) = u_{sy}(y) - u_{cy}(y) = (N_{sy} - a_y N_c) u_{el} \qquad (52)$$

with:

$$a_y = [0 \ 1 \ 0 \ -z_{ay} \ 0 \ 0] \qquad (53)$$

## 2.6.2. Weak formulation

The virtual power principle of the enhanced plate element is expressed as follows:

$$\begin{aligned} & \mathbf{u}_{el}^* \mathbf{F}_{intc}^e + \int_0^{L_x} \epsilon_{sx}^* \sigma_{sx} (\epsilon_{sx}) \epsilon_{sx} dx + \int_0^{L_y} \epsilon_{sy}^* \sigma_{sy} (\epsilon_{sy}) \epsilon_{sy} dy \\ & + \int_0^{L_x} u_{gx}^* \tau_x (u_{gx}) P_x dx + \int_0^{L_y} u_{gy}^* \tau_y (u_{gy}) P_y dy = \mathbf{u}_{el}^* \mathbf{F}_e \end{aligned} \qquad (54)$$

which gives after a simplification with respect to $\mathbf{u}_{el}^* \mathbf{F}_e$, the following expression:

$$\begin{aligned} & \mathbf{F}_{intc}^e + \int_0^{L_x} \mathbf{B}_{sx}^T \sigma_{sx} (\epsilon_{sx}) \epsilon_{sx} dx + \int_0^{L_y} \mathbf{B}_{sy}^T \sigma_{sy} (\epsilon_{sy}) \epsilon_{sy} dy \\ & + \int_0^{L_x} (N_{sx}^T - N_c^T a_x^T) \tau_x (u_{gx}) \mathbf{u}_{el} P_x dx + \int_0^{L_y} (N_{sy}^T - N_c^T a_y^T) \tau_y (u_{gy}) \mathbf{u}_{el} P_y dy = \mathbf{F}_e \end{aligned} \qquad (55)$$

where $\tau_x$ and $\tau_y$ are the bond stresses around the x and y bars of perimeters $P_x$ and $P_y$, respectively. The entirety of the expression on the left-hand side of Eq. (55) denotes for the internal forces vector $\mathbf{F}_e^i$ of the enhanced plate element. Its derivation with respect to $\mathbf{u}_{el}$ gives the tangent operator of the enhanced element:

$$\frac{\partial \mathbf{F}_e^e}{\partial \mathbf{u}_{el}} = k_{sc}^e = k_c^e + k_{sx}^e + k_{sy}^e + k_{gx}^e + k_{gy}^e \qquad (56)$$

$k_c^e$ is the tangent operator of the concrete plate element. $k_{sx}^e$ and $k_{sy}^e$ are the tangent operators of the two steel bars. $k_{gx}^e$ and $k_{gy}^e$ are the tangent operators due to the steel-concrete bond of the two steel bars, calculated as follows:

$$k_{gx}^e = \int_0^{L_x} (N_{sx}^T - N_c^T a_x^T) \frac{\partial \tau_x (u_{gx})}{\partial \mathbf{u}_{el}} P_x dx \qquad (57)$$

$$k_{gy}^e = \int_0^{L_y} (N_{sy}^T - N_c^T a_y^T) \frac{\partial \tau_y (u_{gy})}{\partial \mathbf{u}_{el}} P_y dy \qquad (58)$$

where:

$$\frac{\partial \tau_x (u_{gx})}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_x (u_{gx})}{\partial u_{gx}} \frac{\partial u_{gx}}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_x (u_{gx})}{\partial u_{gx}} (N_{sx} - a_x N_c) \qquad (59)$$

$$\frac{\partial \tau_y (u_{gy})}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_y (u_{gy})}{\partial u_{gy}} \frac{\partial u_{gy}}{\partial \mathbf{u}_{el}} = \frac{\partial \tau_y (u_{gy})}{\partial u_{gy}} (N_{sy} - a_y N_c) \qquad (60)$$

knowing that $\frac{\partial \tau_x (u_{gx})}{\partial u_{gx}}$ and $\frac{\partial \tau_y (u_{gy})}{\partial u_{gy}}$ are the slopes of the bond laws considered for bond stresses around two steel bars. The formulation allows the utilization of two distinct laws.

# 3. Verification tests

## 3.1. Pull-out test

In this subsection, the pull-out experimental test of [55] is modeled.
<figure><img src="image_13.png" /><figcaption>Fig. 9. Normalized pull-out test.</figcaption></figure>

<table><caption>**Table 1**
Material properties for the pull-out model.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Concrete Young's modulus</td><td>28</td><td>GPa</td></tr><tr><td>E<sub>s</sub></td><td>Steel Young's modulus</td><td>200</td><td>GPa</td></tr><tr><td>τ<sub>1</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>22.5</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>1.45</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>10</td><td>mm</td></tr></tbody></table>

### 3.1.1. Test description

The test setup consists of a concrete cube with a single steel reinforcement bar passing through it. The cube's movement is restricted by a metal plate. A Teflon support ensures proper alignment between the steel bar and the direction of movement, preventing bending. Two LVDT (Linear Variable Differential Transformer) are used to measure the steel-concrete slip at the two bar edges. The contact length is 5 times the steel bar's diameter, a value recommended in [56]. The bond stress $\tau$ along the steel bar is estimated as follows:

$$ \tau = \frac{F}{d_a L \pi} \qquad (61) $$

where $F$ is the measured reaction, $L$ is the steel-concrete contact length, and $d_a$ is the steel bar diameter. The experimental bond law is defined as the evolution of the calculated bond stress value $\tau$ with respect to the slip measured at the unloaded edge of the steel bar (see Fig. 9).

### 3.1.2. Material properties

Linear elastic constitutive laws are assigned to both steel and concrete, intentionally focusing on the pull-out nonlinearities at the interface. The only source of non-linearity in this scenario arises from the non-linear expression of the bond law at the steel-concrete interface.

The bond law of [54] is used (see Fig. 6). Table 1 sums up the used material properties. It has to be noted that the bond law parameters are chosen to best describe the experimental bond law curve of [55] (see Fig. 10).

### 3.1.3. Finite elements mesh

The central part of the pull-out test where a bond contact links steel and concrete is represented with a bond length equal to five times the bar diameter. Fig. 11 shows the mesh and the boundary conditions. Three enhanced beam elements are used to construct the mesh. Each enhanced element is an assembly of: two Euler-Bernoulli concrete beams, one three-node steel bar element, and bond stresses.

Monotonic and cyclic load configurations are tested. The boundary conditions of Fig. 11 are adopted for a monotonic loading path. For the cyclic case, boundary conditions are modified every time the sign of the imposed displacement changes, by switching the boundary conditions of nodes 1 and 7 of Fig. 11.

### 3.1.4. Results

Fig. 12 illustrates the imposed displacement and the reaction curves for the monotonic applied load path. Steel and concrete nodes displacements for the different time steps are presented in Fig. 13. The displacements of Fig. 13 show that the steel bar is successfully being pulled with respect to concrete. Even though no experimental result is given in [55] for cyclic pull-out tests, the law of [54] can reproduce the cyclic behavior with no additional input parameters. Fig. 14 shows a cyclic applied displacement path and the resultant reaction curve.
<figure><img src="image_13.png" /><figcaption>Fig. 9. Normalized pull-out test.</figcaption></figure>

<table><caption>**Table 1**
Material properties for the pull-out model.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Concrete Young's modulus</td><td>28</td><td>GPa</td></tr><tr><td>E<sub>s</sub></td><td>Steel Young's modulus</td><td>200</td><td>GPa</td></tr><tr><td>τ<sub>1</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>22.5</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>1.45</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>Input parameter of the bond law (Fig. 6(a))</td><td>10</td><td>mm</td></tr></tbody></table>

### 3.1.1. Test description

The test setup consists of a concrete cube with a single steel reinforcement bar passing through it. The cube's movement is restricted by a metal plate. A Teflon support ensures proper alignment between the steel bar and the direction of movement, preventing bending. Two LVDT (Linear Variable Differential Transformer) are used to measure the steel-concrete slip at the two bar edges. The contact length is 5 times the steel bar's diameter, a value recommended in [56]. The bond stress $\tau$ along the steel bar is estimated as follows:

$$ \tau = \frac{F}{d_a L \pi} \qquad (61) $$

where $F$ is the measured reaction, $L$ is the steel-concrete contact length, and $d_a$ is the steel bar diameter. The experimental bond law is defined as the evolution of the calculated bond stress value $\tau$ with respect to the slip measured at the unloaded edge of the steel bar (see Fig. 9).

### 3.1.2. Material properties

Linear elastic constitutive laws are assigned to both steel and concrete, intentionally focusing on the pull-out nonlinearities at the interface. The only source of non-linearity in this scenario arises from the non-linear expression of the bond law at the steel-concrete interface.

The bond law of [54] is used (see Fig. 6). Table 1 sums up the used material properties. It has to be noted that the bond law parameters are chosen to best describe the experimental bond law curve of [55] (see Fig. 10).

### 3.1.3. Finite elements mesh

The central part of the pull-out test where a bond contact links steel and concrete is represented with a bond length equal to five times the bar diameter. Fig. 11 shows the mesh and the boundary conditions. Three enhanced beam elements are used to construct the mesh. Each enhanced element is an assembly of: two Euler-Bernoulli concrete beams, one three-node steel bar element, and bond stresses.

Monotonic and cyclic load configurations are tested. The boundary conditions of Fig. 11 are adopted for a monotonic loading path. For the cyclic case, boundary conditions are modified every time the sign of the imposed displacement changes, by switching the boundary conditions of nodes 1 and 7 of Fig. 11.

### 3.1.4. Results

Fig. 12 illustrates the imposed displacement and the reaction curves for the monotonic applied load path. Steel and concrete nodes displacements for the different time steps are presented in Fig. 13. The displacements of Fig. 13 show that the steel bar is successfully being pulled with respect to concrete. Even though no experimental result is given in [55] for cyclic pull-out tests, the law of [54] can reproduce the cyclic behavior with no additional input parameters. Fig. 14 shows a cyclic applied displacement path and the resultant reaction curve.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Analytical and experimental bond laws.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Considered mesh and boundary conditions for the pull-out model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Monotonic 1D pull-out: imposed displacement (a) and reaction (b).</figcaption></figure>

### 3.1.5. Discussion

The reaction force curves show that the constitutive expression that links the reaction force *F* to the bond stress *τ* (Eq. (61)) is fulfilled, which validates the enhancement approach.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Analytical and experimental bond laws.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Considered mesh and boundary conditions for the pull-out model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Monotonic 1D pull-out: imposed displacement (a) and reaction (b).</figcaption></figure>

### 3.1.5. Discussion

The reaction force curves show that the constitutive expression that links the reaction force *F* to the bond stress *τ* (Eq. (61)) is fulfilled, which validates the enhancement approach.
<figure><img src="image_17.png" /><figcaption>Fig. 13. Concrete (a) and steel (b) nodes displacements with respect to their longitudinal position, for the used time steps.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Cyclic pull-out test: imposed displacement (a) and reaction (b).</figcaption></figure>

Both steel-concrete slip interpolation methods presented in Section 2.5 were tested, producing the same results, for this specific application. Using three-node steel bars instead of regular truss elements usually allows for better interpolation of the steel-concrete slip field due to their nonlinear shape functions, which are of higher order compared to two-node bars.

Furthermore, the enhancement methodology that assembles distinct finite elements for steel and concrete gives the possibility to impose a steel displacement, independently of the concrete displacement, which is the case for this pull-out test.

## 3.2. Enhanced plate membrane behavior test

This test examines the membrane behavior of a reinforced plate element, which includes an eight-node plate concrete element and a truss steel bar element in the direction of the x axis. The bond between the steel bar and concrete is duly incorporated within the enhanced element.

### 3.2.1. Test description

A concrete square plate with dimensions 10 × 10 m² and a thickness of 0.25 m is considered. The steel bar, with a diameter of 16 mm, is situated in the plate's mid-plane. The focus of this study is on the plate membrane behavior.

### 3.2.2. Material properties

A linear behavior is attributed to both the steel and the concrete, with the properties outlined in Table 2.

### 3.2.3. Finite elements mesh

The mesh is composed of a single enhanced element. The boundary conditions described in Fig. 15 are considered, where an imposed displacement is applied to the edge of the steel bar. A linear steel-concrete bond law is considered.
<figure><img src="image_17.png" /><figcaption>Fig. 13. Concrete (a) and steel (b) nodes displacements with respect to their longitudinal position, for the used time steps.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Cyclic pull-out test: imposed displacement (a) and reaction (b).</figcaption></figure>

Both steel-concrete slip interpolation methods presented in Section 2.5 were tested, producing the same results, for this specific application. Using three-node steel bars instead of regular truss elements usually allows for better interpolation of the steel-concrete slip field due to their nonlinear shape functions, which are of higher order compared to two-node bars.

Furthermore, the enhancement methodology that assembles distinct finite elements for steel and concrete gives the possibility to impose a steel displacement, independently of the concrete displacement, which is the case for this pull-out test.

## 3.2. Enhanced plate membrane behavior test

This test examines the membrane behavior of a reinforced plate element, which includes an eight-node plate concrete element and a truss steel bar element in the direction of the x axis. The bond between the steel bar and concrete is duly incorporated within the enhanced element.

### 3.2.1. Test description

A concrete square plate with dimensions 10 × 10 m² and a thickness of 0.25 m is considered. The steel bar, with a diameter of 16 mm, is situated in the plate's mid-plane. The focus of this study is on the plate membrane behavior.

### 3.2.2. Material properties

A linear behavior is attributed to both the steel and the concrete, with the properties outlined in Table 2.

### 3.2.3. Finite elements mesh

The mesh is composed of a single enhanced element. The boundary conditions described in Fig. 15 are considered, where an imposed displacement is applied to the edge of the steel bar. A linear steel-concrete bond law is considered.
<table><caption>**Table 2**
Steel and concrete Parameters of the enhanced plate membrane test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Concrete Young's modulus</td><td>30</td><td>GPa</td></tr><tr><td>E<sub>s</sub></td><td>Steel Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>ν<sub>c</sub></td><td>Concrete Poisson's ratio</td><td>0.22</td><td>-</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Boundary conditions of the enhanced plate membrane test.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Deformed enhanced plate shape for different stiffness values of the linear bond law: (a) $10^5$ Pa/m, (b) $10^{10}$ Pa/m, (c) $10^{20}$ Pa/m.</figcaption></figure>

### 3.2.4. Result

A parametric study is conducted to examine the variation in stiffness of the bond law, with the aim of observing its impact on the deformed shape of the plate, as shown in Fig. 16.

### 3.2.5. Discussion

The deformed plate shapes observed in Fig. 16 align with theoretical expectations. Indeed, an increase in the slope of the bond law leads to a diminished sliding of the steel with respect to concrete, resulting in a decrease in the relative displacement between the two materials. Consequently, by imposing a displacement at the steel edge, a higher slope of the bond law facilitates a more pronounced carrying effect of the concrete by the steel (resulting in a smaller sliding motion).

## 3.3. Enhanced plate flexural behavior test

This test aims to analyze the flexural behavior of a concrete plate reinforced with steel bars aligned in both the x and y directions.

### 3.3.1. Test description

The considered test is illustrated in Fig. 17. The plate has a thickness of 0.2 m. All the bars are eccentrically positioned in the downward direction along the z axis by a distance of 0.05 m, with respect to the concrete plate mid-plane, and have a diameter of 16 mm.

### 3.3.2. Material properties

A linear behavior is assigned to steel and concrete, with the parameters of Table 3. The bond law of [54] is used, with the parameters of Table 4.
<table><caption>**Table 2**
Steel and concrete Parameters of the enhanced plate membrane test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Concrete Young's modulus</td><td>30</td><td>GPa</td></tr><tr><td>E<sub>s</sub></td><td>Steel Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>v<sub>c</sub></td><td>Concrete Poisson's ratio</td><td>0.22</td><td>-</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>Fig. 15. Boundary conditions of the enhanced plate membrane test.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Deformed enhanced plate shape for different stiffness values of the linear bond law: (a) $10^5$ Pa/m, (b) $10^{10}$ Pa/m, (c) $10^{20}$ Pa/m.</figcaption></figure>

### 3.2.4. Result

A parametric study is conducted to examine the variation in stiffness of the bond law, with the aim of observing its impact on the deformed shape of the plate, as shown in Fig. 16.

### 3.2.5. Discussion

The deformed plate shapes observed in Fig. 16 align with theoretical expectations. Indeed, an increase in the slope of the bond law leads to a diminished sliding of the steel with respect to concrete, resulting in a decrease in the relative displacement between the two materials. Consequently, by imposing a displacement at the steel edge, a higher slope of the bond law facilitates a more pronounced carrying effect of the concrete by the steel (resulting in a smaller sliding motion).

## 3.3. Enhanced plate flexural behavior test

This test aims to analyze the flexural behavior of a concrete plate reinforced with steel bars aligned in both the x and y directions.

### 3.3.1. Test description

The considered test is illustrated in Fig. 17. The plate has a thickness of 0.2 m. All the bars are eccentrically positioned in the downward direction along the z axis by a distance of 0.05 m, with respect to the concrete plate mid-plane, and have a diameter of 16 mm.

### 3.3.2. Material properties

A linear behavior is assigned to steel and concrete, with the parameters of Table 3. The bond law of [54] is used, with the parameters of Table 4.
<table><caption>**Table 3**

Steel and concrete material parameters for the plate bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>b</sub></td><td>Concrete Young's modulus</td><td>30</td><td>GPa</td></tr><tr><td>E<sub>a</sub></td><td>Steel Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>ν<sub>b</sub></td><td>Concrete Poisson's ratio</td><td>0.2</td><td>-</td></tr></tbody></table>

<table><caption>**Table 4**

Bond law parameters for the plate bending test.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>τ<sub>1</sub></td><td>10</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>1</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>8</td><td>mm</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Finite elements mesh of the plate bending test.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** Plate bending test reaction curve (a) and corresponding deformed shape (b).</figcaption></figure>

### 3.3.3. Finite elements mesh

The considered mesh is composed of 10 × 10 enhanced plate elements, composed each of: an eight-node concrete plate element, two steel truss elements along x and y, and steel-concrete bond stresses. The plate is simply supported (in the z direction) along its 4 edges. A displacement in the z direction (downward) is applied at the central node of the plate.

### 3.3.4. Results

Fig. 18(a) shows the plate reaction curve and Fig. 18(b) illustrates the corresponding deformed shape of the plate mid-plane.
<table><caption>**Table 3**

Steel and concrete material parameters for the plate bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>b</sub></td><td>Concrete Young's modulus</td><td>30</td><td>GPa</td></tr><tr><td>E<sub>a</sub></td><td>Steel Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>ν<sub>b</sub></td><td>Concrete Poisson's ratio</td><td>0.2</td><td>-</td></tr></tbody></table>

<table><caption>**Table 4**

Bond law parameters for the plate bending test.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>τ<sub>1</sub></td><td>10</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>1</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>8</td><td>mm</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Finite elements mesh of the plate bending test.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** Plate bending test reaction curve (a) and corresponding deformed shape (b).</figcaption></figure>

### 3.3.3. Finite elements mesh

The considered mesh is composed of 10 × 10 enhanced plate elements, composed each of: an eight-node concrete plate element, two steel truss elements along x and y, and steel-concrete bond stresses. The plate is simply supported (in the z direction) along its 4 edges. A displacement in the z direction (downward) is applied at the central node of the plate.

### 3.3.4. Results

Fig. 18(a) shows the plate reaction curve and Fig. 18(b) illustrates the corresponding deformed shape of the plate mid-plane.
<figure><img src="image_23.png" /><figcaption>Fig. 19. Slip values of x (a) and y (b) steel bars, for the different time steps of the solution.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Configuration of the four-point beam test of [57].</figcaption></figure>

Fig. 19 shows the slip values of the steel bars along the x and y directions, with respect to concrete, for the different time steps of the solution.

### 3.3.5. Discussion

The slip values are equal in both x and y directions, and they are symmetrical with respect to the plate mid-plane central axes. This observation is consistent, considering that the plate geometry, the boundary conditions, and the applied loading are symmetrical. This test validates the out-of-plane behavior of the enhanced plate formulation, which is reinforced with steel bars capable of sliding, in both plate directions, with respect to concrete.

## 4. Numerical application

### 4.1. Four-point beam bending test

A four-point bending test of a reinforced concrete beam conducted in [57] is considered.

#### 4.1.1. Test description

The considered test is described in Fig. 20.

Two distinct modeling options are chosen for this test using enhanced plate and beam elements. The modeling option with enhanced concrete beams defines each enhanced element as a combination of three contributions: one concrete Timoshenko fiber beam element, one steel truss element, and steel-concrete bond stresses. The modeling option with enhanced concrete plates defines each enhanced element as a combination of three contributions: a Mindlin-Reissner concrete plate element, one steel truss element, and steel-concrete bond stresses.

#### 4.1.2. Material properties

An elastic plastic behavior with a linear strain hardening is considered for the steel behavior, with the parameters of Table 5.

For the enhanced concrete beam elements modeling choice, a 1D traction damage law is accorded to the concrete beam fibers. This law is detailed in Appendix C. The parameters of the law are presented in Table 6.
<figure><img src="image_23.png" /><figcaption>Fig. 19. Slip values of x (a) and y (b) steel bars, for the different time steps of the solution.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Configuration of the four-point beam test of [57].</figcaption></figure>

Fig. 19 shows the slip values of the steel bars along the x and y directions, with respect to concrete, for the different time steps of the solution.

### 3.3.5. Discussion

The slip values are equal in both x and y directions, and they are symmetrical with respect to the plate mid-plane central axes. This observation is consistent, considering that the plate geometry, the boundary conditions, and the applied loading are symmetrical. This test validates the out-of-plane behavior of the enhanced plate formulation, which is reinforced with steel bars capable of sliding, in both plate directions, with respect to concrete.

## 4. Numerical application

### 4.1. Four-point beam bending test

A four-point bending test of a reinforced concrete beam conducted in [57] is considered.

#### 4.1.1. Test description

The considered test is described in Fig. 20.

Two distinct modeling options are chosen for this test using enhanced plate and beam elements. The modeling option with enhanced concrete beams defines each enhanced element as a combination of three contributions: one concrete Timoshenko fiber beam element, one steel truss element, and steel-concrete bond stresses. The modeling option with enhanced concrete plates defines each enhanced element as a combination of three contributions: a Mindlin-Reissner concrete plate element, one steel truss element, and steel-concrete bond stresses.

#### 4.1.2. Material properties

An elastic plastic behavior with a linear strain hardening is considered for the steel behavior, with the parameters of Table 5.

For the enhanced concrete beam elements modeling choice, a 1D traction damage law is accorded to the concrete beam fibers. This law is detailed in Appendix C. The parameters of the law are presented in Table 6.
<table><caption>**Table 5**

Steel parameters: four-point flexural beam test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>a</sub></td><td>Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>f<sub>y</sub></td><td>Limit of elasticity</td><td>450</td><td>MPa</td></tr><tr><td>E<sub>h</sub></td><td>Hardening slope</td><td>2100</td><td>MPa</td></tr></tbody></table>

<table><caption>**Table 6**

Concrete nonlinear parameters: enhanced beams model of the four-point bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Young's modulus</td><td>28</td><td>GPa</td></tr><tr><td>f<sub>t</sub></td><td>Tensile strength</td><td>3.2</td><td>MPa</td></tr><tr><td>ε<sub>d0</sub></td><td>Damage threshold</td><td>f<sub>L</sub>/E<sub>c</sub> = 1.1429 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>B</td><td>Parameter of the damage law</td><td>1340</td><td>-</td></tr><tr><td>G<sub>f</sub></td><td>Fracture energy</td><td>150</td><td>N/m</td></tr></tbody></table>

<table><caption>**Table 7**

Concrete nonlinear parameters: enhanced plates model of the four-point bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>f<sub>t</sub></td><td>Tensile strength</td><td>2.6</td><td>MPa</td></tr><tr><td>f<sub>c</sub></td><td>Compressive strength</td><td>56.9</td><td>MPa</td></tr><tr><td>ε<sub>i0</sub></td><td>Revised Mazars' model input</td><td>f<sub>L</sub>/E<sub>c</sub> = 1.1429 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>ε<sub>c0</sub></td><td>Revised Mazars' model input</td><td>f<sub>L</sub>/E<sub>c</sub> = 1.9 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>A<sub>1</sub></td><td>Revised Mazars' model input</td><td>0.99</td><td>-</td></tr><tr><td>B<sub>1</sub></td><td>Revised Mazars' model input</td><td>8000</td><td>-</td></tr><tr><td>A<sub>c</sub></td><td>Revised Mazars' model input</td><td>1.2</td><td>-</td></tr><tr><td>B<sub>c</sub></td><td>Revised Mazars' model input</td><td>400</td><td>-</td></tr><tr><td>β</td><td>Revised Mazars' model input</td><td>1.06</td><td>-</td></tr></tbody></table>

For the enhanced concrete plate elements modeling choice, a revised Mazars' law is considered for concrete [58,59] (see Appendix D) with the parameters of Table 7.

The bond law of [54] is used. In order to determine the parameters of this law, the failure mode of the interface is predicted according to the criterion of [55]. Using this criterion, the ratio $c/d$ is calculated, where c and d represent the concrete cover and the steel bar diameter, respectively. This ratio is then compared to the value of $0.39f_c/f_t - 0.24$. In this case study, one can say:

$$ \frac{c}{d} < 0.39 \frac{f_c}{f_t} - 0.24 : \text{ so concrete cover splitting failure} \qquad (62) $$

Knowing that $f_c$ is the concrete compressive strength. Since the predicted failure mode is by concrete splitting, the parameters $\tau_1$ and $g_1$ of the bond law, that represent the maximum bond stress and the corresponding slip, are estimated as follows:

$$ g_1 = 0.17 \frac{c}{d} = 0.17 \times 2.5 = 0.425 \text{ mm} \qquad (63) $$

$$ \tau_1 = f_t \left( 1.53 \frac{c}{d} + 0.36 \right) = 3.2(1.53 \times 2.5 + 0.36) = 13.392 \text{ Pa} \qquad (64) $$

Since the work of [55] does not provide specifications about the calculation of the third defining parameter $g_3$ of the bond law, the recommendations of [60] are followed for the identification of this parameter. For unconfined concrete around the steel (which is the case for this test that does not include steel stirrups in the beam volume), $g_3$ is estimated to be equal to $1.2g_1$, so for this example it is equal to 0.51 mm.

### 4.1.3. Finite elements mesh

Fig. 21(a) illustrates the various fibers that constitute the concrete mesh of the enhanced fiber beams. The boundary conditions of the enhanced fiber beams are detailed in Fig. 21(b).

The finite elements mesh and boundary conditions of the modeling choice carried out with enhanced plates are illustrated in Fig. 22. This modeling choice focuses on the membrane plate elements behavior. The mesh construction is done using enhanced plates at the level of the steel bar. The rest of the structure is meshed with standard concrete plate elements with no enhancement.

### 4.1.4. Results

Fig. 23 illustrates the reaction curves for both modeling choices (enhanced beams and plates) with: a perfect steel-concrete bond, and a non-linear interface behavior (imperfect bond). These curves are compared to the experimental reaction curve. These reaction curves are compared to the numerical curves of [25] where this same test is modeled.
<table><caption>**Table 5**

Steel parameters: four-point flexural beam test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>a</sub></td><td>Young's modulus</td><td>210</td><td>GPa</td></tr><tr><td>f<sub>y</sub></td><td>Limit of elasticity</td><td>450</td><td>MPa</td></tr><tr><td>E<sub>h</sub></td><td>Hardening slope</td><td>2100</td><td>MPa</td></tr></tbody></table>

<table><caption>**Table 6**

Concrete nonlinear parameters: enhanced beams model of the four-point bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Young's modulus</td><td>28</td><td>GPa</td></tr><tr><td>f<sub>t</sub></td><td>Tensile strength</td><td>3.2</td><td>MPa</td></tr><tr><td>ε<sub>0</sub></td><td>Damage threshold</td><td>1.1429 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>B</td><td>Parameter of the damage law</td><td>1340</td><td>-</td></tr><tr><td>G<sub>f</sub></td><td>Fracture energy</td><td>150</td><td>N/m</td></tr></tbody></table>

<table><caption>**Table 7**

Concrete nonlinear parameters: enhanced plates model of the four-point bending test.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>f<sub>t</sub></td><td>Tensile strength</td><td>2.6</td><td>MPa</td></tr><tr><td>f<sub>c</sub></td><td>Compressive strength</td><td>56.9</td><td>MPa</td></tr><tr><td>ε<sub>i0</sub></td><td>Revised Mazars' model input</td><td>1.1429 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>ε<sub>c0</sub></td><td>Revised Mazars' model input</td><td>1.9 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>A<sub>i</sub></td><td>Revised Mazars' model input</td><td>0.99</td><td>-</td></tr><tr><td>B<sub>i</sub></td><td>Revised Mazars' model input</td><td>8000</td><td>-</td></tr><tr><td>A<sub>c</sub></td><td>Revised Mazars' model input</td><td>1.2</td><td>-</td></tr><tr><td>B<sub>c</sub></td><td>Revised Mazars' model input</td><td>400</td><td>-</td></tr><tr><td>β</td><td>Revised Mazars' model input</td><td>1.06</td><td>-</td></tr></tbody></table>

For the enhanced concrete plate elements modeling choice, a revised Mazars' law is considered for concrete [58,59] (see Appendix D) with the parameters of Table 7.

The bond law of [54] is used. In order to determine the parameters of this law, the failure mode of the interface is predicted according to the criterion of [55]. Using this criterion, the ratio $\frac{c}{d}$ is calculated, where c and d represent the concrete cover and the steel bar diameter, respectively. This ratio is then compared to the value of 0.39 $\frac{f_c}{f_t}$ - 0.24. In this case study, one can say:

$$ \frac{c}{d} < 0.39 \frac{f_c}{f_t} - 0.24 : \text{ so concrete cover splitting failure} \qquad (62) $$

Knowing that $f_c$ is the concrete compressive strength. Since the predicted failure mode is by concrete splitting, the parameters $\tau_1$ and $g_1$ of the bond law, that represent the maximum bond stress and the corresponding slip, are estimated as follows:

$$ g_1 = 0.17 \frac{c}{d} = 0.17 \times 2.5 = 0.425 \text{ mm} \qquad (63) $$

$$ \tau_1 = f_t \left( 1.53 \frac{c}{d} + 0.36 \right) = 3.2(1.53 \times 2.5 + 0.36) = 13.392 \text{ Pa} \qquad (64) $$

Since the work of [55] does not provide specifications about the calculation of the third defining parameter $g_3$ of the bond law, the recommendations of [60] are followed for the identification of this parameter. For unconfined concrete around the steel (which is the case for this test that does not include steel stirrups in the beam volume), $g_3$ is estimated to be equal to 1.2$g_1$, so for this example it is equal to 0.51 mm.

### 4.1.3. Finite elements mesh

Fig. 21(a) illustrates the various fibers that constitute the concrete mesh of the enhanced fiber beams. The boundary conditions of the enhanced fiber beams are detailed in Fig. 21(b).

The finite elements mesh and boundary conditions of the modeling choice carried out with enhanced plates are illustrated in Fig. 22. This modeling choice focuses on the membrane plate elements behavior. The mesh construction is done using enhanced plates at the level of the steel bar. The rest of the structure is meshed with standard concrete plate elements with no enhancement.

### 4.1.4. Results

Fig. 23 illustrates the reaction curves for both modeling choices (enhanced beams and plates) with: a perfect steel-concrete bond, and a non-linear interface behavior (imperfect bond). These curves are compared to the experimental reaction curve. These reaction curves are compared to the numerical curves of [25] where this same test is modeled.
<figure><img src="image_25.png" /><figcaption>Fig. 21. Concrete fibers mesh (a) and boundary conditions: four-point bending test model with enhanced concrete beams.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 22. Finite elements mesh and boundary conditions: four-point bending test model with enhanced concrete plates.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 23. Reaction curves of the four-point bending beam test with enhanced beams (a) and plates (b).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 21. Concrete fibers mesh (a) and boundary conditions: four-point bending test model with enhanced concrete beams.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 22. Finite elements mesh and boundary conditions: four-point bending test model with enhanced concrete plates.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 23. Reaction curves of the four-point bending beam test with enhanced beams (a) and plates (b).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 24. Numerical damage field distributions compared to the experimental cracking pattern.</figcaption></figure>

Considering the nonlinear behavior of the interface does not significantly impact the reaction curve. Nevertheless, accounting for a nonlinear interface behavior leads to a more accurate representation of the experimental concrete cracking pattern, as depicted in Fig. 24. Specifically, in the vicinity of the steel, the simulations that incorporate a nonlinear interface behavior show improved detection of the concrete cracks localization and spacing. A maximum of 15 and 16 iterations are needed for beam and plate models respectively, while considering an imperfect bond condition. For perfect bond simulations, a maximum of 20 and 35 iterations is respectively needed. These iterations numbers are less than those typically seen in other models discussed in the literature, such as [25]. Indeed, the number of iterations needed in [25] is between 200 and 1100 iterations.

### 4.1.5. Discussion

Fig. 24 reveals a notable difference in the damage patterns with and without considering the steel-concrete interface behavior for the enhanced beams modeling choice. This contrast is subdued when using enhanced plate elements. Indeed, the assumption of a beam plane section that remains plane after deformation, adopted for concrete beam elements, spreads the damage within the concrete beam volume. Considering the nonlinear interface behavior helps overcome this limitation of beam elements. The kinematics of plates are less constrained compared to beam kinematics, resulting in a less significant difference in damage patterns with and without considering the interface nonlinear behavior.

The local description of the concrete cracking process is characterized by a steel stress localization and a discontinuity in the slip field, as shown in Figs. 25 and 26 that present the steel stress and the slip values for the resolution time steps. The sliding in these figures is calculated at the nodes of the steel truss elements. The steel stresses are calculated at the Gauss points of these elements. These figures show the progressive apparition of the concrete cracks, detect by the time steps where a steel stress concentration and a steel-concrete slip occur at the crack position, near the steel reinforcement.

## 5. Conclusive remarks and perspectives

This paper proposes a modeling method for steel-concrete interface behavior in semi-global beam and plate finite elements. It involves creating an enhanced element gathering concrete and steel elements, along with steel-concrete bond stresses, allowing the use of existing finite elements library components by defining additional steel degrees of freedom. The static equilibrium of the enhanced element is performed in the framework of principle of virtual work. The methodology allows the construction of a variety of enhanced elements configurations, knowing that for configurations where inner nodes are defined a static condensation technique is applied to the enhanced element.

The paper's modeling technique draws inspiration from prior literature that focused on modeling interface behavior in fiber beam elements [25,26]. In contrast, this paper presents a versatile approach that leverages existing classical finite elements for plates and beams. These elements are assembled to form a single enhanced reinforced concrete element, simplifying both formulation and implementation.

Three distinct examples of enhanced elements assemblies are presented with: one beam concrete element and one truss steel element, two beam concrete elements and one three-node steel bar element, and one concrete plate element with two steel bar elements. This paper presents three validation examples testing the behavior of enhanced beams and plate elements (membrane and out-of-plane plate behaviors). These tests are followed by an application example where a four-point beam bending test is independently modeled with enhanced beam and plate elements. The numerical results show a good agreement with experimental and numerical results of the literature. The different shown applications allow a quantitative evaluation of steel-concrete slip and bond stress values.

The perspectives of future research is toward testing more capabilities of the enhancement approach, by upgrading its current version. More specifically, using an assembly of multilayered concrete plate elements will allow the study of the out-of-plane
<figure><img src="image_28.png" /><figcaption>Fig. 24. Numerical damage field distributions compared to the experimental cracking pattern.</figcaption></figure>

Considering the nonlinear behavior of the interface does not significantly impact the reaction curve. Nevertheless, accounting for a nonlinear interface behavior leads to a more accurate representation of the experimental concrete cracking pattern, as depicted in Fig. 24. Specifically, in the vicinity of the steel, the simulations that incorporate a nonlinear interface behavior show improved detection of the concrete cracks localization and spacing. A maximum of 15 and 16 iterations are needed for beam and plate models respectively, while considering an imperfect bond condition. For perfect bond simulations, a maximum of 20 and 35 iterations is respectively needed. These iterations numbers are less than those typically seen in other models discussed in the literature, such as [25]. Indeed, the number of iterations needed in [25] is between 200 and 1100 iterations.

### 4.1.5. Discussion

Fig. 24 reveals a notable difference in the damage patterns with and without considering the steel-concrete interface behavior for the enhanced beams modeling choice. This contrast is subdued when using enhanced plate elements. Indeed, the assumption of a beam plane section that remains plane after deformation, adopted for concrete beam elements, spreads the damage within the concrete beam volume. Considering the nonlinear interface behavior helps overcome this limitation of beam elements. The kinematics of plates are less constrained compared to beam kinematics, resulting in a less significant difference in damage patterns with and without considering the interface nonlinear behavior.

The local description of the concrete cracking process is characterized by a steel stress localization and a discontinuity in the slip field, as shown in Figs. 25 and 26 that present the steel stress and the slip values for the resolution time steps. The sliding in these figures is calculated at the nodes of the steel truss elements. The steel stresses are calculated at the Gauss points of these elements. These figures show the progressive apparition of the concrete cracks, detect by the time steps where a steel stress concentration and a steel-concrete slip occur at the crack position, near the steel reinforcement.

## 5. Conclusive remarks and perspectives

This paper proposes a modeling method for steel-concrete interface behavior in semi-global beam and plate finite elements. It involves creating an enhanced element gathering concrete and steel elements, along with steel-concrete bond stresses, allowing the use of existing finite elements library components by defining additional steel degrees of freedom. The static equilibrium of the enhanced element is performed in the framework of principle of virtual work. The methodology allows the construction of a variety of enhanced elements configurations, knowing that for configurations where inner nodes are defined a static condensation technique is applied to the enhanced element.

The paper's modeling technique draws inspiration from prior literature that focused on modeling interface behavior in fiber beam elements [25,26]. In contrast, this paper presents a versatile approach that leverages existing classical finite elements for plates and beams. These elements are assembled to form a single enhanced reinforced concrete element, simplifying both formulation and implementation.

Three distinct examples of enhanced elements assemblies are presented with: one beam concrete element and one truss steel element, two beam concrete elements and one three-node steel bar element, and one concrete plate element with two steel bar elements. This paper presents three validation examples testing the behavior of enhanced beams and plate elements (membrane and out-of-plane plate behaviors). These tests are followed by an application example where a four-point beam bending test is independently modeled with enhanced beam and plate elements. The numerical results show a good agreement with experimental and numerical results of the literature. The different shown applications allow a quantitative evaluation of steel-concrete slip and bond stress values.

The perspectives of future research is toward testing more capabilities of the enhancement approach, by upgrading its current version. More specifically, using an assembly of multilayered concrete plate elements will allow the study of the out-of-plane
<figure><img src="image_29.png" /><figcaption>Fig. 25. Local cracking description of the four-point beam test: enhanced beams modeling choice.</figcaption></figure>

nonlinear behavior of concrete. In addition, considering large displacements can be done by assembling co-rotational beams or plates. Assembling higher-order layered kinematic elements can be examined, to create one enhanced multiscale element [31]. Finally, considering steel-concrete interactions in the normal direction with respect to the steel can be considered, by integrating normal stresses between steel and concrete (see [61]) in the enhanced elements.

CRediT authorship contribution statement

**Maryam Trad:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Visualization, Writing - original draft, Writing - review & editing, Project administration, Resources, Software, Supervision. **Ibrahim Bitar:** Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Conceptualization. **Stéphane Grange:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing. **Benjamin Richard:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing.
<figure><img src="image_29.png" /><figcaption>Fig. 25. Local cracking description of the four-point beam test: enhanced beams modeling choice.</figcaption></figure>

nonlinear behavior of concrete. In addition, considering large displacements can be done by assembling co-rotational beams or plates. Assembling higher-order layered kinematic elements can be examined, to create one enhanced multiscale element [31]. Finally, considering steel-concrete interactions in the normal direction with respect to the steel can be considered, by integrating normal stresses between steel and concrete (see [61]) in the enhanced elements.

CRediT authorship contribution statement

**Maryam Trad:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Visualization, Writing - original draft, Writing - review & editing, Project administration, Resources, Software, Supervision. **Ibrahim Bitar:** Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Conceptualization. **Stéphane Grange:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing. **Benjamin Richard:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing.
**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 will be made available on request.

**Acknowledgments**

The authors wish to express their most grateful thanks to IRSN, France for its financial and technical support. This work was supported and partially funded by EDF, France.

<figure><img src="image_30.png" /><figcaption>Fig. 26. Local cracking description of the four-point beam test: enhanced plates modeling choice.</figcaption></figure>
**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 will be made available on request.

**Acknowledgments**

The authors wish to express their most grateful thanks to IRSN, France for its financial and technical support. This work was supported and partially funded by EDF, France.

<figure><img src="image_30.png" /><figcaption>Fig. 26. Local cracking description of the four-point beam test: enhanced plates modeling choice.</figcaption></figure>
## Appendix A. Static condensation

The degrees of freedom vector of one enhanced element $u_{el}$ is described as follows:

$$ u_{el} = [u_e \quad u_l]^T \qquad (A.1) $$

$u_e$ et $u_l$ refer to the external (boundary) and internal degrees of freedom, respectively. The linearization of Eq. (15) gives:

$$ k_{sc} du_{el} = dF_l^e \qquad (A.2) $$

$k_{sc}$ represents the stiffness matrix of the enhanced element. The vector $du_{el}$ represents the incremental values of the total degrees of freedom of the enhanced element (internal and external degrees). $dF_l^e$ is the vector of incremental values of the resisting forces of this element, where $dF_l^e = [df_e \quad df_l]^T$. $df_e$ represents the incremental values of resisting forces at the external nodes of the enhanced element. On the other hand, $df_l$ represents the incremental values of resisting forces at the internal nodes. At convergence, the static condensation imposes that $df_l$ is equal to zero (internal equilibrium). Thus, Eq. (A.2) is expressed in an expanded form as follows:

$$ \begin{bmatrix} k_{ee} & k_{el} \\ k_{le} & k_{ll} \end{bmatrix} \begin{bmatrix} du_e \\ du_l \end{bmatrix} = \begin{bmatrix} df_e \\ 0 \end{bmatrix} \qquad (A.3) $$

Consequently, static condensation provides a link between the incremental vectors $du_e$ and $df_e$ such that:

$$ (k_{ee} - k_{el} k_{ll}^{-1} k_{le}) du_e = df_e \qquad (A.4) $$

In a more condensed form, Eq. (A.4) is expressed as follows

$$ k_{en} du_e = df_e \qquad (A.5) $$

$k_{en}$ is the condensed stiffness matrix of the enhanced element.

Enhanced beam and plate elements that require static condensation involve the use of a local Newton-Raphson algorithm (algorithm 1) at the level of each enhanced element, in addition to the Newton-Raphson standard global algorithm used to solve the overall reinforced concrete structure equilibrium. $conv_i$ used in algorithm 1 is as a convergence indicator, and $k_i$ is a local iterations counter.

A sub-structured resolution is performed where the local algorithm solves the inner equilibrium equation of one enhanced element ($f_r = 0$). It is called by the global algorithm whenever the internal forces vector or the rigidity matrix of one enhanced element is calculated.

## Algorithm 1 Local resolution algorithm of one enhanced element

$conv_i \leftarrow 0$

$k_i \leftarrow 0$

$u_{el}^T = [u_e \quad u_l]^T$. $u_e$ is provided by the global algorithm. It remains constant in this local one. An initial estimate of $u_l$ is its value at the last converged time step.

while $conv_i = 0$ do

Compute $F_l^e(u_{el})$ and $k_{sc}$. $k_{sc}$ is the rigidity matrix of the enhanced element equal to $\frac{dF_e}{du_{el}}$.

Compute $f_l = F_l^e$ (the part of $F_l^e$ at the internal nodes).

Compute $R_{k_i} = -f_i$.

if $R_{k_i} < tolerance$ then

$conv_i \leftarrow 1$

Static condensation: $k_{en} = k_{ee} - k_{el} k_{ll}^{-1} k_{le}$ (see Eq. (A.4)) and $f_e = F_e^e$ (the part of $F_e^e$ at the external nodes). $k_{en}$ and $f_e$ represent the condensed stiffness matrix and the internal forces vector of the enhanced element to be returned to the global algorithm.

else

$u_{i_{k_i+1}} = u_{k_i} + k_{ll}^{-1} R_{k_i}$

$k_i \leftarrow k_i + 1$

## Appendix B. Enhanced element formulation: two concrete beams and one steel bar

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as

$u_{el} = [u_1 \quad v_1 \quad \theta_1 \quad u_{s1} \quad u_2 \quad v_2 \quad \theta_2 \quad u_{s2} \quad u_i \quad v_i \quad \theta_i \quad u_{si}]^T$, such that:

$$ u_{el} = [u_e \quad u_l]^T \qquad (B.1) $$

with:

$$ u_e = [u_1 \quad v_1 \quad \theta_1 \quad u_{s1} \quad u_2 \quad v_2 \quad \theta_2 \quad u_{s2}]^T, u_l = [u_i \quad v_i \quad \theta_i \quad u_{si}]^T \qquad (B.2) $$
## Appendix A. Static condensation

The degrees of freedom vector of one enhanced element $u_{el}$ is described as follows:

$$ u_{el} = [u_e \quad u_l]^T \qquad (A.1) $$

$u_e$ et $u_l$ refer to the external (boundary) and internal degrees of freedom, respectively. The linearization of Eq. (15) gives:

$$ k_{sc} du_{el} = dF_{i}^{e} \qquad (A.2) $$

$k_{sc}$ represents the stiffness matrix of the enhanced element. The vector $du_{el}$ represents the incremental values of the total degrees of freedom of the enhanced element (internal and external degrees). $dF_{i}^{e}$ is the vector of incremental values of the resisting forces of this element, where $dF_{i}^{eT} = [df_{e} \ df_{l}]^{T}$. $df_{e}$ represents the incremental values of resisting forces at the external nodes of the enhanced element. On the other hand, $df_{l}$ represents the incremental values of resisting forces at the internal nodes. At convergence, the static condensation imposes that $df_{l}$ is equal to zero (internal equilibrium). Thus, Eq. (A.2) is expressed in an expanded form as follows:

$$ \begin{bmatrix} k_{ee} & k_{el} \\ k_{le} & k_{ll} \end{bmatrix} \begin{bmatrix} du_e \\ du_l \end{bmatrix} = \begin{bmatrix} df_e \\ 0 \end{bmatrix} \qquad (A.3) $$

Consequently, static condensation provides a link between the incremental vectors $du_e$ and $df_e$ such that:

$$ (k_{ee} - k_{el} k_{ll}^{-1} k_{le}) du_e = df_e \qquad (A.4) $$

In a more condensed form, Eq. (A.4) is expressed as follows

$$ k_{en} du_e = df_e \qquad (A.5) $$

$k_{en}$ is the condensed stiffness matrix of the enhanced element.

Enhanced beam and plate elements that require static condensation involve the use of a local Newton-Raphson algorithm (algorithm 1) at the level of each enhanced element, in addition to the Newton-Raphson standard global algorithm used to solve the overall reinforced concrete structure equilibrium. $conv_i$ used in algorithm 1 is as a convergence indicator, and $k_i$ is a local iterations counter.

A sub-structured resolution is performed where the local algorithm solves the inner equilibrium equation of one enhanced element ($f_r = 0$). It is called by the global algorithm whenever the internal forces vector or the rigidity matrix of one enhanced element is calculated.

## Algorithm 1 Local resolution algorithm of one enhanced element

$conv_i \leftarrow 0$

$k_i \leftarrow 0$

$u_{el}^T = [u_e \quad u_l]^T$. $u_e$ is provided by the global algorithm. It remains constant in this local one. An initial estimate of $u_l$ is its value at the last converged time step.

while $conv_i = 0$ do

Compute $F_i^e(u_{el})$ and $k_{sc}$. $k_{sc}$ is the rigidity matrix of the enhanced element equal to $\frac{dF_i^e}{du_{el}}$.

Compute $f_i = F_i^e$ (the part of $F_i^e$ at the internal nodes).

Compute $R_{k_i} = -f_i$.

if $R_{k_i} < tolerance$ then

$conv_i \leftarrow 1$

Static condensation: $k_{en} = k_{ee} - k_{el} k_{ll}^{-1} k_{le}$ (see Eq. (A.4)) and $f_e = F_e^e$ (the part of $F_e^e$ at the external nodes). $k_{en}$ and $f_e$ represent the condensed stiffness matrix and the internal forces vector of the enhanced element to be returned to the global algorithm.

else

$$ \begin{gathered} u_{i(k_i+1)} = u_{k_i} + k_{ll}^{-1} R_{k_i} \\ k_i \leftarrow k_i + 1 \end{gathered} $$

## Appendix B. Enhanced element formulation: two concrete beams and one steel bar

The elementary degrees of freedom vector $u_{el}$ of the enhanced element is defined as $u_{el} = [u_1 \ v_1 \ \theta_1 \ u_{s1} \ u_2 \ v_2 \ \theta_2 \ u_{s2} \ u_i \ v_i \ \theta_i \ u_{si}]^T$, such that:

$$ u_{el} = [u_e \quad u_l]^T \qquad (B.1) $$

with:

$$ u_e = [u_1 \ v_1 \ \theta_1 \ u_{s1} \ u_2 \ v_2 \ \theta_2 \ u_{s2}]^T, \quad u_l = [u_i \ v_i \ \theta_i \ u_{si}]^T \qquad (B.2) $$
$u_e$ and $u_l$ are the external and the internal degrees of freedom, respectively.

The steel and concrete degrees of freedom and strain vectors are assessed similarly to the previous subsection, utilizing the shape functions of the steel and concrete elements and their derivatives.

### B.0.1. Steel-concrete slip evaluation

Eq. (25) can be here used to evaluate the slip, where:

$$u_s(x) = N_b u_s \qquad (B.3)$$

knowing that $u_s = I_s u_{el}$, which gives:

$$u_s(x) = N_b I_s u_{el} = N_s u_{el} \qquad (B.4)$$

where:

$$I_s = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}, u_s = \begin{bmatrix} u_{a1} \\ u_{al} \\ u_{a2} \end{bmatrix} \qquad (B.5)$$

$N_b$ is the shape functions matrix of the three-node steel bar element. Two methods are here proposed for the calculation of $u_c(x)$. The first one is as follows:

$$u_b(x) = \underbrace{\begin{bmatrix} 1 & -y \end{bmatrix}}_{a_s} \begin{bmatrix} u(x) \\ \theta(x) \end{bmatrix} \qquad (B.6)$$

$u(x)$ and $\theta(x)$ are the concrete beam longitudinal displacement and section rotation. For a longitudinal position $0 \le x < l/2$, the shape functions of concrete beam 1 are used to evaluate $u(x)$ and $\theta(x)$. The shape functions of beam 2 are otherwise used (see Fig. 7(a)).

This first method involves a node-to-beam association for slip calculation. Hence, a second method, which bypasses this identification step, is proposed. It consists of calculating $u_f$, which is the vector of concrete longitudinal displacement near the steel bar at the beam nodes:

$$u_f = \begin{bmatrix} 1 & -y & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -y & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -y \end{bmatrix} \begin{bmatrix} u_1 \\ \theta_1 \\ u_l \\ \theta_l \\ u_2 \\ \theta_2 \end{bmatrix} \qquad (B.7)$$

so:

$$u_f = \underbrace{\begin{bmatrix} 1 & 0 & -y & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & -y & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & -y & 0 & 0 & 0 & 0 \end{bmatrix}}_{I_f} u_{el} \qquad (B.8)$$

Then the concrete longitudinal displacement at the longitudinal position $x$ is interpolated. The shape functions of a three-node bar element considered for the steel bar are used:

$$u_b(x) = N_b(x)u_f = N_b(x)I_f u_{el} \qquad (B.9)$$

Finally:

$$u_g(x) = N_s u_{el} - N_b I_f u_{el} = (N_s - N_b I_f) u_{el} \qquad (B.10)$$

### B.0.2. Weak formulation

The virtual power principle is expressed as follows:

$$u_{el}^* T \underbrace{\mathbf{A}}_{e=1}^2 \mathbf{F}_{intc}^e + \int_0^{L_s} \varepsilon_s^* \sigma_s (\varepsilon_s) S_s dx + \int_0^{L_s} u_g^* \tau (u_g) Pdx = \underbrace{u_{el}^* T \mathbf{F}_e}_{I_f} \qquad (B.11)$$

which gives:

$$\underbrace{\mathbf{A}}_{e=1}^2 \mathbf{F}_{intc}^e + \int_0^{L_s} \mathbf{B}_s^T \sigma_s (\varepsilon_s) S_s dx + \int_0^{L_s} (\mathbf{N}_s^T - \mathbf{I}_f^T \mathbf{N}_b^T) \tau (u_g) Pdx = \mathbf{F}_e \qquad (B.12)$$

$\mathbf{B}_s$ matrix is deduced by deriving the interpolation matrix $N_s$ with respect to the enhanced element elementary degrees of freedom vector $u_{el}$. $\mathbf{F}_{intc}^e$ is an elementary internal forces vector of one of the two assembled concrete beams.
$u_e$ and $u_l$ are the external and the internal degrees of freedom, respectively.

The steel and concrete degrees of freedom and strain vectors are assessed similarly to the previous subsection, utilizing the shape functions of the steel and concrete elements and their derivatives.

### B.0.1. Steel-concrete slip evaluation

Eq. (25) can be here used to evaluate the slip, where:

$$u_s(x) = N_b u_s \qquad (B.3)$$

knowing that $u_s = I_s u_{el}$, which gives:

$$u_s(x) = N_b I_s u_{el} = N_s u_{el} \qquad (B.4)$$

where:

$$I_s = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}, u_s = \begin{bmatrix} u_{a1} \\ u_{al} \\ u_{a2} \end{bmatrix} \qquad (B.5)$$

$N_b$ is the shape functions matrix of the three-node steel bar element. Two methods are here proposed for the calculation of $u_c(x)$. The first one is as follows:

$$u_b(x) = \underbrace{\begin{bmatrix} 1 & -y \end{bmatrix}}_{a_s} \begin{bmatrix} u(x) \\ \theta(x) \end{bmatrix} \qquad (B.6)$$

$u(x)$ and $\theta(x)$ are the concrete beam longitudinal displacement and section rotation. For a longitudinal position $0 \le x < l/2$, the shape functions of concrete beam 1 are used to evaluate $u(x)$ and $\theta(x)$. The shape functions of beam 2 are otherwise used (see Fig. 7(a)).

This first method involves a node-to-beam association for slip calculation. Hence, a second method, which bypasses this identification step, is proposed. It consists of calculating $u_f$, which is the vector of concrete longitudinal displacement near the steel bar at the beam nodes:

$$u_f = \begin{bmatrix} 1 & -y & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -y & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -y \end{bmatrix} \begin{bmatrix} u_1 \\ \theta_1 \\ u_l \\ \theta_l \\ u_2 \\ \theta_2 \end{bmatrix} \qquad (B.7)$$

so:

$$u_f = \underbrace{\begin{bmatrix} 1 & 0 & -y & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & -y & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & -y & 0 & 0 & 0 & 0 \end{bmatrix}}_{I_f} u_{el} \qquad (B.8)$$

Then the concrete longitudinal displacement at the longitudinal position $x$ is interpolated. The shape functions of a three-node bar element considered for the steel bar are used:

$$u_b(x) = N_b(x)u_f = N_b(x)I_f u_{el} \qquad (B.9)$$

Finally:

$$u_g(x) = N_s u_{el} - N_b I_f u_{el} = (N_s - N_b I_f) u_{el} \qquad (B.10)$$

### B.0.2. Weak formulation

The virtual power principle is expressed as follows:

$$u_{el}^* T \underbrace{\mathbf{A}}_{e=1}^2 \mathbf{F}_{intc}^e + \int_0^{L_s} \varepsilon_s^* \sigma_s (\varepsilon_s) S_s dx + \int_0^{L_s} u_g^* \tau (u_g) P dx = \underbrace{u_{el}^* T \mathbf{F}_e}_{I_f} \qquad (B.11)$$

which gives:

$$\underbrace{\mathbf{A}}_{e=1}^2 \mathbf{F}_{intc}^e + \int_0^{L_s} \mathbf{B}_s^T \sigma_s (\varepsilon_s) S_s dx + \int_0^{L_s} (\mathbf{N}_s^T - \mathbf{I}_f^T \mathbf{N}_b^T) \tau (u_g) P dx = \mathbf{F}_e \qquad (B.12)$$

$\mathbf{B}_s$ matrix is deduced by deriving the interpolation matrix $N_s$ with respect to the enhanced element elementary degrees of freedom vector $u_{el}$. $\mathbf{F}_{intc}^e$ is an elementary internal forces vector of one of the two assembled concrete beams.
The entirety of the expression on the left-hand side of Eq. (B.12) denotes for the internal forces vector $F_i^e$ vector of the enhanced beams element. Its derivation with respect to $u_{el}$ gives the tangent operator of the enhanced element:

$$ \frac{\partial F_i^e}{\partial u_{el}} = k_{sc}^e = k_c^e + k_s^e + k_g^e \qquad (B.13) $$

$k_c^e$ is the assembly of the two concrete beams tangent operators. $k_s^e$ is the steel three-node bar element tangent operator. $k_g^e$ is the bond tangent operator evaluated as follows:

$$ k_g^e = \int_0^l (N_s^T - I_f^T N_b^T) \frac{\partial \tau(u_g)}{\partial u_{el}} P dx \qquad (B.14) $$

with:

$$ \frac{\partial \tau(u_g)}{\partial u_{el}} = \underbrace{\frac{\partial \tau(u_g)}{\partial u_g}}_{k_{int}} \frac{\partial u_g}{\partial u_{el}} = k_{int}(N_s - N_b I_f) \qquad (B.15) $$

where $k_{int}$ is the bond law slope. It is to be noted here that the second slip interpolation method of paragraph Appendix B.0.1 is chosen in Eq. (B.14). Choosing the first method remains possible. For that, the integration of this equation would be divided into the sum of two parts: an integration for $0 \le x < \frac{l}{2}$ and another one for $\frac{l}{2} \le x \le l$.

## Appendix C. Regularized 1D damage law

The 1D damage law used in 4.1 defines a Mazars' damage criterion [62]. Once the traction strain $\epsilon_{eq}$ reaches the value of $\epsilon_{d0}$, concrete starts to damage. An exponential evolution of the damage variable $D_t$ is described as follows:

$$ D_t = 1 - \frac{\epsilon_{d0}}{\epsilon_{eq}} \exp(B_t(\epsilon_{d0} - \epsilon_{eq})) \qquad (C.1) $$

A Hillerborg regularization is applied to this law [59] [63]. With no other dissipative mechanisms, the area under the total stress-strain curve represents the volumetric cracking dissipated energy. The multiplication of this area by the width $h$ of the cracking localization area represents the fracture energy $G_f$.

$$ G_f = h \int_0^\infty \sigma d\epsilon \qquad (C.2) $$

For the simple 1D case it is possible to express $G_f$ as:

$$ G_f = h \int_0^\infty E (1 - D_t) \epsilon d\epsilon = h \int_0^\infty E \frac{\epsilon_{d0}}{\epsilon_{eq}} \exp(B_t (\epsilon_{d0} - \epsilon_{eq})) \epsilon d\epsilon = h \frac{E \epsilon_{d0}^2}{2} + h \frac{E \epsilon_{d0}}{B_t} \qquad (C.3) $$

Which leads to the regularized expression of $B_t$:

$$ B_t = \frac{h E \epsilon_{d0}}{G_f - h \left( \frac{\epsilon_{d0}^2 E}{2} \right)} \qquad (C.4) $$

$B_t$ is calculated for each finite element. One simple way to calculate $h$ is to suppose that it is equal to the finite element size.

## Appendix D. Revised and regularized Mazars' constitutive law

The revised Mazars's behavior law [58] defines a single damage variable $D$ such that:

$$ D = 1 - \frac{(1 - A)Y_0}{Y} - A \exp(-B(Y - Y_0)) \qquad (D.1) $$

$Y$ is an internal variable introduced to the initial Mazars' law [62] in order to improve its shear behavior:

$$ Y = rY_t + (1 - r)Y_c \qquad (D.2) $$

such as:

$$ r = \frac{\sum_{i=1}^{3} \langle \bar{\sigma}_i \rangle_{+}}{\sum_{i=1}^{3} |\bar{\sigma}_i|} \qquad (D.3) $$

$\langle \bar{\sigma}_i \rangle_+$ and $|\bar{\sigma}_i|$ represent respectively the positive part and the absolute value of the principal stress $\bar{\sigma}_i$. The effective stresses matrix $\bar{\sigma}$ takes into account the damaged state of the material such that:

$$ \bar{\sigma} = \frac{\sigma}{(1 - D)} \qquad (D.4) $$
The entirety of the expression on the left-hand side of Eq. (B.12) denotes for the internal forces vector $F_i^e$ vector of the enhanced beams element. Its derivation with respect to $u_{el}$ gives the tangent operator of the enhanced element:

$$ \frac{\partial F_i^e}{\partial u_{el}} = k_{sc}^e = k_c^e + k_s^e + k_g^e \qquad (B.13) $$

$k_c^e$ is the assembly of the two concrete beams tangent operators. $k_s^e$ is the steel three-node bar element tangent operator. $k_g^e$ is the bond tangent operator evaluated as follows:

$$ k_g^e = \int_0^l (N_s^T - I_f^T N_b^T) \frac{\partial \tau(u_g)}{\partial u_{el}} P dx \qquad (B.14) $$

with:

$$ \frac{\partial \tau(u_g)}{\partial u_{el}} = \underbrace{\frac{\partial \tau(u_g)}{\partial u_g}}_{k_{int}} \frac{\partial u_g}{\partial u_{el}} = k_{int}(N_s - N_b I_f) \qquad (B.15) $$

where $k_{int}$ is the bond law slope. It is to be noted here that the second slip interpolation method of paragraph Appendix B.0.1 is chosen in Eq. (B.14). Choosing the first method remains possible. For that, the integration of this equation would be divided into the sum of two parts: an integration for $0 \le x < l/2$ and another one for $l/2 \le x \le l$.

## Appendix C. Regularized 1D damage law

The 1D damage law used in 4.1 defines a Mazars' damage criterion [62]. Once the traction strain $\epsilon_{eq}$ reaches the value of $\epsilon_{d0}$, concrete starts to damage. An exponential evolution of the damage variable $D_t$ is described as follows:

$$ D_t = 1 - \frac{\epsilon_{d0}}{\epsilon_{eq}} \exp(B_t(\epsilon_{d0} - \epsilon_{eq})) \qquad (C.1) $$

A Hillerborg regularization is applied to this law [59] [63]. With no other dissipative mechanisms, the area under the total stress-strain curve represents the volumetric cracking dissipated energy. The multiplication of this area by the width $h$ of the cracking localization area represents the fracture energy $G_f$.

$$ G_f = h \int_0^\infty \sigma d\epsilon \qquad (C.2) $$

For the simple 1D case it is possible to express $G_f$ as:

$$ G_f = h \int_0^\infty E (1 - D_t) \epsilon d\epsilon = h \int_0^\infty E \frac{\epsilon_{d0}}{\epsilon_{eq}} \exp(B_t (\epsilon_{d0} - \epsilon_{eq})) \epsilon d\epsilon = h \frac{E \epsilon_{d0}^2}{2} + h \frac{E \epsilon_{d0}}{B_t} \qquad (C.3) $$

Which leads to the regularized expression of $B_t$:

$$ B_t = \frac{h E \epsilon_{d0}}{G_f - h \left( \frac{\epsilon_{d0}^2 E}{2} \right)} \qquad (C.4) $$

$B_t$ is calculated for each finite element. One simple way to calculate $h$ is to suppose that it is equal to the finite element size.

## Appendix D. Revised and regularized Mazars' constitutive law

The revised Mazars's behavior law [58] defines a single damage variable $D$ such that:

$$ D = 1 - \frac{(1 - A)Y_0}{Y} - A \exp(-B(Y - Y_0)) \qquad (D.1) $$

$Y$ is an internal variable introduced to the initial Mazars' law [62] in order to improve its shear behavior:

$$ Y = rY_t + (1 - r)Y_c \qquad (D.2) $$

such as:

$$ r = \frac{\sum_{i=1}^{3} \langle \bar{\sigma}_i \rangle_{+}}{\sum_{i=1}^{3} |\bar{\sigma}_i|} \qquad (D.3) $$

$\langle \bar{\sigma}_i \rangle_+$ and $|\bar{\sigma}_i|$ represent respectively the positive part and the absolute value of the principal stress $\bar{\sigma}_i$. The effective stresses matrix $\bar{\sigma}$ takes into account the damaged state of the material such that:

$$ \bar{\sigma} = \frac{\sigma}{(1 - D)} \qquad (D.4) $$
<figure><img src="image_31.png" /><figcaption>Fig. D.1. Principe du calcul de l'énergie de fissuration.</figcaption></figure>

where $\sigma$ is the stress matrix of the undamaged material. The variables $Y_t$ and $Y_c$ of Eq. (D.2) are associated with the loading surfaces. Two surfaces for tension and compression are defined such that:

$$f_t = \varepsilon_t - Y_t \qquad (D.5)$$

and:

$$f_c = \varepsilon_c - Y_c \qquad (D.6)$$

with:

$$\varepsilon_t = \frac{I_\varepsilon}{2(1-2\nu)} + \frac{\sqrt{J_\varepsilon}}{2(1+\nu)} \qquad (D.7)$$

and:

$$\varepsilon_c = \frac{I_\varepsilon}{5(1-2\nu)} + \frac{6\sqrt{J_\varepsilon}}{5(1+\nu)} \qquad (D.8)$$

such that:

$$I_\varepsilon = \varepsilon_1 + \varepsilon_2 + \varepsilon_3 \qquad (D.9)$$

$$J_\varepsilon = 0.5 [(\varepsilon_1 - \varepsilon_2)^2 + (\varepsilon_2 - \varepsilon_3)^2 + (\varepsilon_3 - \varepsilon_1)^2] \qquad (D.10)$$

where $\varepsilon_1$, $\varepsilon_2$, and $\varepsilon_3$ represent the principal values of the strain tensor $\varepsilon$. $Y_t$ and $Y_c$ are defined as follows:

$$Y_t = \text{Sup} [\varepsilon_{0t}, \text{max}(\varepsilon_t)] , \quad Y_c = \text{Sup} [\varepsilon_{0c}, \text{max}(\varepsilon_c)] \qquad (D.11)$$

The symbol “Sup” denotes the greater value of the considered set. $\varepsilon_{0t}$ and $\varepsilon_{0c}$ are the associated threshold values to $\varepsilon_t$ and $\varepsilon_c$, respectively.

The parameters A and B of Eq. (D.1) define the shape of the law:

$$\left\{ \begin{array}{l} A = A_t (2r^2(1-2k) - r(1-4k)) + A_c (2r^2 - 3r + 1) \\ B = r(r^2-2r+2) B_t + (1-r(r^2-2r+2)) B_c \end{array} \right. \qquad (D.12)$$

For a purely tensile behavior, A is equal to $A_t$, and B is equal to $B_t$. In pure compression, A is equal to $A_c$ and B is equal to $B_c$.

The energy regularization of the tensile behavior of this revised Mazars' law involves the identification of the regularized values of $A_t$ and $B_t$. Eq. (C.2) defines the cracking energy as the crack localization width $h$ multiplied by the area under the $\sigma - \varepsilon$ curve. In 1D, it is possible to draw a triangle with the same area for the revised regularized Mazars' law, as shown in Fig. D.1.

From Fig. D.1, one can deduce the following:

$$h = 2 \frac{G_f}{f_t^2} \left( \frac{1}{E} - \frac{1}{E_t} \right)^{-1} \qquad (D.13)$$

$h$ can be considered equal to the size of the finite element. $G_f$, $E_l$, and $f_t$ are material parameters. Therefore, the value of $E_t$ can be deduced using Eq. (D.13). The regularization process involves plotting the equivalent-area triangle of the constitutive law and determining suitable values for $A_t$ and $B_t$ of the regularized law. This is achieved by plotting the tensile behavior of the law
<figure><img src="image_31.png" /><figcaption>Fig. D.1. Principe du calcul de l'énergie de fissuration.</figcaption></figure>

where $\sigma$ is the stress matrix of the undamaged material. The variables $Y_t$ and $Y_c$ of Eq. (D.2) are associated with the loading surfaces. Two surfaces for tension and compression are defined such that:

$$f_t = \varepsilon_t - Y_t \qquad (D.5)$$

and:

$$f_c = \varepsilon_c - Y_c \qquad (D.6)$$

with:

$$\varepsilon_t = \frac{I_{\varepsilon}}{2(1 - 2\nu)} + \frac{\sqrt{J_{\varepsilon}}}{2(1 + \nu)} \qquad (D.7)$$

and:

$$\varepsilon_c = \frac{I_{\varepsilon}}{5(1 - 2\nu)} + \frac{6\sqrt{J_{\varepsilon}}}{5(1 + \nu)} \qquad (D.8)$$

such that:

$$I_{\varepsilon} = \varepsilon_{1} + \varepsilon_{2} + \varepsilon_{3} \qquad (D.9)$$

$$J_{\varepsilon} = 0.5 \left[ (\varepsilon_{1} - \varepsilon_{2})^{2} + (\varepsilon_{2} - \varepsilon_{3})^{2} + (\varepsilon_{3} - \varepsilon_{1})^{2} \right] \qquad (D.10)$$

where $\varepsilon_1$, $\varepsilon_2$, and $\varepsilon_3$ represent the principal values of the strain tensor $\varepsilon$. $Y_t$ and $Y_c$ are defined as follows:

$$Y_t = \mathrm{Sup} [\varepsilon_{0t}, \max(\varepsilon_t)] , \quad Y_c = \mathrm{Sup} [\varepsilon_{0c}, \max(\varepsilon_c)] \qquad (D.11)$$

The symbol “Sup” denotes the greater value of the considered set. $\varepsilon_{0t}$ and $\varepsilon_{0c}$ are the associated threshold values to $\varepsilon_t$ and $\varepsilon_c$, respectively.

The parameters A and B of Eq. (D.1) define the shape of the law:

$$\left\{ \begin{array}{l} A = A_t (2r^2(1-2k) - r(1-4k)) + A_c (2r^2 - 3r + 1) \\ B = r(r^2-2r+2) B_t + (1-r(r^2-2r+2)) B_c \end{array} \right. \qquad (D.12)$$

For a purely tensile behavior, A is equal to $A_t$, and B is equal to $B_t$. In pure compression, A is equal to $A_c$ and B is equal to $B_c$.

The energy regularization of the tensile behavior of this revised Mazars' law involves the identification of the regularized values of $A_t$ and $B_t$. Eq. (C.2) defines the cracking energy as the crack localization width $h$ multiplied by the area under the $\sigma - \varepsilon$ curve. In 1D, it is possible to draw a triangle with the same area for the revised regularized Mazars' law, as shown in Fig. D.1.

From Fig. D.1, one can deduce the following:

$$h = 2 \frac{G_f}{f_t^2} \left( \frac{1}{E} - \frac{1}{E_t} \right)^{-1} \qquad (D.13)$$

$h$ can be considered equal to the size of the finite element. $G_f$, $E_l$, and $f_t$ are material parameters. Therefore, the value of $E_t$ can be deduced using Eq. (D.13). The regularization process involves plotting the equivalent-area triangle of the constitutive law and determining suitable values for $A_t$ and $B_t$ of the regularized law. This is achieved by plotting the tensile behavior of the law
for each pair of values ($A_t$, $B_t$) and calculating the area under the law curve. The objective is to find the values of $A_t$ and $B_t$ that give the closest area to the area of the initially plotted triangle.

The value of $A_t$ defines the asymptote of the law of large strain values. $A_t = 1$ is associated corresponds to an asymptote of a stress equal to zero. This value of $A_t$ gives the best representation of the experimental concrete behavior, but it may lead to numerical convergence problems. In this case, a lower value is used. The value of $A_t$ is between 0 and 1.

## References

[1] N. Herrmann, H.S. Müller, S. Michel-Ponnelle, B. Masson, M. Herve, The PACE-1450 experiment-investigations regarding crack and leakage behaviour of a pre-stressed concrete containment, in: High Tech Concrete: Where Technology and Engineering Meet: Proceedings of the 2017 Fib Symposium, Held in Maastricht, The Netherlands, 12-14, Springer International Publishing, 2017, pp. 1487-1495.

[2] D. Ngo, A.C. Scordelis, Finite element analysis of reinforced concrete beams, ACI J. 6 (3) (1967) 152-163, http://dx.doi.org/10.14359/7551.

[3] A. Lin, C.P. Ostertag, Multiscale pull-out resistance of steel reinforcing bar embedded in hybrid fiber reinforced concrete (HYFRC), IOP Conf. Ser.: Mater. Sci. Eng. 246 (1) (2017) 012022, http://dx.doi.org/10.1088/1757-899X/246/1/012022.

[4] R. Tepfers, A Theory of Bond Applied to Overlapped Tensile Reinforcement Splices for Deformed Bars, Chalmers University of Technology, Publication 73:2, Göteborg, Sweden, 1973, p. 328.

[5] M.T.G. Barbosa, EDS. Sánchez Filho, T.M.D. Oliveira, WJD. Santos, Analysis of the relative rib area of reinforcing bars pull out tests, Mater. Res. 11 (2008) 453-457, http://dx.doi.org/10.1590/S1516-14392008000400013.

[6] A. Daoud, M. Lorrain, Influence de la position des armatures sur l'adhérence des bétons autoplacants: interprétation par analyse d'image, Mater. Struct. Mat. Constr. 36 (2003) 231-237, http://dx.doi.org/10.1007/BF02479616.

[7] B.H. Tran, Y. Berthaud, F. Ragueneau, Essais PIAF: Pour identifier l'adhérence et le frottement, in: CFM 2007-18ème Congrès Français de Mécanique, 2007.

[8] A. Tixier, C. Rospars, F. Dufour, A. Khadour, M. Quiertant, B. Masson, Numerical modeling to analyse optical fiber measurements along a steel-concrete interface, in: Computational Modeling of Fracture and Failure of Materials and Structures Conference, 2013.

[9] ACI Committee 408, A guide for determination of bond strength in beam specimens, J. Proc. 61 (2) (1964) 129-136, http://dx.doi.org/10.14359/7769.

[10] RILEM, Essai portant sur l'adhérence des armatures du béton I. essai par flexion, Mat. Constr. 3 (15) (1970) 169-174, http://dx.doi.org/10.1007/BF02478967.

[11] A.H. Nilson, Nonlinear analysis of reinforced concrete by the finite element method, J. Proc. 65 (9) (1968) 757-766, http://dx.doi.org/10.14359/7510.

[12] H.W. Reinhardt, J. Blaauwendraad, E. Vos, Prediction of bond between steel and concrete by numerical analysis, Mat. Constr. 17 (4) (1984) 311-320, http://dx.doi.org/10.1007/BF02479089.

[13] B. Richard, F. Ragueneau, C. Cremona, L. Adelaide, J.L. Tailhana, A three-dimensional steel/concrete interface model including corrosion effects, Eng. Fract. Mech. 77 (6) (2010) 951-973, http://dx.doi.org/10.1016/j.engfracmech.2010.01.017.

[14] M. Abbas, B. Bary, L. Jason, A 3D mesoscopic frictional cohesive zone model for the steel-concrete interface, Int. J. Mech. Sci. (2022) 107819, http://dx.doi.org/10.1016/j.ijmecsci.2022.107819.

[15] C. Mang, L. Jason, L. Davenne, A new bond slip model for reinforced concrete structures in monotonic and cyclic loadings, Eng. Comput. 32 (7) (2015) 1934-1958, http://dx.doi.org/10.1108/EC-11-2014-0234.

[16] C. Turgut, L. Jason, L. Davenne, Structural-scale modeling of the active confinement effect in the steel-concrete bond for reinforced concrete structures, Finite Elem. Anal. Des. 172 (2020) 103386, http://dx.doi.org/10.1016/j.finel.2020.103386.

[17] A. Sellier, A. Millard, A homogenized formulation to account for sliding of non-meshed reinforcements during the cracking of brittle matrix composites: Application to reinforced concrete, Eng. Fract. Mech. 213 (2019) 182-196, http://dx.doi.org/10.1016/j.engfracmech.2019.04.008.

[18] F. Wang, S. Lavarenne, T. Chaudat, D. Combescure, T. Payen, B. Fouré, Simulation analysis of static and shaking table tests on rc columns with insufficient lap splices, in: SMIRT 19th Conference, 2007.

[19] C. Christelle, H. Dumontet, F. Voldoire, Homogenised constitutive model coupling damage and debonding for reinforced concrete structures under cyclic solicitations, Int. J. Solids Struct. 50 (24) (2013) 3861-3874, http://dx.doi.org/10.1016/j.ijsolstr.2013.07.021.

[20] B. Richard, F. Ragueneau, L. Adelaïde, C. Cremona, A multi-fiber approach for modeling corroded reinforced concrete structures, Eur. J. Mech. A Solids 30 (6) (2011) 950-961, http://dx.doi.org/10.1016/j.euromechsol.2011.06.002.

[21] C. Combescure, Formulation d'un Modèle Homogénéisé de Plaque en Béton Armé Pour des Applications Sismiques, Université Pierre et Marie Curie, 2013.

[22] M.H. Aguilera, Un Modèle Global Homogénéisé Pour la Fissuration des Voiles en Béton Armé Sous Chargements Sismiques, École Centrale de Nantes, 2016.

[23] R. Benjamin, R. Frédéric, A. Lucas, C. Christian, A multi-fiber approach for modeling corroded reinforced concrete structures, Eur. J. Mech. A Solids 30 (6) (2011) 950-961, http://dx.doi.org/10.1016/j.euromechsol.2011.06.002.

[24] S. Limkatanyu, E. Spacone, Reinforced concrete frame element with bond interfaces. I: Displacement-based, force-based, and mixed formulations, J. Struct. Eng. 128 (3) (2002) 346-355, http://dx.doi.org/10.1061/(ASCE)0733-9445(2002)128:3(346).

[25] B. Yousefi, M.R. Esfahani, M. Tavakkolizadehm, A multi-fiber approach with directional stiffness matrix in reinforced concrete structures, Eng. Comput. 37 (7) (2020) 2411-2437, http://dx.doi.org/10.1108/EC-09-2019-0424.

[26] S. Abtahi, Y. Li, Efficient modeling of steel bar slippage effect in reinforced concrete structures using a newly implemented nonlinear element, Comput. Struct. 279 (2023) 106958, http://dx.doi.org/10.1016/j.comstruc.2022.106958.

[27] J.C. Yu, J.T. Wang, J.W. Pan, N. Guo, C.H. Zhang, A dynamic FEM-DEM multiscale modeling approach for concrete structures, Eng. Fract. Mech. 278 (2023) 109031, http://dx.doi.org/10.1016/j.engfracmech.2022.109031.

[28] M. Petrolo, E. Carrera, M. Cinefra, E. Zappino, Finite Element Analysis of Structures Through Unified Formulation, John Wiley & Sons, 2014.

[29] J. Shen, A. Pagani, M.R.T. Arruda, E. Carrera, Exact component-wise solutions for 3D free vibration and stress analysis of hybrid steel-concrete composite beams, Thin-Walled Struct. 174 (2022) 109094, http://dx.doi.org/10.1016/j.tws.2022.109094.

[30] R. Wang, Z. Fang, M. Lezgy-Nazargah, H. Khosravi, Nonlinear analysis of reinforced concrete slabs using a quasi-3D mixed finite element formulation, Eng. Struct. 294 (2023) 116781, http://dx.doi.org/10.1016/j.engstruct.2023.116781.

[31] J. Shen, MT. Arruda, A. Pagani, Concrete damage analysis based on higher-order beam theories using fracture energy regularization, Mech. Adv. Mater. Struct. (2022) 1-15, http://dx.doi.org/10.1080/15376494.2022.2098430.

[32] 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.

[33] D.R.J. Owen, E. Hinton, Finite Elements Plasticity: Theory and Practice, Pineridge Press Limited, 1980.

[34] P. Kotronis, F. Ragueneau, J. Mazars, A simplified modelling strategy for R/C walls satisfying PS92 and EC8 design, Eng. Struct. 27 (8) (2005) 1197-1208, http://dx.doi.org/10.1016/j.engstruct.2005.03.003.

[35] S. Grange, L. Botrugno, P. Kotronis, C. Tamagnini, The effects of soil-structure interaction on a reinforced concrete viaduct, Earthq. Eng. Struct. Dyn. 41 (11) (2011) 1549-1568.
for each pair of values ($A_t$, $B_t$) and calculating the area under the law curve. The objective is to find the values of $A_t$ and $B_t$ that give the closest area to the area of the initially plotted triangle.

The value of $A_t$ defines the asymptote of the law of large strain values. $A_t = 1$ is associated corresponds to an asymptote of a stress equal to zero. This value of $A_t$ gives the best representation of the experimental concrete behavior, but it may lead to numerical convergence problems. In this case, a lower value is used. The value of $A_t$ is between 0 and 1.

## References

[1] N. Herrmann, H.S. Müller, S. Michel-Ponnelle, B. Masson, M. Herve, The PACE-1450 experiment-investigations regarding crack and leakage behaviour of a pre-stressed concrete containment, in: High Tech Concrete: Where Technology and Engineering Meet: Proceedings of the 2017 Fib Symposium, Held in Maastricht, The Netherlands, 12-14, Springer International Publishing, 2017, pp. 1487-1495.

[2] D. Ngo, A.C. Scordelis, Finite element analysis of reinforced concrete beams, ACI J. 6 (3) (1967) 152-163, http://dx.doi.org/10.14359/7551.

[3] A. Lin, C.P. Ostertag, Multiscale pull-out resistance of steel reinforcing bar embedded in hybrid fiber reinforced concrete (HYFRC), IOP Conf. Ser.: Mater. Sci. Eng. 246 (1) (2017) 012022, http://dx.doi.org/10.1088/1757-899X/246/1/012022.

[4] R. Tepfers, A Theory of Bond Applied to Overlapped Tensile Reinforcement Splices for Deformed Bars, Chalmers University of Technology, Publication 73:2, Göteborg, Sweden, 1973, p. 328.

[5] M.T.G. Barbosa, EDS. Sánchez Filho, T.M.D. Oliveira, WJD. Santos, Analysis of the relative rib area of reinforcing bars pull out tests, Mater. Res. 11 (2008) 453-457, http://dx.doi.org/10.1590/S1516-14392008000400013.

[6] A. Daoud, M. Lorrain, Influence de la position des armatures sur l'adhérence des bétons autoplacants: interprétation par analyse d'image, Mater. Struct. Mat. Constr. 36 (2003) 231-237, http://dx.doi.org/10.1007/BF02479616.

[7] B.H. Tran, Y. Berthaud, F. Ragueneau, Essais PIAF: Pour identifier l'adhérence et le frottement, in: CFM 2007-18ème Congrès Français de Mécanique, 2007.

[8] A. Tixier, C. Rospars, F. Dufour, A. Khadour, M. Quiertant, B. Masson, Numerical modeling to analyse optical fiber measurements along a steel-concrete interface, in: Computational Modeling of Fracture and Failure of Materials and Structures Conference, 2013.

[9] ACI Committee 408, A guide for determination of bond strength in beam specimens, J. Proc. 61 (2) (1964) 129-136, http://dx.doi.org/10.14359/7769.

[10] RILEM, Essai portant sur l'adhérence des armatures du béton I. essai par flexion, Mat. Constr. 3 (15) (1970) 169-174, http://dx.doi.org/10.1007/BF02478967.

[11] A.H. Nilson, Nonlinear analysis of reinforced concrete by the finite element method, J. Proc. 65 (9) (1968) 757-766, http://dx.doi.org/10.14359/7510.

[12] H.W. Reinhardt, J. Blaauwendraad, E. Vos, Prediction of bond between steel and concrete by numerical analysis, Mat. Constr. 17 (4) (1984) 311-320, http://dx.doi.org/10.1007/BF02479089.

[13] B. Richard, F. Ragueneau, C. Cremona, L. Adelaide, J.L. Tailhana, A three-dimensional steel/concrete interface model including corrosion effects, Eng. Fract. Mech. 77 (6) (2010) 951-973, http://dx.doi.org/10.1016/j.engfracmech.2010.01.017.

[14] M. Abbas, B. Bary, L. Jason, A 3D mesoscopic frictional cohesive zone model for the steel-concrete interface, Int. J. Mech. Sci. (2022) 107819, http://dx.doi.org/10.1016/j.ijmecsci.2022.107819.

[15] C. Mang, L. Jason, L. Davenne, A new bond slip model for reinforced concrete structures in monotonic and cyclic loadings, Eng. Comput. 32 (7) (2015) 1934-1958, http://dx.doi.org/10.1108/EC-11-2014-0234.

[16] C. Turgut, L. Jason, L. Davenne, Structural-scale modeling of the active confinement effect in the steel-concrete bond for reinforced concrete structures, Finite Elem. Anal. Des. 172 (2020) 103386, http://dx.doi.org/10.1016/j.finel.2020.103386.

[17] A. Sellier, A. Millard, A homogenized formulation to account for sliding of non-meshed reinforcements during the cracking of brittle matrix composites: Application to reinforced concrete, Eng. Fract. Mech. 213 (2019) 182-196, http://dx.doi.org/10.1016/j.engfracmech.2019.04.008.

[18] F. Wang, S. Lavarenne, T. Chaudat, D. Combescure, T. Payen, B. Fouré, Simulation analysis of static and shaking table tests on rc columns with insufficient lap splices, in: SMIRT 19th Conference, 2007.

[19] C. Christelle, H. Dumontet, F. Voldoire, Homogenised constitutive model coupling damage and debonding for reinforced concrete structures under cyclic solicitations, Int. J. Solids Struct. 50 (24) (2013) 3861-3874, http://dx.doi.org/10.1016/j.ijsolstr.2013.07.021.

[20] B. Richard, F. Ragueneau, L. Adelaïde, C. Cremona, A multi-fiber approach for modeling corroded reinforced concrete structures, Eur. J. Mech. A Solids 30 (6) (2011) 950-961, http://dx.doi.org/10.1016/j.euromechsol.2011.06.002.

[21] C. Combescure, Formulation d'un Modèle Homogénéisé de Plaque en Béton Armé Pour des Applications Sismiques, Université Pierre et Marie Curie, 2013.

[22] M.H. Aguilera, Un Modèle Global Homogénéisé Pour la Fissuration des Voiles en Béton Armé Sous Chargements Sismiques, École Centrale de Nantes, 2016.

[23] R. Benjamin, R. Frédéric, A. Lucas, C. Christian, A multi-fiber approach for modeling corroded reinforced concrete structures, Eur. J. Mech. A Solids 30 (6) (2011) 950-961, http://dx.doi.org/10.1016/j.euromechsol.2011.06.002.

[24] S. Limkatanyu, E. Spacone, Reinforced concrete frame element with bond interfaces. I: Displacement-based, force-based, and mixed formulations, J. Struct. Eng. 128 (3) (2002) 346-355, http://dx.doi.org/10.1061/(ASCE)0733-9445(2002)128:3(346).

[25] B. Yousefi, M.R. Esfahani, M. Tavakkolizadehm, A multi-fiber approach with directional stiffness matrix in reinforced concrete structures, Eng. Comput. 37 (7) (2020) 2411-2437, http://dx.doi.org/10.1108/EC-09-2019-0424.

[26] S. Abtahi, Y. Li, Efficient modeling of steel bar slippage effect in reinforced concrete structures using a newly implemented nonlinear element, Comput. Struct. 279 (2023) 106958, http://dx.doi.org/10.1016/j.comstruc.2022.106958.

[27] J.C. Yu, J.T. Wang, J.W. Pan, N. Guo, C.H. Zhang, A dynamic FEM-DEM multiscale modeling approach for concrete structures, Eng. Fract. Mech. 278 (2023) 109031, http://dx.doi.org/10.1016/j.engfracmech.2022.109031.

[28] M. Petrolo, E. Carrera, M. Cinefra, E. Zappino, Finite Element Analysis of Structures Through Unified Formulation, John Wiley & Sons, 2014.

[29] J. Shen, A. Pagani, M.R.T. Arruda, E. Carrera, Exact component-wise solutions for 3D free vibration and stress analysis of hybrid steel-concrete composite beams, Thin-Walled Struct. 174 (2022) 109094, http://dx.doi.org/10.1016/j.tws.2022.109094.

[30] R. Wang, Z. Fang, M. Lezgy-Nazargah, H. Khosravi, Nonlinear analysis of reinforced concrete slabs using a quasi-3D mixed finite element formulation, Eng. Struct. 294 (2023) 116781, http://dx.doi.org/10.1016/j.engstruct.2023.116781.

[31] J. Shen, MT. Arruda, A. Pagani, Concrete damage analysis based on higher-order beam theories using fracture energy regularization, Mech. Adv. Mater. Struct. (2022) 1-15, http://dx.doi.org/10.1080/15376494.2022.2098430.

[32] 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.

[33] D.R.J. Owen, E. Hinton, Finite Elements Plasticity: Theory and Practice, Pineridge Press Limited, 1980.

[34] P. Kotronis, F. Ragueneau, J. Mazars, A simplified modelling strategy for R/C walls satisfying PS92 and EC8 design, Eng. Struct. 27 (8) (2005) 1197-1208, http://dx.doi.org/10.1016/j.engstruct.2005.03.003.

[35] S. Grange, L. Botrugno, P. Kotronis, C. Tamagnini, The effects of soil-structure interaction on a reinforced concrete viaduct, Earthq. Eng. Struct. Dyn. 41 (11) (2011) 1549-1568.
[36] I. Bitar, N. Benkemoun, P. Kotronis, S. Grange, A multifiber Timoshenko beam with embedded discontinuities, Eng. Fract. Mech. 214 (2019) 339–364, http://dx.doi.org/10.1016/j.engfracmech.2019.03.032.

[37] I. Bitar, P. Kotronis, N. Benkemoun, S. Grange, A generalized Timoshenko beam with embedded rotation discontinuity, Finite Elem. Anal. Des. 150 (2018) 34–50, http://dx.doi.org/10.1016/j.finel.2018.07.002.

[38] C. Lejouad, S. Capdevielle, B. Richard, F. Ragueneau, A numerical model to predict the effect of corrosion on the dynamic behavior of reinforced concrete beams, in: 12th Canadian Conference of Earthquake Engineering, 2019.

[39] N. Khoder, S. Grange, Y. Sieffert, Enhancement of a two-dimensional multifibre beam element in the case of reinforced concrete structures for taking into account the lateral confinement of concrete due to stirrups, Eur. J. Environ. Civ. Eng. 23 (5) (2019) 564–585, http://dx.doi.org/10.1080/19648189.2018.1446364.

[40] N. Viet-Tung, C. Jean-François, S. Karam, A model for thick laminates and sandwich plates, Compos. Sci. Technol. 65 (3-4) (2005) 475–489, http://dx.doi.org/10.1016/j.compsitech.2004.09.014.

[41] P. Zoghipour, K. Hasim, A. Kefal, M. Yildiz, Implementation of shear-locking-free triangular refined zigzag element for structural analysis of multilayered plates with curvilinear fibers, Compos. Struct. 305 (2023) 116462, http://dx.doi.org/10.1016/j.compstruct.2022.116462.

[42] E. Carrera, R. Augello, A. Pagani, D. Scano, Refined multilayered beam, plate and shell elements based on Jacobi polynomials, Compos. Struct. 304 (2023) 116275, http://dx.doi.org/10.1016/j.compstruct.2022.116275.

[43] S. Grange, D. Bertrand, Implicit coupling of heterogeneous and asynchronous time-schemes using a primal approach based on velocity continuity at the subdomain interface, Finite Elem. Anal. Des. 196 (2021) 103604, http://dx.doi.org/10.1016/j.finel.2021.103604.

[44] R. Eligethausen, E.P. Popov, V.V. Bertero, Local bond stress-slip relationships of deformed bars under generalized excitations, 1982, http://dx.doi.org/10.18419/opus-415.

[45] M.H. Harajli, Development/splice strength of reinforcing bars embedded in plain and fiber reinforced concrete, ACI Mater. J. 91 (5) (1994) 511–520, http://dx.doi.org/10.14359/4163.

[46] P. Desnerck, G. De Schutter, L. Taerwe, A local bond stress-slip model for reinforcing bars in self-compacting concrete, fracture mechanics of concrete and concrete strucures-assesment, in: Durability, Monitoring and Retrofitting of Concrete Structures, Korea Concrete Institute, 2010.

[47] D.Z. Yankelevsky, A two-phase one dimensional model for steel-concrete interaction, Comput. Struct. 65 (6) (1997) 781–794, http://dx.doi.org/10.1016/S0045-7949(97)00075-8.

[48] D.Z. Yankelevsky, M.A. Adin, D.N. Farhey, Mathematical model for bond-slip behavior under cyclic loading, Struct. J. 89 (6) (1992) 692–698, http://dx.doi.org/10.14359/4143.

[49] H.G. Kwak, S.P. Kim, Bond-slip behavior under monotonic uniaxial loads, Eng. Struct. 23 (3) (2001) 298–309, http://dx.doi.org/10.1016/S0141-0296(00)00008-0.

[50] B. Banholzer, W. Brameshuber, W. Jung, Analytical simulation of pull-out tests—the direct problem, Cem. Concr. Compos. 27 (1) (2005) 93–101, http://dx.doi.org/10.1016/j.cemconcomp.2004.01.006.

[51] S. Khalfallah, M. Ouchenane, A numerical simulation of bond for pull-out tests: the direct problem, Asian J. Civ. Eng. (Buil. Hous.) 8 (5) (2007) 491–505.

[52] M. Haskett, D.J. Oehlers, MS. Mohamed Ali, Local and global bond characteristics of steel reinforcing bars, Eng. Struct. 30 (2) (2008) 376–383, http://dx.doi.org/10.1016/j.engstruct.2007.04.007.

[53] G.M. Verderame, G. De Carlo, P. Ricci, G. Fabbrocino, Cyclic bond behaviour of plain bars, part II: analytical investigation, Constr. Build. Mater. 23 (12) (2009) 3512–3522, http://dx.doi.org/10.1016/j.conbuildmat.2009.07.001.

[54] J. Murcia-Delso, A. Stavridis, B. Shing, Modeling the bond-slip behavior of confined large-diameter reinforcing bars, in: III ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering, 2011.

[55] A. Torre-Casanova, L. Jason, L. Davenne, X. Pinelli, Confinement effects on the steel-concrete bond strength and pull-out failure, Eng. Fract. Mech. 97 (2013) 92–104, http://dx.doi.org/10.1016/j.engfracmech.2012.10.013.

[56] RILEM, Essais portant sur l'adhérence des armatures du béton :2 : Essai par traction, Mater. Struct. 3 (1970) 175–178, http://dx.doi.org/10.1007/BF02478968.

[57] R.I. Gilbert, S. Nejadi, An Experimental Study of Flexural Cracking in Reinforced Concrete Members under Short Term Loads, University of New South Wales, School of Civil and Environmental Engineering, 2004.

[58] J. Mazars, F. Hamon, S. Grange, A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings, Mater. Struct. 48 (2015) 3779–3793, http://dx.doi.org/10.1617/s11527-014-0439-8.

[59] A. Hillerborg, M. Modéer, P.E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cem. Concr. Res. 6 (6) (1976) 773–781, http://dx.doi.org/10.1016/0008-8846(76)90007-7.

[60] Fib model code for concrete structures, 2010, http://dx.doi.org/10.35789/fib.BULL.0055.

[61] J. Murcia-Delso, P. Benson Shing, Bond-slip model for detailed finite-element analysis of reinforced concrete structures, J. Struct. Eng. 141 (4) (2015) 04014125, http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0001070.

[62] J. Mazars, A description of micro and macroscale damage of concrete structure, Eng. Fract. Mech. 25 (1986) 729–737, http://dx.doi.org/10.1016/0013-7944(86)90036-6.

[63] A. Nassima, M. Matallah, Méthode de Régularisation Énergétique Basée Sur l'Approche Crack Band: Limites d'Application et Sources d'Erreurs, Rencontres Universitaires de Génie Civil, 2015.
[36] I. Bitar, N. Benkemoun, P. Kotronis, S. Grange, A multifiber Timoshenko beam with embedded discontinuities, Eng. Fract. Mech. 214 (2019) 339–364, http://dx.doi.org/10.1016/j.engfracmech.2019.03.032.

[37] I. Bitar, P. Kotronis, N. Benkemoun, S. Grange, A generalized Timoshenko beam with embedded rotation discontinuity, Finite Elem. Anal. Des. 150 (2018) 34–50, http://dx.doi.org/10.1016/j.finel.2018.07.002.

[38] C. Lejouad, S. Capdevielle, B. Richard, F. Ragueneau, A numerical model to predict the effect of corrosion on the dynamic behavior of reinforced concrete beams, in: 12th Canadian Conference of Earthquake Engineering, 2019.

[39] N. Khoder, S. Grange, Y. Sieffert, Enhancement of a two-dimensional multifibre beam element in the case of reinforced concrete structures for taking into account the lateral confinement of concrete due to stirrups, Eur. J. Environ. Civ. Eng. 23 (5) (2019) 564–585, http://dx.doi.org/10.1080/19648189.2018.1446364.

[40] N. Viet-Tung, C. Jean-François, S. Karam, A model for thick laminates and sandwich plates, Compos. Sci. Technol. 65 (3-4) (2005) 475–489, http://dx.doi.org/10.1016/j.compsitech.2004.09.014.

[41] P. Zoghipour, K. Hasim, A. Kefal, M. Yildiz, Implementation of shear-locking-free triangular refined zigzag element for structural analysis of multilayered plates with curvilinear fibers, Compos. Struct. 305 (2023) 116462, http://dx.doi.org/10.1016/j.compstruct.2022.116462.

[42] E. Carrera, R. Augello, A. Pagani, D. Scano, Refined multilayered beam, plate and shell elements based on Jacobi polynomials, Compos. Struct. 304 (2023) 116275, http://dx.doi.org/10.1016/j.compstruct.2022.116275.

[43] S. Grange, D. Bertrand, Implicit coupling of heterogeneous and asynchronous time-schemes using a primal approach based on velocity continuity at the subdomain interface, Finite Elem. Anal. Des. 196 (2021) 103604, http://dx.doi.org/10.1016/j.finel.2021.103604.

[44] R. Eligethausen, E.P. Popov, V.V. Bertero, Local bond stress-slip relationships of deformed bars under generalized excitations, 1982, http://dx.doi.org/10.18419/opus-415.

[45] M.H. Harajli, Development/splice strength of reinforcing bars embedded in plain and fiber reinforced concrete, ACI Mater. J. 91 (5) (1994) 511–520, http://dx.doi.org/10.14359/4163.

[46] P. Desnerck, G. De Schutter, L. Taerwe, A local bond stress-slip model for reinforcing bars in self-compacting concrete, fracture mechanics of concrete and concrete strucures-assesment, in: Durability, Monitoring and Retrofitting of Concrete Structures, Korea Concrete Institute, 2010.

[47] D.Z. Yankelevsky, A two-phase one dimensional model for steel-concrete interaction, Comput. Struct. 65 (6) (1997) 781–794, http://dx.doi.org/10.1016/S0045-7949(97)00075-8.

[48] D.Z. Yankelevsky, M.A. Adin, D.N. Farhey, Mathematical model for bond-slip behavior under cyclic loading, Struct. J. 89 (6) (1992) 692–698, http://dx.doi.org/10.14359/4143.

[49] H.G. Kwak, S.P. Kim, Bond-slip behavior under monotonic uniaxial loads, Eng. Struct. 23 (3) (2001) 298–309, http://dx.doi.org/10.1016/S0141-0296(00)00008-0.

[50] B. Banholzer, W. Brameshuber, W. Jung, Analytical simulation of pull-out tests—the direct problem, Cem. Concr. Compos. 27 (1) (2005) 93–101, http://dx.doi.org/10.1016/j.cemconcomp.2004.01.006.

[51] S. Khalfallah, M. Ouchenane, A numerical simulation of bond for pull-out tests: the direct problem, Asian J. Civ. Eng. (Buil. Hous.) 8 (5) (2007) 491–505.

[52] M. Haskett, D.J. Oehlers, MS. Mohamed Ali, Local and global bond characteristics of steel reinforcing bars, Eng. Struct. 30 (2) (2008) 376–383, http://dx.doi.org/10.1016/j.engstruct.2007.04.007.

[53] G.M. Verderame, G. De Carlo, P. Ricci, G. Fabbrocino, Cyclic bond behaviour of plain bars, part II: analytical investigation, Constr. Build. Mater. 23 (12) (2009) 3512–3522, http://dx.doi.org/10.1016/j.conbuildmat.2009.07.001.

[54] J. Murcia-Delso, A. Stavridis, B. Shing, Modeling the bond-slip behavior of confined large-diameter reinforcing bars, in: III ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering, 2011.

[55] A. Torre-Casanova, L. Jason, L. Davenne, X. Pinelli, Confinement effects on the steel-concrete bond strength and pull-out failure, Eng. Fract. Mech. 97 (2013) 92–104, http://dx.doi.org/10.1016/j.engfracmech.2012.10.013.

[56] RILEM, Essais portant sur l'adhérence des armatures du béton :2 : Essai par traction, Mater. Struct. 3 (1970) 175–178, http://dx.doi.org/10.1007/BF02478968.

[57] R.I. Gilbert, S. Nejadi, An Experimental Study of Flexural Cracking in Reinforced Concrete Members under Short Term Loads, University of New South Wales, School of Civil and Environmental Engineering, 2004.

[58] J. Mazars, F. Hamon, S. Grange, A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings, Mater. Struct. 48 (2015) 3779–3793, http://dx.doi.org/10.1617/s11527-014-0439-8.

[59] A. Hillerborg, M. Modéer, P.E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cem. Concr. Res. 6 (6) (1976) 773–781, http://dx.doi.org/10.1016/0008-8846(76)90007-7.

[60] Fib model code for concrete structures, 2010, http://dx.doi.org/10.35789/fib.BULL.0055.

[61] J. Murcia-Delso, P. Benson Shing, Bond-slip model for detailed finite-element analysis of reinforced concrete structures, J. Struct. Eng. 141 (4) (2015) 04014125, http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0001070.

[62] J. Mazars, A description of micro and macroscale damage of concrete structure, Eng. Fract. Mech. 25 (1986) 729–737, http://dx.doi.org/10.1016/0013-7944(86)90036-6.

[63] A. Nassima, M. Matallah, Méthode de Régularisation Énergétique Basée Sur l'Approche Crack Band: Limites d'Application et Sources d'Erreurs, Rencontres Universitaires de Génie Civil, 2015.