# A multiscale steel-concrete interface model for structural applications

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

► To cite this version:

Maryam Trad, Ibrahim Bitar, Stéphane Grange, Benjamin Richard. A multiscale steel-concrete interface model for structural applications. Structures, 2024, 68, pp.107137. 10.1016/j.istruc.2024.107137. hal-04688327

HAL Id: hal-04688327

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

Submitted on 7 Oct 2024

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

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

<figure><img src="image_1.png" /></figure>
# A multiscale steel-concrete interface model for structural applications

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

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

## ARTICLE INFO

### Keywords:

Steel-concrete interface model

Bond model

Reinforced concrete structures

Macro-element

Multiscale model

Static condensation

## ABSTRACT

In this paper, a new multiscale macro-element formulation for the steel-concrete interface modeling is proposed. This element allows for the representation of the behavior of steel and the interface zone surrounding it. It can also model the interfacial bond stresses in between. Compared to conventional interface models of the literature, which employ separate mesh elements for the steel and the interface, utilizing a macro-element to model both components simplifies the creation of a reinforced concrete structure mesh. Moreover, the macro-element equilibrium is solved using a sub-structuring method that aims to reduce the computational cost. At the global level, it is considered as a four-node element linked to two-dimensional and three-dimensional concrete elements. At the local level, an assembly of multiple three-node bar elements with bond stresses is performed. An inner mesh discretization is therefore possible at the local level independently of the global level. The coupling between the two modeling scales is done using a static condensation technique. The formulation of the macro-element is presented in this paper. A selection of numerical examples is provided. The presented applications demonstrate the robustness of the proposed interface model and its capacity to reproduce the experimental behavior of reinforced concrete structural elements.

## 1. Introduction

The robust characterization of the mechanical behavior for reinforced concrete structures when external loads beyond the design level are applied can be a complex goal to achieve since different types of non-linearities can be involved. In fact, to ensure a numerical robustness of the non-linear approaches used for reinforced concrete structures modeling, it is mandatory to consider the different energy dissipation sources: viscous dissipation, numerical dissipation related to the time integration scheme, and especially the material dissipation. The interaction between concrete and steel reinforcement that occurs at the steel-concrete interface contributes to the material dissipation part. Indeed, the functioning of reinforced concrete relies on the stress transfer between steel and concrete through the interface between both materials. As soon as the first cracks appear in concrete, the steel-concrete interface transmits the internal loads to the steel. The consideration of this steel-concrete interface in numerical modeling has therefore a significant influence on the realistic representation of the cracking process in reinforced concrete structures [1]. According to [2], the energy dissipation along the steel-concrete interface may account for up to 15% of the total material energy dissipation.

Several numerical models are proposed in the literature within different frameworks to describe the interface between steel and concrete. A large variety of models define specific finite elements that link steel

and concrete elements (spring elements [3,4], two-dimensional (2D) and three-dimensional (3D) elements [5,6], and joint elements [7-11]). The interface behavior is considered as internal stresses between steel and concrete in other methodological frameworks [12,13]. A kinematic enhancement with the X-FEM method is considered in [14] as an alternative approach. It is important to note that these various models are most often incorporated in fine and detailed 2D and 3D analyses, to improve the prediction of the local and global behavior of a structural element. A description of the steel-concrete interface at the scale of an industrial building is challenging and highly time consuming while considering these approaches.

More recent works are proposed in order to consider the interface behavior for large-scaled structures, such as homogenization techniques of the steel, the interface, and the concrete behaviors [15]. In parallel, kinematic enhancement approaches within fiber reinforced concrete elements [16,17] are proposed. Nevertheless, using homogenization techniques of reinforced concrete modeling may be challenging to provide a detailed localized description of the concrete cracking around the steel reinforcement [18,19]. In parallel, even though fiber elements make it possible to study large-scaled structures with a relatively short calculation time, full 3D simulations provide more refined and exact results with respect to fiber elements [20].

<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 at: Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-EXP/SES/LMAPS, B.P. 17 - 92262 Fontenay-aux-Roses Cedex, France. E-mail address: maryamtrad100@hotmail.com (M. Trad).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Concept of studying a reinforced concrete structure with the macro-element approach.</figcaption></figure>

The main objective of this work is to address the challenges of the existing interface models of the literature by improving the way the steel-concrete interface is modeled in the context of finite element computations, for large-scale structures. Moreover, multiscale modeling strategies have proven their numerical efficiency for a large variety of applications such as concrete damage [21] and steel-concrete composite components [22–24]. For these reasons, a multiscale macro-element formulation with internal degrees of freedom has been adopted. The element formulation is initially considered in [25] to represent a rigid-inclusion embedded in a soil volume. The formulation of [25] is developed further here and enriched so it can be used for the steel-concrete interface modeling. This element represents the behavior of steel, an interface zone surrounding the steel bar, and concrete-steel bond stresses. Hence, the reinforced concrete structure geometry is discretized into concrete elements and macro-elements. The steel and the interface behavior are considered at the level of the macro-elements.

This work is proposed in the framework of 2D and 3D reinforced concrete models, where distinct constitutive laws are used for steel and concrete. It presents an alternative interface modeling approach of the coupling elements method of [26] where the construction of the finite elements mesh may be challenging since it needs the identification of the finite elements crossed by steel bars to be defined as coupling elements.

This paper is organized as follows. Section 2 recalls the principles of the finite elements resolution using the macro-element. Section 3 presents model applications for pull-out tests, a reinforced concrete tie rod test, a beam-end test, and two beam bending test examples. The numerical behavior of the macro-element formulation is compared to some existing interface models. Indeed, a higher robustness is considered for the macro-element model. The paper ends with conclusions and perspectives in Section 4.

## 2. Theoretical background

In a 2D or 3D framework, with a perfect steel-concrete bond, one can mesh the reinforced concrete structure with 2D/3D concrete elements and 1D steel elements. The interest of this work is to model the interface behavior in a non-intrusive way with respect to perfect bond simulations. For that, 1D macro-elements that represent steel and interface behaviors with bond stresses (see Fig. 1) are defined instead of 1D steel elements. Unlike coupling elements methods that require identifying and replacing concrete elements intersected by steel bars with coupling elements [26], this macro-element technique eliminates the need for such identification steps, providing a practical and non-intrusive way to simulate bond behavior.

Two scales are considered when using the macro-element formulation: a global scale of the whole structure, and a local scale at the level of the macro-element. Two complementary Newton-Raphson resolution algorithms are adopted at the two distinct levels. The boundary value problems are here detailed for the global and the local levels. The link between the two levels is ensured by means of a static condensation technique that is also presented in this section.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Ω domain.</figcaption></figure>

### 2.1. Global boundary value problem

The global boundary value problem consists in considering a solid body that occupies a volume $Ω ⊂ R^n$ (with $n = 1, 2, 3$ the problem's dimension and $⊂$ the proper subset symbol). $∂Ω ⊂ R^n$ is the boundary of $Ω$. Density forces $f$ are applied to $Ω$. $∂Ω$ is composed of two complementary parts $∂Ω_f ⊂ ∂Ω$ and $∂Ω_u ⊂ ∂Ω$. Surface forces $t$ are applied to $∂Ω_f$, and displacement values $U_d$ are imposed along $∂Ω_u$ (see Fig. 2).

#### 2.1.1. Strong formulation

Under the small perturbations assumption, solving the boundary value problem consists in finding the admissible displacement field $U$ and the statically admissible stress tensor $\sigma$ at a position $x$ such that:

$$ 
\left\{ 
\begin{array}{l}
\nabla \cdot \sigma + \mathbf{f} = 0 \quad \forall \mathbf{x} \in \Omega \quad : \text{Equilibrium equations} \\
\sigma = \mathcal{F}(\epsilon) \quad \forall \mathbf{x} \in \Omega \quad : \text{Stress-strain (constitutive) equations} \\
\epsilon = \nabla^s U \quad \forall \mathbf{x} \in \Omega \quad : \text{Strain-displacement (kinematic) equations} \\
\sigma \mathbf{n}_f = \mathbf{t} \quad \forall \mathbf{x} \in \partial \Omega_f \quad : \text{Neumann boundary conditions} \\
U = U_d \quad \forall \mathbf{x} \in \partial \Omega_u \quad : \text{Dirichlet boundary conditions}
\end{array}
\right. 
\quad (1) 
$$

Where:

* $\nabla^s$: is the symmetric part of the gradient operator $\nabla$.

* $\mathbf{n}_f$: is a normal vector to $∂Ω_f$.

* $\mathcal{F}$: is a constitutive relationship.

* $\forall$: is the universal quantification symbol.

* $\in$: is the set membership symbol.

#### 2.1.2. Weak formulation

The virtual power principle leads to the variational formulation of the problem which is formulated as follows:

$$ \int_{\Omega} \varepsilon^* : \sigma dV = \int_{\Omega} U^{*T} f dV + \int_{\partial \Omega_f} U^{*T} t dS + \int_{\partial \Omega_u} U_d^{*T} (\sigma \mathbf{n}_u) dS \quad (2) $$

Where:

* $U^*$: is a virtual degrees of freedom field, equal to $U_d^*$ along $∂Ω_u$.

* $\mathbf{n}_u$: is a normal vector to $∂Ω_u$.

* $T$: is the transpose operator.

By choosing $U^*$ as a kinematically admissible displacement field equal to zero along $∂Ω_u$, Eq. (2) can be rewritten as follows:

$$ \int_{\Omega} \varepsilon^* : \sigma dV = \int_{\Omega} U^{*T} f dV + \int_{\partial \Omega_f} U^{*T} t dS \quad (3) $$

Such as:

* $\int_{\Omega} \varepsilon^* : \sigma dV$: is the virtual power of all the internal forces.

* $\int_{\Omega} f U^{*T} dV + \int_{\partial \Omega_f} t U^{*T} dS$: is the virtual power of the external forces.

Eq. (3) states the principle of virtual power where the power $P_{ext}^*$ of the external forces applied to a system is equal to the power $P_{int}^*$ of its internal forces:

$$ P_{int}^* = P_{ext}^* \quad (4) $$
### 2.1.3. Discretization

A finite elements discretization of $\Omega$ is considered. $\Omega$ is discretized into $N_e$ elements $\Omega_e$. $\partial\Omega$ is then divided into $N_{\Omega_e}$ parts called $\partial\Omega_e$. It is important to note here that when studying a reinforced concrete structure, $\Omega_e$ can be either a concrete element $\Omega_{ec}$ or a macro-element $\Omega_{em}$ that represents steel and steel-concrete interface, such that $\Omega_e = \Omega_{ec} \cup \Omega_{em}$ and $\Omega_{ec} \cap \Omega_{em} = \emptyset$. $\cup$ and $\cap$ are the union and the intersection operators. $\emptyset$ stands for the empty set. For concrete elements, the displacement $U_i$ of a point $i$ located inside an element can be deduced from the nodal displacement field $U$ according to the following equation:

$$U_i = N U \quad (5)$$

Where:

* $N$ is the matrix holding the set of interpolation functions attributed to the finite elements used to construct the mesh.

* $U$ is the vector of all the degrees of freedom of the mesh nodes.

The strain field $\epsilon$ is derived from the generalized displacements field $U$ as follows:

$$\epsilon = \nabla^s U_i = B U \quad (6)$$

Where $B$ denotes the strain-displacement matrix, holding the set of the first derivatives of the shape functions $N$, and $\nabla^s$ stands for the symmetric part of the gradient operator.

The virtual internal forces power is a combination of two contributions $P_{int}^*$ and $P_{intm}^*$ for concrete-type and macro-element-type elements. Hence, the expression of the total power of internal forces $P_{int}^*$ is expressed as follows:

$$P_{int}^* = P_{intc}^* + P_{intm}^* \quad (7)$$

The virtual power principle of Eq. (4) can be expressed as follows:

$$U^{*T} F_{int}(U) = U^{*T} F_{ext} \quad \forall U^{*} \quad (8)$$

The internal forces vector $F_{int}$ is a combination of the internal forces of concrete-type and macro-element-type elements:

$$F_{int} = A_{c=1}^{N_e} \left[ \int_{\Omega_{ec}} B^T \sigma(\epsilon) dV \cup f_b \right] \quad (9)$$

Where the symbol $A$ denotes the assembly operator, and $\sigma$ is the concrete stress tensor that depends of the strain $\epsilon$. The expression of $\sigma(\epsilon)$ is given by the concrete constitutive law. The vector $F_{ext}$ of the external forces is expressed as:

$$F_{ext} = A_{e=1}^{N_e} \int_{\Omega_e} fdV + A_{e=1}^{N_{\Omega_e}} \int_{\partial\Omega_e} tdS \quad (10)$$

To solve Eq. (8) (after discretization in time) a classical Newton Raphson algorithm can be used. The linearization of the internal forces vector $F_{int}(U)$ of Eq. (8) gives:

$$dF_{int} = kdU \quad (11)$$

Where:

$$k = \frac{\partial F_{int}}{\partial U} = A_{e=1}^{N_e} \left[ \int_{\Omega_{ec}} B^T C B dV \cup k_{em} \right] \quad (12)$$

$dU$ is the vector of the incremental values of the total degrees of freedom of the studied structure and $dF_{int}$ is the incremental internal forces vector. $C$ is the concrete constitutive law matrix. $k$ is the stiffness matrix of the studied structure. $\int_{\Omega_{ec}} B^T C B dV$ is an elementary stiffness matrix of a concrete element, and $k_{em}$ is an elementary stiffness matrix of a macro-element. The calculation of $k_{em}$ requires a local equilibrium at the level of the macro-element, which is detailed in the next subsection.

Imposed displacements along $\partial\Omega_u$ can be expressed as follows:

$$LU = U_d \quad (13)$$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Interaction between the steel and concrete interface domains through internal stresses $\tau_{int}$.</figcaption></figure>

Where $L$ is the matrix of blocking, and $U_d$ is the vector of the imposed values of all the degrees of freedom. The finite elements problem consists in resolving both Eqs. (11) and (13). For that, the well-known double Lagrange multipliers method can be used to enforce the Dirichlet boundary conditions [27], as shown in Appendix A.

## 2.2. Local boundary value problem

Each macro-element is an assembly of multiple biphasic elements. The formulation of one biphasic element is studied here in the first place.

### 2.2.1. Study of a biphasic element

Let $H$ be the length of the element. $(E_s, S_s)$ and $(E_i, S_i)$ are the Young's modulus and the cross-sections for the steel bar and the interface domains respectively. Two independent virtual longitudinal strain fields, $\epsilon_s^*$ and $\epsilon_i^*$, are associated to the steel and to the interface. The friction stress at the interface is considered as an external contribution as follows:

$$\tau_i (y_s - y_i) = -\tau_s (y_s - y_i) \quad (14)$$

Where $\tau_i$ is the bond stress applied on the interface, while $\tau_s$ is applied on the steel. The friction is a linear or a nonlinear function of the differential displacement $(y_s - y_i)$. $y_s$ and $y_i$ are the longitudinal displacements along the steel and the interface zones, respectively. Additional external forces may be considered, and the associated virtual power is here denoted $P_{ext}^*$ (see Fig. 3).

The principal of virtual power is then expressed as:

$$\int_{0}^{H} (\epsilon_{s}^{*} \sigma_{s}(\epsilon_{s}) S_{s} + \epsilon_{i}^{*} E_{i} \epsilon_{i} S_{i}) dx + \int_{0}^{H} (y_{s}^{*} - y_{i}^{*}) \tau_{i} (y_{s} - y_{i}) P dx = P_{ext}^{*} \quad (15)$$

Where $P$ is the steel bar perimeter. It is important to state here that in the current version of the macro-element model, and in the sake of simplicity, a linear constitutive law is associated to the interface zone. The interface nonlinearities are concentrated at the level of the friction stresses between the two distinct zones that can have a nonlinear evolution with respect to the steel-interface slip. For this reason, the stress at the interface, denoted $\sigma_i$, is expressed in Eq. (31) as $E_i \epsilon_i$ (a linear function). However, considering a nonlinear constitutive law for the interface zone remains possible. On the other hand, steel can have a linear or a nonlinear behavior depending on the studied structure and the applied loads. For that, the steel stress $\sigma_s$ is a function of the steel strain $\epsilon_s$ defined by the steel constitutive law.

The first term of Eq. (31) defines the virtual internal power as:

$$P_{int}^* = \int_{0}^{H} (\epsilon_s^* \sigma_s(\epsilon_s) S_s + \epsilon_i^* E_i \epsilon_i S_i) dx + \int_{0}^{H} (y_s^* - y_i^*) \tau_i (y_s - y_i) P dx \quad (16)$$

### 2.2.2. Discretization

Eq. (32) can be discretized. Three-node bar elements are chosen for the discretization. It is important to note that this choice imposes
that the external forces should be along the macro-element direction. However, it is possible to set a more complex formulation in which beam elements are used instead of bar elements, so density forces can be applied between the steel and the interface parts of the biphasic element perpendicularly to its direction. The current version of the macro-element focuses on the longitudinal steel-concrete behavior. No normal degrees of freedom are defined for the inner nodes of macro-elements. The longitudinal displacements along the steel and the interface zone are deduced respectively as:

$$ \left\{ \begin{array}{l} y_s(x) = N(x)u_s \\ y_i(x) = N(x)u_i \end{array} \right. \qquad (17) $$

$u_s$ and $u_i$ are the steel and the interface nodal longitudinal displacements. Here are the nodes positions for a three-node bar element: $x_{nodes} = \{0; H/2; H\}$; for $x \in [0 : H]$. The matrix of interpolation functions $N$ is expressed as:

$$ N = \left[ \frac{1}{2} \left( \frac{4x^2}{H^2} - \frac{6x}{H} + 2 \right) - \frac{4x^2}{H^2} + \frac{4x}{H} - \frac{2x^2}{H^2} - \frac{x}{H} \right] \qquad (18) $$

A geometric transformation is used in order to simplify the $N$ matrix, and to express it in the reference space as a function of the variable $\xi \in [-1 : 1]$. The new variable $\xi$ is related to $x$ according to the following equation:

$$ x = \frac{1 + \xi}{2} H \qquad (19) $$

The matrix of shape functions $N$ expressed in terms of $\xi$ as follows:

$$ N = \left[ \frac{1}{2}\xi(\xi - 1) \quad 1 - \xi^2 \quad \frac{1}{2}\xi(\xi + 1) \right] \qquad (20) $$

The $B$ matrix of the derivatives of the interpolation functions is used to calculate $\frac{\partial y_s}{\partial x}$ and $\frac{\partial y_i}{\partial x}$ with respect to $u_s$ and $u_i$. Thus:

$$ B = \frac{\partial N}{\partial \xi} \frac{\partial \xi}{\partial x} = \left[ \xi - \frac{1}{2} \quad -2\xi \quad \xi + \frac{1}{2} \right] \frac{2}{H} \qquad (21) $$

Therefore, Eq. (32) is expressed as follows:

$$ P_{int}^* = \int_0^H \mathbf{u}_s^{*T} \mathbf{B}^T S_s \sigma_s (\mathbf{B} \mathbf{u}_s) + \mathbf{u}_i^{*T} \mathbf{B}^T E_i S_i \mathbf{B} \mathbf{u}_i dx \\ + \int_0^H (\mathbf{u}_i^* - \mathbf{u}_s^*)^T \mathbf{N}^T \tau_i (y_i - y_s) P dx \qquad (22) $$

Since $\mathbf{u}_i^*$ and $\mathbf{u}_s^*$ are defined as two different independent virtual displacement vectors, the internal forces field $p^{el}$ for one biphasic element composed of two three-node bar elements representing the steel and the interface domain with the internal friction stresses is formulated as:

$$ p^{el} = \left[ \begin{array}{l} \int_0^H \mathbf{B}^T S_s \sigma_s (\mathbf{B} \mathbf{u}_s) dx \\ \int_0^H \mathbf{B}^T E_i S_i \mathbf{B} \mathbf{u}_i dx \\ + \int_0^H \left[ \frac{-\mathbf{N}^T}{\mathbf{N}^T} \right] \frac{\partial \tau_i}{\partial (y_i - y_s)} [-\mathbf{N} \quad \mathbf{N}] P dx \end{array} \right] \qquad (23) $$

On the other hand, the subtraction of the steel and the interface displacements $y_s(x)$ and $y_s(x)$ (see Eq. (17)) gives:

$$ y_i - y_s = [-N \quad N] \left[ \begin{array}{c} u_s \\ u_i \end{array} \right] = [-N \quad N] u^{el} \qquad (24) $$

After introducing Eq. (24) into Eq. (23), the derivative of Eq. (23) with respect to the elementary degrees of freedom vector $u^{el}$ returns:

$$ \frac{\partial p^{el}}{\partial u^{el}} = k_{bp}^{el} = \int_{0}^{H} \left[ \begin{array}{cc} \mathbf{B}^{T} C_s S_s \mathbf{B} & 0 \\ 0 & \mathbf{B}^{T} E_i S_i \mathbf{B} \end{array} \right] dx \\ + \int_{0}^{H} \left[ \frac{-\mathbf{N}^{T}}{\mathbf{N}^{T}} \right] \frac{\partial \tau_i}{\partial (y_i - y_s)} [-\mathbf{N} \quad \mathbf{N}] P dx \qquad (25) $$

$k_{bp}^{el}$ is the elementary stiffness matrix of one biphasic element. $C_s$ is the steel constitutive law matrix. $\frac{\partial \tau_i}{\partial (y_i - y_s)}$ is calculated due to the expression of a steel-interface constitutive bond law which links the stress value $\tau_i$ to the steel-concrete slip $y_s - y_i$. Pull-out tests are usually adopted to characterize the relation between steel-concrete bond stress and the corresponding slip. Several analytical expressions for constitutive bond laws are proposed in the literature [3,4,13] [28–37]. The choice of the bond law will be detailed and justified for each application example in the next section.

## 2.3. Coupling between the global and the local levels

The macro-element can be seen at different scales: a global and a local one. At the global level, it is considered as a four-node element. An inner discretization can be done at the local level. More precisely, each macro-element is an assembly of a set of biphasic elements.

### 2.3.1. Connection between the macro-elements and the concrete elements meshes

At the global level, 2D or 3D elements can be associated to concrete. Two distinct mesh configurations are possible: coincident or non-coincident meshes for concrete and macro-elements (see Fig. 4). In case of coincident meshes, at global level, both nodes of the interface part of the macro and concrete elements are superimposed and thus the displacements are considered equal. At local level, in the normal direction of the macro-element, linear kinematic relations are imposed so that the displacements of the steel nodes are equal to the displacements of the interface nodes. In the longitudinal direction, however, the steel nodes are free to move independently of the interface node displacements. So the steel-interface slip can occur in the longitudinal direction.

In case of a 3D simulation with non-coincident meshes, at global level, the displacements of the interface nodes are linked to displacements of the surrounding concrete nodes, denoted by the subscript $c_i$, with the shape functions $N_i$. When using eight-node elements for concrete mesh these kinematic relations are defined such that:

$$ \left\{ \begin{array}{l} u_{\text{interface},x} = \sum_{i=1}^{8} N_i u_{c_i,x} \\ u_{\text{interface},y} = \sum_{i=1}^{8} N_i u_{c_i,y} \\ u_{\text{interface},z} = \sum_{i=1}^{8} N_i u_{c_i,z} \end{array} \right. \qquad (26) $$

For 2D simulations where four-node elements constitute the concrete mesh, four shape functions $N_i$ are defined instead of defining eight shape functions. Steel nodes displacements are similarly connected to surrounding concrete nodes in the normal directions with respect to the macro-element. No connections are imposed along the longitudinal direction, so longitudinal interface-steel slip can occur.

One possible approach to introduce these displacements relations is to enrich the matrix $k_{tot}$ of Eq. (A.5) with new kinematic relations (Eq. (26)) by applying the double Lagrange multipliers methodology. This Lagrange method, even though it is widely used to add kinematic relations to a finite elements resolution system, adds additional unknowns to be calculated (the Lagrange multipliers) [38]. Using this approach is convenient when considering Dirichlet boundary conditions (which are technically kinematic relations). In contrast, linking steel and interface zones to concrete adds kinematic relations for all the steel and the interface nodes. An alternative method is proposed here to incorporate steel/interface-concrete kinematic relations in the resolution. This approach based on a kinematic projection principle is here described.

The kinematic projection methodology consists of classifying the total degrees of freedom $U_{tot}$ (see Eq. (A.3) of Appendix A) into two complementary vectors such as:

$$ U_{tot} = [ U_c \quad U_i ]^T \qquad (27) $$

Where:

* $U_c$ is the vector of the dependent degrees of freedom.

* $U_i$ is the vector of the independent degrees of freedom.

The kinematic relations impose that the values of the displacements of the steel and the interface nodes are dependent of the values of the displacements of the surrounding concrete nodes, which allows to
<figure><img src="image_6.png" /><figcaption>Fig. 4. Global and local levels.</figcaption></figure>

define $U_c$ as the vector of the displacements of the interface nodes in the three dimensions of the space and the displacements of the steel nodes in the normal directions with respect to the steel bars directions.

By deriving the kinematic relations equations, the vectors of the incremental values of the total degrees of freedom and the independent ones can be linked with a kinematic projection matrix $P$:

$$ \delta U_{tot} = P \delta U_i \quad (28) $$

The projection matrix here $P$ holds the sets of the derivatives of the kinematic relations. Multiplying Eq. (A.2) of Appendix A by the transpose of the vector $\delta U_{tot}$ gives:

$$ \delta U_i^T P^T F(U_{tot}) = \delta U_i^T P^T F_{tot} \quad (29) $$

So,

$$ P^T F(U_{tot}) = P^T F_{tot} \quad (30) $$

Eq. (30) is the equilibrium equation to be resolved. The residue $R$ associated with the equilibrium equation is defined as:

$$ R = P^T F_{tot} - P^T F(U_{tot}) \quad (31) $$

The iterative resolution of the global Newton Raphson algorithm aims to minimize the value of the residue $R$. It incorporates the derivation of the residue calculated as:

$$ \frac{\partial R}{\partial U_{du_{b,i}}} = \frac{\partial (P^T F_{tot})}{\partial U_i} - \frac{\partial (P^T F(U_{tot}))}{\partial U_i} \quad (32) $$

The residue $R$ is derived with respect to $U_i$. The advantage of this kinematic projection approach is that only $U_i$ is calculated at each resolution iteration. $U_c$ is deduced using the kinematic relations. Eq. (32) is developed as:

$$ \frac{\partial R}{\partial U_i} = - \left[ P^T \frac{\partial F(U_{tot})}{\partial U_{tot}} \frac{\partial U_{tot}}{\partial U_i} \right] \quad (33) $$

Knowing that the linear expressions of the used kinematic relations impose that $\frac{\partial P}{\partial U_i}$ is equal to zero. $\frac{\partial U_{tot}}{\partial U_i}$ is equal to the projection matrix $P$ (see Eq. (A.1) of Appendix A), so:

$$ \frac{\partial R}{\partial U_i} = - \left[ P^T \frac{\partial F(U_{tot})}{\partial U_{tot}} P \right] \quad (34) $$

$\frac{\partial R}{\partial U_i}$ is the tangent operator of the nonlinear global Newton Raphson resolution (where the Dirichlet boundary conditions and the kinematics relations that link the steel/interface displacements to the concrete nodes displacements are taken into account).

Fig. 5 illustrates the way the study of a reinforced concrete structure is designed.

The studied structure is discretized into two types of finite elements: concrete elements and macro-elements. Each macro-element is locally a set of assembled biphasic elements. One biphasic element is a combination of two three-node bar elements and shear interfacial stresses in between. The two three nodes bars represent steel and interface zones.

It is important to highlight the advantage of the macro-element methodology regarding the mesh construction. Indeed, in a 3D framework, with a perfect bond steel-concrete behavior, the reinforced concrete structure geometry is discretized with 2D/3D concrete elements and steel with 1D elements. The steel nodes are linked to concrete nodes with kinematic relations. The interest of the macro-element formulation is about modeling the interface behavior within this framework in a non-intrusive way. For that, one defines macro-elements instead of 1D steel elements. The main difference with respect to the perfect bond case is about defining two nodes at each longitudinal position of the steel bars (interface and steel nodes). These nodes have the same coordinates at the initial configuration of the structure.

### 2.3.2. Inner resolution of the macro-element equilibrium

In the present work, the coupling between the two levels is done by adopting a static condensation technique. The aim of this technique is to decrease the total number of degrees of freedom to be computed at the global scale. Indeed, in the global resolution, four degrees of freedom are considered for each macro-element. These degrees represent the longitudinal displacements of the four outer nodes of the inner local discretization. $u_{tot}^T = [u_r \ u_b]^T$ is the total degrees of freedom vector of the biphasic elements constituting one macro-element. $r$ and $b$ subscripts refer to internal and external degrees of freedom. $u_b$ represents the degrees of freedom of the four outer nodes, and $u_r$ is the internal degrees of freedom vector of the inner nodes.

The resistant forces vector $f_{tot}$ of a macro-element is an assembly of the elementary internal forces vectors $p^{el}$ (see Eq. (23)). Let $f_r$ be the forces vector at the level of the inner nodes of a macro-element, and $f_b$ be the forces at the level of the outer nodes (see Eq. (9)). Hence, it is possible to define $f_{tot}$ such that $f_{tot}^T = [f_r \ f_b]^T$. The inner resolution at the local level of the macro-element consists in resolving the following equation:

$$ f_r = 0 \quad (35) $$

Eq. (35) states the inner equilibrium at the level of the inner discretization of the macro-element.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Reinforced concrete structure study.</figcaption></figure>

Let $k_{bp}$ be the assembly of the elementary stiffness matrices $k_{bp}^{el}$ of biphasic elements constituting one macro-element (see Eq. (25)). It is possible to state that:

$$df_{tot} = k_{bp}du_{tot} \qquad (36)$$

Where $du_{tot}$ represents the vector of the incremental values of the total degrees of freedom of the inner discretization of the macro-element. $du_{tot}^T = [ du_r du_b ]^T$ where $du_r$ and $du_b$ are the vectors of the incremental values of $u_r$ and $u_b$. According to these degrees of freedom subscripts r and b (internal and external degrees), the matrix $k_{bp}$ is composed of four parts: $k_{rr}$, $k_{rb}$, $k_{br}$, and $k_{bb}$. In order to calculate the algorithmic tangent operator, it is necessary to differentiate $f_b$ according to $u_b$. Hence, the differentiation of the internal efforts vector $f_b$ gives:

$$ \begin{bmatrix} k_{bb} & k_{br} \\ k_{rb} & k_{rr} \end{bmatrix} \begin{bmatrix} du_b \\ du_r \end{bmatrix} = \begin{bmatrix} df_b \\ 0 \end{bmatrix} \qquad (37) $$

A static condensation provides a link between the incremental vectors $du_b$ and $df_b$ as:

$$(k_{bb} - k_{br}k_{rr}^{-1}k_{rb}) du_b = df_b \qquad (38)$$

In a more condensed form, Eq. (38) is written as:

$$k_{em}du_b = df_b \qquad (39)$$

Where $k_{em}$ is the condensed elementary stiffness matrix of one macro-element (see Eq. (12)).

The independence between the mesh sizes at the local and the global level of the proposed approach makes it possible to refine the study at the interface local level. Refining the mesh at the level of the macro-element has a slight numerical cost in comparison to refining the 3D concrete mesh. The local mesh density is a user choice. In addition, the local discretization accelerates the global resolution of the studied structure equilibrium. This resolution acceleration is due to the sub-structuring resolution adopted methodology using a global and a local resolution algorithms.

## 2.4. Calibration of material parameters

In order to define the biphasic element, the behavior laws for the steel and the interface zones have to be defined in addition to the bond law that defines the bond stress. The steel behavior law can

<table><caption>Table 1
Main calibration parameters of a biphasic element.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>i</sub></td><td>Steel Young's modulus</td><td>GPa</td></tr><tr><td>S<sub>i</sub></td><td>Steel bar section</td><td>m<sup>2</sup></td></tr><tr><td>E<sub>f</sub></td><td>Interface Young's modulus</td><td>GPa</td></tr><tr><td>S<sub>f</sub></td><td>Interface zone section</td><td>m<sup>2</sup></td></tr></tbody></table>

follow a linear law or an elastic-plastic law depending on the studied structure. For structures where steel bars do not plastify, considering a linear behavior law for steel is satisfactory. The bond laws relations and parameters are usually identified based on pull-out experimental tests. For the interface zone, and for the sake of simplicity, it is supposed here to have the same material parameters as concrete. For the same reason, a linear behavior law is considered for this zone. However, changing the interface zone parameters and constitutive law remains possible. If a nonlinear constitutive law is being considered for the interface zone, the term $E_i S_i B u_i$ of Eq. (22) has to be replaced by $S_i \sigma(B u_i)$ where $\sigma(B u_i)$ is the stress in the interface linked to the strain $B u_i$ via the adopted nonlinear law. The thickness to be associated to the interface zone is a parameter of each biphasic element. This thickness parameter has a physical interpretation since the experimental behavior of the interface shows that a specific concrete zone around the steel is mostly damaged, and this zone represents the interface. Sensitivity studies will be carried out in the next section to study the effect of the interface thickness on the structural behavior of reinforced concrete elements.

Table 1 recalls the main biphasic element parameters. Additional parameters can be added to describe non-linear constitutive laws for steel an interface zones. These parameters depend in this case on the chosen constitutive laws. In addition, a nonlinear constitutive bond law is usually used. The parameters of this bond law depend (in terms of number of parameters and physical interpretation of each parameter) on the chosen bond law.

## 2.5. Global resolution algorithm

The finite element resolution includes two Newton Raphson algorithms: one at the global level of the whole studied structure, and
another one at the local level of the macro-element. The global Algorithm 1 is called for each time step n.

**Algorithm 1 Global resolution algorithm**

$U_{tot,n+1} \leftarrow U_{tot,n}$ (=0 for initial time step), see equation (A.3) for the definition of $U_{tot}$

$conv \leftarrow 0$

$k \leftarrow 0$

$U_{tot,k} = U_{tot,n+1}$
while $conv = 0$ do

Compute $R_k = F_{tot,n} - F_k(U_{tot,k})$ (see equation (A.2)). At this step, the local algorithm 2 is called to calculate the contribution of each macro-element in $F_k$ (represented by the vector $f_b$ in algorithm 2).

Compute $U_{tot,k+1} = U_{tot,k} + (k_{tot,k})^{-1} R_k$ (see equation (A.5)). At this step, the local algorithm 2 is called to calculate the contribution of each macro-element in $k_{tot}$ (represented by the matrix $k_{em}$ in algorithm 2).

Compute $R_{k+1} = F_{tot,n} - F_{k+1}(U_{tot,k+1})$

if $R_{k+1} < tolerance$ then

$conv \leftarrow 1$

$U_{tot,n} = U_{tot,k+1}$

else

$k \leftarrow k + 1$

End if

End while

It is important to note that the global Algorithm 1 is a classical Newton Raphson algorithm. When calculating the stiffness and the internal forces of the macro-elements, the local Algorithm 2 is used (which also has a classical Newton Raphson structure). In other terms, the macro-element can be easily implemented in any finite elements resolution code with no change in its structure.

## 2.6. Local resolution algorithm

The local Algorithm 2 is at the level of each macro-element. $k_i$ is the index used for the local iterations.

**Algorithm 2 Local resolution algorithm**

$conv_i \leftarrow 0$

$k_i \leftarrow 0$

$u_{tot}^T = [u_r \ u_b]^T$. $u_b$ is given by the global algorithm. It remains constant in this local one. $u_r$ is firstly estimated equal to its value for the previous converged time step

while $conv_i = 0$ do

Compute the elementary quantities $p^{el}$ and $\frac{\partial p^{el}}{\partial u^{el}}$ for each biphasic element of the macro element using equations (23) and (25)

Compute $f_{tot} = A_{el=1}^{N_{bp}} p^{el}$ and $k_{bp} = A_{el=1}^{N_{bp}} \frac{\partial p^{el}}{\partial u^{el}}$

Compute $f_r = f_{tot,r}$

Compute $R_{k_i} = -f_r$

if $R_{k_i} < tolerance$ then

$conv_i \leftarrow 1$

Static condensation: $k_{em} = k_{bb} - k_{br}k_{rr}^{-1}k_{rb}$ (see equation (38)) and $f_b = f_{tot,b}$. $k_{em}$ and $f_b$ are the rigidity matrix and the internal forces vector of the macro-element sent back to be used in the global algorithm

else

$$u_{r_{k_i+1}} = u_{r_{k_i}} + k_r^{-1} R_{k_i} \\ k_i \leftarrow k_i + 1$$

End if

End while

<figure><img src="image_8.png" /><figcaption>Fig. 6. Normalized pull-out test.</figcaption></figure>

## 3. Applications and results

The macro-element interface model is here adopted for three distinct application examples: a pull-out test, a tie-rod example, and a beam test.

### 3.1. 1D pull-out test model

In this subsection, the pull-out experimental test of [13] is modeled using macro-elements, as a validation test of the macro-element formulation. For this 1D model, the macro-element model itself represents the steel, the concrete, and the steel-concrete bond stresses. In other terms, the interface part of the macro-element represents the whole concrete volume. The loading configurations studied are monotonic and cyclic.

The test geometry is composed of a concrete cube crossed by a single steel reinforcement bar. The translation of the concrete cube is blocked by a metal plate. The specimen also rests on a Teflon support to ensure a proper alignment between the steel bar and the direction of the imposed displacement, thus preventing the reinforcement from bending. The concrete-steel slip is measured by a first LVDT (Linear Variable Differential Transformer) sensor located at one edge of the reinforcement and a second LVDT sensor is located at the other edge of the steel bar in order to access to a second relative displacement value between steel and concrete. The contact length is equal to five times the steel bar diameter which is denoted $d_a$. This contact length value is recommended in [39]. The adhesion value $\tau$ is supposed to be constant along the steel bar and is estimated as follows:

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

Where F is the measured reaction and L is the steel-concrete contact length. 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. 6).

#### 3.1.1. Material properties

For this 1D model, linear elastic constitutive laws are considered for steel and concrete (the interface zone of the macro-element here). There is only one source of non-linearity in the problem, which arises from the non-linear expression of the bond stress law at the steel-concrete interface.

The bond law of [40] is here used (see Appendix B). This law was initially proposed to model the bond-slip behavior of confined large-diameter steel reinforcing bars. However, it can be used for the general case of any bar diameter, and applied for confined or unconfined bars, by calibrating its parameters in an appropriate manner. This law is chosen for its simplicity. It can be applied for monotonic and cyclic loading configurations, and is defined by three parameters: 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. Table 2 sums up the material properties used for the 1D pull-out model.

It has to be noted that the bond law parameters are chosen to best describe the experimental bond law curve of [13] (see Fig. 7).
<table><caption>**Table 2**
Material properties for the 1D 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. B.1(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. B.1(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. B.1(a))</td><td>10</td><td>mm</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Analytical and experimental adhesion laws.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** 1D pull-out model.</figcaption></figure>

### 3.1.2. 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. One macro-element represents the test. This macro-element is discretized into three biphasic elements. Fig. 8 shows the mesh and the boundary conditions. For this 1D example, the cross-section of the interface part of the macro-element is taken equal to the section of the concrete cube specimen of the pull-out test, which is equal to $15d_a \times 15d_a$, where $d_a$ refers to the steel bar diameter equal to 12 mm.

### 3.1.3. Results

Figs. 9 and 10 illustrate the imposed displacement and the reaction curves for monotonic and cyclic loads. Even though no experimental result is given in [13] for cyclic pull-out tests, the law of [40] can reproduce the cyclic behavior with no additional input parameters. For the cyclic case boundary conditions are modified between the two configurations of Fig. 11 every time the imposed displacement is equal to zero.

### 3.1.4. Discussion

The reaction force curves show that the constitutive expression that links the reaction force F to the bond stress $\tau$ (Eq. (40)) is fulfilled.

The advantage of the inner discretization is that it allows to have an access to internal variables such as the displacements of the internal nodes and the internal stresses in steel and concrete. For instance, Figs. 12 and 13 show the forces in the concrete and steel elements and

the displacements of the inner and the outer concrete and steel nodes and for all the time steps.

The 1D pull-out model is a simplified yet representative model of the pull-out test. Nevertheless, this model does not represent a detailed 3D description of the test and does not show the distribution of the damage in the concrete volume. Hence, a 3D pull-out model is studied and presented in the next paragraph.

### 3.2. 3D pull-out test model

The pullout of test [13] modeled in Section 3.1 within a 1D configuration is here modeled in a 3D one. For this 3D model, steel bar elements, an interface zone with a certain thickness, and bond stresses, and are attached to 3D cubic concrete elements.

In the previous 1D case study of Section 3.1, the thickness of the interface was simply equal to the thickness of the concrete. On the contrary, 3D volumetric elements are used here for concrete. This is why the first step consists in a sensitivity analysis on the thickness value of the interface part of the macro-element.

#### 3.2.1. Material properties

The steel linear behavior law and the bond law are identical to the ones used for the 1D pull-out model. Choosing a linear behavior for steel is a convenient choice since the steel bar does not yield in the experiment of [13]. For this 3D model, a regularized Mazars' nonlinear damage law is considered for concrete [41,42]. The used properties for concrete are summarized in Table C.1.

#### 3.2.2. Finite elements mesh

The mesh shown in Fig. 14(a) is chosen. Fig. 14(b) shows the boundary conditions detailed on a side view of the mesh cross section shown in Fig. 14(a). A displacement is imposed on an external steel node of the macro-elements. Three macro-elements connect steel and concrete. For this 3D simulation, no discretization has been done in each macro-element. Interface and steel nodes of the macro-elements have the exact same coordinates before the imposed displacement is applied.

#### 3.2.3. Results

Fig. 15 illustrates the reaction curves for a monotonic loading scheme. A sensitivity study is performed with respect to the interface thickness. Fig. 15(b) shows the reaction curves for different thickness values. Changing this thickness slightly changes the initial stiffness of the structure.

Fig. 16 shows a cyclic imposed displacement path and the corresponding reaction curve, where the interface thickness is equal to the steel bar diameter. Boundary conditions are successfully modified between the two configurations shown in Fig. 17 each time the direction of the imposed displacement is changed.

The advantage of the 3D simulation with respect to the 1D one is the capability to visualize the damage field in the concrete volume. Fig. 18 shows the distribution of the concrete damage field at different loading states denoted A1, A2, and A3 (see Fig. 16) of the cyclic application.

### 3.2.4. Discussion

A change of the global stiffness of the specimen due to the change of the interface thickness is expected since changing the thickness of the interface changes its stiffness. At the global level of the structure, this results in a slight change of the global initial stiffness. Among the different tested thickness values for the monotonic case, a value of 12 mm gives the closest curve to the curve resulting from the 1D simulation. This value is equivalent to the diameter of the steel bar. In addition, reaction curves are so close to each other and almost identical for thickness values between 3 mm and 12 mm. The main reason why for the 3D pull-out model no significant influence of the
<figure><img src="image_11.png" /><figcaption>Fig. 9. Monotonic 1D pull-out; imposed displacement (a) and reaction (b).</figcaption></figure>

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

<figure><img src="image_13.png" /><figcaption>Fig. 10. Cyclic 1D pull-out; imposed displacement (a) and reaction (b).</figcaption></figure>

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

<figure><img src="image_15.png" /><figcaption>Fig. 11. Configurations (a) and (b) for the cyclic 1D pull-out model.</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>Fig. 12. Monotonic 1D pull-out; concrete (a) and steel (b) forces.</figcaption></figure>

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

interface thickness is observed on the reaction curve is that the applied load is a longitudinal one.

The concrete damage is localized in a specific zone around the pulled steel bar. Indeed, for the 3D simulation, the nonlinearities are divided into two types: the nonlinearities of the interface model carried

out by the macro-element, and the concrete damage. Hence, since an important part of the total nonlinearities is considered at the interface level, only a small concrete area surrounding the steel is damaged. This result justifies the choice to consider a linear constitutive law for steel and concrete in the 1D pull-out model. In other terms, the 1D pull-out
<figure><img src="image_19.png" /><figcaption>Fig. 13. Monotonic 1D pull-out; concrete (a) and steel (b) nodes displacements.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 14. 3D pull-out test mesh (a) and (b) boundary conditions.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 15. Reaction curves (b) for the monotonic 3D pull-out model.</figcaption></figure>

model can be considered as a sufficiently detailed and representative model of the pull-out test.

### 3.3. Tie rod test: linear concrete behavior model

Tie-rod tests are usually studied in order to understand the influence of the steel–concrete interface on the cracks distribution in concrete. The experimental test of [43] is considered. This test was previously studied in [10,11,13]. The aim of the tie-rod model here is to test the

possibility to reproduce similar results using macro-elements to connect steel and concrete.

The steel bar is tied at one end, and a horizontal displacement is imposed at its other end. It has a diameter of 10 mm, while the concrete has a cross-section of 0.1 × 0.1 m² (see Fig. 19).

Here a linear behavior is considered for concrete, with the aim of focusing the study on the consideration of a linear or non-linear steel–concrete bond behavior. Hence, two interface behavior options are tested: the first one with a perfect steel–concrete bond and the
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 16. Cyclic 3D pull-out; imposed displacement (a) and reaction (b).</figcaption></figure>

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

<figure><img src="image_25.png" /><figcaption>Fig. 17. Configurations (a) and (b) for the cyclic 3D pull-out model.</figcaption></figure>

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

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

<figure><img src="image_28.png" /><figcaption>Fig. 18. Damage fields for different points of the reaction curve of the 3D pull-out cyclic application.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 19. Presentation of the tie-rod test.</figcaption></figure>

second one considering a non-linear interface model. For a perfect bond, the macro-element is not used but steel bar elements are perfectly connected to concrete elements with kinematic relations. On the contrary, when the interface non-linearities are taken into account in the simulation, the macro-element formulation is used.

### 3.3.1. Material properties

The concrete's linear constitutive law is characterized by a Young's modulus of 30.4 GPa and a Poisson's ratio of 0.22.

It is important to note that steel only yields at the end of the tie-rod test. But since the interest here is to study the cracking phase of concrete in order to assess the influence of the steel-concrete interface on the cracking process, there is no need to define an elastic-plastic law for steel. A linear behavior is here used for the steel bar. The steel Young's modulus is equal to 200 GPa.
<figure><img src="image_30.png" /><figcaption>Fig. 20. Bond law.</figcaption></figure>

The interface bond law is identical to the bond law used in [10] when modeling this tie-rod example. This law illustrated in Fig. 20. The different slip-bond stress values indicated in Fig. 20 are linked with linear expressions.

### 3.3.2. Finite elements mesh

A major objective of this linear concrete simulation choice is to compare the results to the those of the same simulation done by [10]. This comparison allows to validate the representativity of the macro-element formulation. In other words, this comparison is done to test the possibility to reproduce similar results by replacing the interface model of [10] with the macro-elements approach. Therefore, the mesh is chosen identical to the one used in [10] for this type of simulations (linear behavior for concrete). It is generated to have 5 × 5 elements in the concrete cross-section, and 58 elements along the tie-rod axis. When modeling a perfect steel-concrete bond, 58 steel elements are used along the contact length (equal to the total concrete length of 1.15 m), and two additional steel elements are linked to the edges of the first 58 elements. When considering a nonlinear model for the interface, 58 macro-elements are used for the steel-concrete contact length. Two steel bar elements are linked to the extremities of the steel part of the macro-elements. The two mesh options are illustrated in Fig. 21.

### 3.3.3. Results

The length of the two steel elements added outside the concrete volume is not indicated in [10] although changing it slightly changes the stiffness of the tie-rod specimen. The length of each one of these two steel elements is here taken equal to 0.05 m. This value is deduced by modeling the tie-rod with a linear behavior for concrete and a perfect steel-concrete bond. The reaction-displacement curve is then plotted and the value of 0.05 m is chosen since it gives the same reaction curve of [10], for the exact simulation option. When macro-elements are used to describe the non-linear behavior of the interface, each macro-element is discretized in two biphasic elements. The macro-element formulation is introduced while keeping the choice of a linear behavior for concrete. Hence, for this case, the nonlinearity is only related to the non-linear bond law (Fig. 20).

The convergence criterions are identical for this simulation and the one of [10]. The only difference between the two simulations is the way the interface is modeled (here the macro-element approach is adopted and in [10] four-node elements connect steel and concrete). This allows to compare the number of iterations needed for convergence and to assess the numerical efficiency of the two interface modeling approaches.

[10] provides the reaction curves when concrete is described by a linear elastic constitutive law using the interface models of [10,13] for the nonlinear interface simulation case. The models of [10,13]

have different formulations but are based on the same principle of introducing bond stresses in a finite element reinforced concrete calculation. So, it is expected to have the same reaction curve using the models of [10,13] (see Fig. 22). Using the macro-element approach gives a close but a slightly different reaction curve. The macro-element formulation incorporates the behavior of the interface zone, which is not the case for the models of [10,13]. For this reason, it is expected not to have the same reaction curve using the macro-element approach and the models of [10,13]. It is important to note that, as shown in Fig. 22, changing the thickness of the interface zone of the macro-element can slightly change the stiffness of the specimen. However, the curves associated to different thickness values are close to each other and almost identical.

The number of iterations needed for the convergence at each time step is given in [10] for the simulation's assumptions of a perfect bond and nonlinear interface models of [10,13]. Choosing the same number of the total time steps, the number of iterations of the different simulations can be compared to the number of iterations using a macro-element interface model. A 5 mm interface thickness is chosen for the macro-element simulation for the comparison (Fig. 23).

### 3.3.4. Discussion

Fig. 23 demonstrates the numerical robustness of the macro-element formulation with respect to the interface models of [10,13]. It can be seen that the number of iterations at each time step is stable and does not exceed 3 iterations. Table 3 sums up the total number of iterations for the different simulations cases.

The ratio of the total number of iterations using the different interface modeling methods to the total number of iterations with a perfect bond is equal to 1.445 using the model of [13] and 1.495 using the model of [10]. For the macro-element simulation, this ratio becomes equal to 1.045 that is to say approximately 33% less.

The reason for which the usage of the macro-elements implies a better numerical robustness is that a local algorithm is used to establish the inner equilibrium at the level of each macro-element. The local algorithm is called each time the stiffness matrix is calculated or updated, and when calculating the internal forces. Here the macro-element local algorithm is called 519 times:

* A first time when calculating the initial stiffness of the structure.

* 100 times at the beginning of each time step when updating the stiffness matrix.

* 209 × 2 times at each iteration when updating the stiffness matrix and when calculating the internal forces to estimate the residue.

The local iterations number is equal to one for the most of the times the local algorithm is called. The maximum observed iterations number is equal to three. The macro-elements are numbered from 1 to 58 in an ascending order with respect to their longitudinal position. Fig. 24 shows the local iterations number for macro-elements 1, 30, and 40.

The local iterations are realized at the level of the macro-elements, which means that they are less costly than the global iterations at the level of the whole structure. For this tie test, the computational time of one local iteration is around 1% of the computational time of one global iteration. Table 4 compares the computational time ratios of the nonlinear interface simulations with respect to a perfect bond simulation.

This example demonstrates the advantage of inner discretization that accelerates the convergence. However, the difference between the computational time of an inner iteration and a local one becomes more important while modeling large scale structures since a global iteration is at the level of the whole structure while a local one is at the level of one macro-element.

### 3.4. Tie rod test: nonlinear concrete behavior model

The tie test of [43] is here modeled by considering the concrete damage.
<figure><img src="image_31.png" /><figcaption>Fig. 21. 3D mesh view (a) and (b) 2D mesh cut-section view-Linear concrete behavior simulations.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Fig. 22. Reaction curves for linear concrete simulations.</figcaption></figure>

<table><caption>Table 3
Total number of iterations.</caption><thead><tr><th>Interface model</th><th>Perfect bond</th><th>Model of [13]</th><th>Model of [10]</th><th>Macro-element</th></tr></thead><tbody><tr><td>Total number of iterations</td><td>200</td><td>289</td><td>299</td><td>209</td></tr></tbody></table>

<table><caption>Table 4
Total computational time ratios.</caption><thead><tr><th>Interface model</th><th>Model of [13]</th><th>Model of [10]</th><th>Macro-element</th></tr></thead><tbody><tr><td>Computational time ratio with respect to a perfect bond simulation</td><td>3.4375</td><td>2.5</td><td>1.426</td></tr></tbody></table>

### 3.4.1. Material properties

A regularized Mazars' damage law [41,42] is used for the nonlinear concrete simulations. Concrete parameters are presented in Table C.1. The steel linear behavior law and the interface bond law are identical to the laws used for the linear concrete simulations.

### 3.4.2. Finite elements mesh

For the nonlinear concrete behavior simulations, the mesh is generated to have 3 × 3 elements in the concrete cross-section, and 35 elements along the tie-rod axis. Each macro-element is discretized into 5 biphasic elements, and a thickness of 2 mm is here used for the interface. An aleatory strain-based damage threshold distribution is considered for concrete elements (Fig. 25). This distribution follows an average Gaussian law with a mean value equal to the ratio of the tensile strength of the concrete to its Young's modulus (value indicated in Table C.1) and a coefficient of variation of 5%. An isotropic correlation is used with a correlation length equal to the size of the concrete

elements. The turning bands method is used to generate the damage threshold distribution [44].

### 3.4.3. Results

Fig. 26 shows the reaction curves for the nonlinear concrete behavior simulations. The numbers one to five indicated in this figure designate the successive apparition of the different concrete cracks for the macro-element simulation.

By comparing the damage distributions of Fig. 27 (for a nonlinear interface behavior simulation and a perfect bond simulation), one can deduce that when introducing a nonlinear behavior at the interface, a smaller number of cracks occurs in concrete, with a higher mean distance between the different cracks.

An average concrete strain is experimentally calculated using the relative displacement on a length equal to 1 m in the center of the tie-rod [43] (see Fig. 28). Fig. 29 shows the evolution of the steel stress at the end of the rod (ratio of the applied force to the steel section)
<figure><img src="image_33.png" /><figcaption>Fig. 23. Iterations number with respect to the pseudotime.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 24. Local iterations number.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 25. Aleatory distribution of the damage threshold in the tie-rod specimen.</figcaption></figure>

as a function of the average concrete strain. The experimental and numerical curves are also compared in Fig. 29.

### 3.4.4. Discussion

When comparing the damage patterns of the two simulations (perfect bond and nonlinear interface using macro-elements) to the experimental cracking pattern, one can deduce that modeling a nonlinear

<figure><img src="image_36.png" /><figcaption>Fig. 26. Reaction curves for the nonlinear concrete behavior simulations.</figcaption></figure>

interface behavior gives a better estimation of the experimental cracking features. The mean value of the spacing between cracks is 16 cm experimentally and 18 cm for the macro-element simulation, while it is about 11.2 cm for the perfect bond simulation. From the experimental side, the total number of cracks is 5. For the macro-element simulation 5 cracks are developed while for the perfect bond simulation the development of total of 10 cracks is observed (Fig. 30).

The use of a non-linear bond model allows to observe steel stress drops at the cracking moment (Fig. 29). The amplitude of these drops is related, according to [43], to the experimental device. Since the experimental device is not modeled here, it is not of interest to compare the experimental-numerical drop amplitudes. The numerical curve is judged to be a satisfactory representation of the experimental curve (Fig. 29).

The main result of this tie-rod study revealed that a better representation of the experimental concrete cracking pattern (in terms of number of cracks and cracks spacing) is done when taking into account the nonlinear behavior of the steel-concrete interface. In addition, this application showed the numerical robustness of the macro-element formulation and demonstrated the capability of the proposed approach to reproduce the experimental behavior of the interface.

### 3.5. Beam-end test

Beam-end tests can be adopted as an alternative interface characterization setup to the pull-out tests [45–48]. Unlike the pull-out test where the concrete cube is in compression, the beam-end test represents the bond condition in a beam. Concrete around the rebar bond area is subjected to a constant shear force and a bending moment. Moreover, the concrete cover in a pull-out test is usually higher than the minimum cover in practical constructions (except the case where the bar diameter is relatively small). The beam-end test of [48] is modeled in the following.

The specimen dimensions and setup are shown in Fig. 31. The pulled steel bar of a 16 mm diameter is in direct contact with concrete for a specific distance A2A3 (see Fig. 31(b)) of 80 mm. This distance is here equal to 7 times the pulled steel bar diameter.

#### 3.5.1. Material properties

Steel bars are supposed to have an elastic behavior with a Young's modulus of 200 GPa, which is compatible with the experimental observations since steel does not yield in the studied beam-end test. A regularized Mazars' law is used for concrete. Table C.1 summarizes the concrete material parameters.
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /><figcaption>Fig. 27. Damage patterns for a nonlinear interface model using the macro-element (a) and a perfect steel-concrete bond (b).</figcaption></figure>

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

$$ \epsilon = \frac{u_B - u_A}{1 m} $$

<figure><img src="image_40.png" /><figcaption>Fig. 29. Experimental and numerical steel stress-concrete strain curves.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. 30. (a) Experimental cracks pattern, (b) damage field distribution for the macro-element simulation, (c) and for the perfect bond simulation.</figcaption></figure>

The bond law of [40] Appendix B is used. It is interesting here to estimate the bond law parameters $\tau_1$, $g_1$, and $g_3$ before comparing the simulation results with the experimental ones. Different empirical expressions are recommended in the literature for these three parameters. Table 5 presents different empirical expressions for the bond stress $\tau_1$ which is the maximum bond stress of the initial global envelop curve of the bond law.

One single expression is given in [49–51], while two distinct expressions are given in [13] for pull-out and splitting failure types of the

interface. In [52], several expressions are given as follows (Eq. (41)):

$$ \left\{ \begin{array}{ll} 2.5\sqrt{f_c} & : \text{ pull-out failure with a good bond condition} \\ 1.25\sqrt{f_c} & : \text{ pull-out failure with all other bond conditions} \\ 7.0\left(\frac{f_c}{20}\right)^{0.25} & : \text{ splitting failure with a good bond condition} \\ & \text{forunconfined concrete} \\ 8.0\left(\frac{f_c}{20}\right)^{0.25} & : \text{ splitting failure with a good bond condition} \\ & \text{forconcrete confinement due to stirrups} \\ 5.0\left(\frac{f_c}{20}\right)^{0.25} & : \text{ splitting failure with all other bond conditions} \\ & \text{forunconfined concrete} \\ 5.5\left(\frac{f_c}{20}\right)^{0.25} & : \text{ splitting failure with all other bond conditions} \\ & \text{forconcrete confinement due to stirrups} \end{array} \right. \qquad (41) $$

$f_c$ and $f'_c$ are respectively concrete cubic and cylindrical compressive strength values. c and d subscripts in Table 5 refer respectively to the concrete cover and the steel bar diameter.

Some proposals are given for the slip $g_1$, which is the slip value corresponding to the bond stress $\tau_1$ (Table 6).

It is recommended to attribute the value of the clear spacing $c_{clear}$ between the steel ribs to the bond law parameter $g_3$ in [50,53]. [52] recommends the same value for the pull-out failure, a value of $1.2g_1$ for a splitting failure under unconfined concrete conditions, and $0.5c_{clear}$ for a splitting failure with a concrete confinement due to stirrups.

For the calculation of $\tau_1$ and $g_1$, the recommendations of [13] can be a convenient choice since distinct values are given for the two failure mechanisms, and yet these recommendations are simple compared to the recommendations of [52]. In contrast, no indications are given in [13] to create a dependence between the slip $g_3$ and the failure mechanism, which should be the case according to [54] (see Fig. 32). Thus, the recommendations of [52] are adopted. These choices give an identification methodology for the parameters of the used bond law, which is simple yet dependant of the failure mechanism. Predicting the mechanism that is more likely to occur can be done using the criterion of [13] (Eq. (42)). For the studied beam-end test, a splitting failure is expected. Table 7 sums up the bond law parameters.

$$ \left\{ \begin{array}{l} \frac{c}{d} < 0.39 \frac{f_c}{f_i} - 0.24: \text{splitting failure} \\ \frac{c}{d} \ge 0.39 \frac{f_c}{f_i} - 0.24: \text{pull-out failure} \end{array} \right. \qquad (42) $$

### 3.5.2. Finite elements mesh

Fig. 33(a) shows the concrete mesh. Fig. 33(b) details the configuration of the steel bar elements and the macro-elements, where macro-elements connect steel and concrete. No inner discretization is done here for the macro-elements; each macro-element is composed of one biphasic element. For the sake of simplicity, the part A3A4 of the pulled steel bar (see Fig. 31) is not explicitly modeled in the numerical simulation. It is important to note here that the interface
<figure><img src="image_42.png" /><figcaption>Fig. 31. Beam-end 3D specimen (a) and longitudinal section (b).</figcaption></figure>

<table><caption>**Table 5**
Empirical expressions for $\tau_1$.</caption><tr><td>τ<sub>1</sub> (MPa)</td><td>f<sub>t</sub>(1.53 ε<sub>d</sub> + 0.36): splitting failure0.6 f<sub>c</sub>: pull-out failure</td><td>2.305 f<sub>c</sub><sup>0.555</sup></td><td>1.163 f<sub>c</sub><sup>0.75</sup></td><td>2.57 f<sub>c</sub><sup>0.5</sup></td><td>Eq. (41)</td></tr><tr><td>Reference</td><td>[13]</td><td>[49]</td><td>[50]</td><td>[51]</td><td>[52]</td></tr></table>

<table><caption>**Table 6**
Empirical expressions for $g_1$.</caption><tr><td>g<sub>1</sub> (mm)</td><td>0.17 ε<sub>d</sub>: splitting failure 1: pull-out failure</td><td>0.07d</td><td>1: pull-out failure with a good bond condition 1.8: pull-out failure with other bond conditions Value not specified for splitting failure</td></tr><tr><td>Reference</td><td>[13]</td><td>[50]</td><td>[52]</td></tr></table>

<figure><img src="image_43.png" /><figcaption>Fig. 32. Bond law shape for different failure mechanisms.</figcaption></figure>

<table><caption>**Table 7**
Bond law parameters for the beam-end model.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>τ<sub>1</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>10.7</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>0.425</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>5</td><td>mm</td></tr></tbody></table>

and the steel nodes of the macro-elements have the same position at the initial configuration even though distinct nodes positions are illustrated in Fig. 33(b). A thickness of 8 mm is here associated to the interface zone of the macro-elements.

Concrete and steel meshes are chosen to be coincident. The displacements of the steel bars beside the pulled one are equal to the displacements of the coincident concrete nodes.

The boundary conditions indicated in Fig. 31 are applied. One steel bar element is used between the steel nodes A1 and A2. The displacements of node A1 and the totality of the macro-elements nodes are blocked in the normal directions y and z. The longitudinal displacements of the interface part of the macro-elements is imposed to be equal to the displacements of concrete nodes, while the steel part is free to slide with respect to the interface.

An additional boundary condition is imposed to the concrete node P indicated in Fig. 33(a). The normal displacement of this node in the y direction is blocked. This boundary condition is added to prevent the displacement of the whole concrete specimen in the y direction.

### 3.5.3. Results

Fig. 34 shows the evolution of the bond stresses with respect to the steel-concrete slip. This evolution is compared to the experimental range. The slip between steel and concrete is calculated at the longitudinal position of the node A3.

The concrete damage field at the time step for which the slip value is equal to 5 mm is shown in Fig. 35.

### 3.5.4. Discussion

Fig. 34 shows that the adopted approach to identify the bond law parameters gives a satisfactory agreement with the experimental results. Since the damage reaches the concrete face (Fig. 35), the damage field demonstrates that a splitting failure mode is numerically well reproduced. Fig. 36 compares the experimental crack pattern to the numerical damage field distribution. It shows that the position of point A where surface cracks with different directions are initiated is numerically well estimated. Experimental longitudinal and 45 degrees inclined cracks are represented in the numerical damage field. However, the numerical damage distribution does not represent the totality of the experimental cracks. This is due to the fact that in the numerical simulation the nonlinearities of the beam-end test are carried out at the
<figure><img src="image_44.png" /><figcaption>Fig. 33. Beam-end test concrete mesh (a) and steel bars configuration (b).</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Fig. 34. Numerical and experimental curves for the bond stress values with respect to the steel-concrete slip.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. 36. Damage field in the concrete beam-end specimen compared to the surface experimental cracks pattern.</figcaption></figure>

<caption>Table 8</caption>

Steel parameters: four-point flexural beam test.

<table><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>

<figure><img src="image_47.png" /><figcaption>Fig. 35. Damage field in the concrete beam-end specimen.</figcaption></figure>

level of the interface by the macro-elements and in concrete elements as a damage field. In other terms, the concrete damage does not represent the whole nonlinearities.

It is convenient to note here that a sensitivity study on the choice of the nonlinear damage concrete model may be interesting to reduce the disparity between the numerical and the experimental results. This

study focuses anyway on the choice of the bond law parameters, that define its shape.

## 3.6. Four-point flexural beam test

The macro-element formulation is used to model a four-point reinforced concrete beam bending test carried out by [55]. The considered test is described in Fig. 37.

### 3.6.1. Material properties

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

A revised Mazars' law is considered for concrete [56], with the parameters of Table 9. A Hillerborg regularization is applied to this law [42].

The bond law of [40] (see Appendix B) is used. Similarly to the numerical application of Section 3.5, the bond law parameters are identified as follows:

* the prediction of the interface failure type is done using Eq. (42). A splitting failure is predicted;

* $\tau_1$ and $g_1$ are estimated by following the recommendations of [13];

* $g_3$ is estimated by following the recommendations of [52];

Table 10 sums up the bond law parameters.
<figure><img src="image_48.png" /><figcaption>Fig. 37. Configuration of the four-point beam test of [55].</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>Fig. 38. Mesh an boundary conditions of the four-point beam test.</figcaption></figure>

<table><caption>**Table 9**

Concrete nonlinear parameters: four-point bending test.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>f<sub>i</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>&epsilon;<sub>10</sub></td><td>Revised Mazars' model input</td><td>f<sub>i</sub>/E<sub>c</sub> = 1.1429 &times; 10<sup>-4</sup></td><td>-</td></tr><tr><td>&epsilon;<sub>03</sub></td><td>Revised Mazars' model input</td><td>f<sub>i</sub>/E<sub>c</sub> = 1.9 &times; 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>2</sub></td><td>Revised Mazars' model input</td><td>1.2</td><td>-</td></tr><tr><td>B<sub>2</sub></td><td>Revised Mazars' model input</td><td>400</td><td>-</td></tr><tr><td>&beta;</td><td>Revised Mazars' model input</td><td>1.06</td><td>-</td></tr></tbody></table>

<table><caption>**Table 10**

Bond law parameters for the beam model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>&tau;<sub>1</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>13.392</td><td>MPa</td></tr><tr><td>g<sub>1</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>0.425</td><td>mm</td></tr><tr><td>g<sub>3</sub></td><td>Input parameter of the bond law (Fig. B.1(a))</td><td>0.51</td><td>mm</td></tr></tbody></table>

<figure><img src="image_50.png" /><figcaption>Fig. 39. Reaction curves of the four-point beam test.</figcaption></figure>

### 3.6.2. Finite elements mesh

A plane stress configuration is chosen for the beam test model. The considered mesh and the boundary conditions are detailed in Fig. 38. Each macro-element of Fig. 38 is composed of one biphasic element. The steel part of the macro-elements has an equivalent section of the two 16 mm steel bars of Fig. 37. A thickness of 1.6 mm is associated to the interface part of the macro-elements. This thickness is equal to 0.1d, where d is the steel bars diameter of 16 mm.

### 3.6.3. Results

Fig. 39 presents the numerical reaction curves with a perfect bond state and a nonlinear interface behavior, compared to the experimental curve. It is shown in this figure that the consideration of the nonlinear interface behavior does not provide a considerable difference of the reaction curve of the four-point bending test. Nevertheless, taking into account a nonlinear interface behavior in the numerical simulation provides a better representation of the experimental cracking pattern, as shown in Fig. 40, especially in the zoomed areas of this figure.

### 3.6.4. Discussion

The consideration of the interface behavior helps localizing the damage around the steel reinforcement, and gives a better representation of the experimental cracking pattern.

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 Fig. 41 that presents the steel stress and the slip values for the resolution time steps.

## 3.7. Three-point flexural beam test

In order to evaluate the performance of the macro-element model in the context of a structural reinforced concrete element, the three-point beam bending test of [1] is modeled.

The studied reinforced concrete beam has a total length of 1.5 m and a span of 1.4 m. Its cross-section measures 0.15 m in width and 0.22 m in height. In the tension zone of the beam, two longitudinal steel
<figure><img src="image_51.png" /><figcaption>Fig. 40. Damage distribution of the four-point beam test with perfect and imperfect bond behaviors compared to the experimental cracking pattern.</figcaption></figure>

<table><caption>**Table 11**
Bond law parameters: three-point flexural beam test.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>τ<sub>1</sub></td><td>12.59</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>

<table><caption>**Table 12**
Steel parameters: three-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>o</sub></td><td>Young's modulus</td><td>200</td><td>GPa</td></tr><tr><td>f<sub>y</sub></td><td>Limit of elasticity</td><td>500</td><td>MPa</td></tr><tr><td>E<sub>h</sub></td><td>Hardening slope</td><td>4000</td><td>MPa</td></tr></tbody></table>

<figure><img src="image_52.png" /><figcaption>Fig. 41. Local cracking description of the four-point beam test.</figcaption></figure>

bars with a diameter of 14 mm each are used. These bars are anchored to the concrete at both ends. Additionally, two longitudinal bars with a diameter of 8 mm each are placed in the compressed zone of the beam. Along beam length, stirrups with a diameter of 6 mm connect the longitudinal bars. The geometric configuration of the beam test conducted by [57] is illustrated in Fig. 42.

### 3.7.1. Material properties

The three-point bending test is here modeled by assigning nonlinear behaviors to the steel, the interface, and the concrete.

A regularized Mazar's nonlinear damage law [41,42] is considered for concrete. Table C.1 sums up the adopted material properties for concrete.

The bond law of [2] is used, with the values indicated in Table 11 of its defining parameters.

An elastic plastic behavior with a linear strain hardening is assigned to the steel with the indicated parameters of Table 12.

### 3.7.2. Finite elements mesh

Because of the beam symmetry, a quarter of its geometry is modeled. Three mesh sizes are tested, as shown in Fig. 43.

Macro-elements are used to represent the longitudinal bottom steel with its interface behavior with respect to concrete. Macro-elements and concrete meshes are coincident. The interface part of the macro-elements has the same displacement as concrete in all directions. The steel part has the same displacements as concrete in the normal directions (via kinematic relationships) but not in its longitudinal direction. Along this longitudinal direction, the steel can slide relatively to the interface, the interface being perfectly attached to the concrete. Each macro-element is discretized into 10 biphasic elements. Fig. 44 shows the configuration of the steel bars and the macro-elements.

Fig. 45(a) shows the boundary conditions applied to the concrete volume. The conditions imposed on surfaces S1 and S2 of this figure establish the double symmetry conditions of the beam. The boundary
<figure><img src="image_53.png" /><figcaption>Fig. 42. Configuration of the three-point beam test of [1].</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Fig. 43. Different mesh sizes tested: coarse mesh (a), intermediate mesh (b) and fine mesh (c).</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Fig. 44. Configuration of the steel elements and the macro-elements of the three-point beam bending model within the 3D concrete volume.</figcaption></figure>

conditions applied to the steel elements and the macro-elements are illustrated in Fig. 45(b).

Anchoring the bottom steel bar extremity to the concrete is achieved through a kinematic relation that requires the external steel node to exhibit the same displacement as the corresponding concrete node in all directions (Fig. 45(b)). The steel and interface nodes belonging to surfaces S1 and S2 of Fig. 45(a) have the same boundary conditions applied to these surfaces.

It is important to highlight that the macro-elements represent the bottom steel behavior with its interface with respect to concrete. The top longitudinal steel and the stirrups are assumed to have a perfect bond to concrete. This assumption is based on the fact that concrete cracking is initiated in its tension zone around the bottom steel reinforcement. For this reason, macro-elements are used to model this bottom steel and its interaction with concrete. A thickness of 0.08d is associated to the interface part of the macro-elements, d being the diameter of the bottom longitudinal reinforcement.

### 3.7.3. Results

Fig. 46 shows the numerical reaction curves, using the intermediate mesh size, compared to the experimental reaction curve. Different values of the interface thickness for macro-elements are tested. This figure

validates the macro-element model with respect to the experimental behavior of the studied structure. In addition, a slight sensitivity is observed with respect to the interface thickness parameter. Fig. 47 shows the reaction curves due to the different mesh sizes. The curves associated the intermediate and fine mesh sizes are very close to each other. The curve of the coarse mesh is slightly stiffer. This figure demonstrates that, using the macro-element model, a representative structural behavior can result from the numerical simulation without the need to use a significantly refined mesh.

### 3.7.4. Discussion

The advantage of the nonlinear interface behavior modeling is that it provides a detailed description of the local behavior at the interface between steel and concrete. Indeed, it is possible to have an access to the local values of slip and bond stresses. Damage occurs around the beam mid-span, where the bending moment has its maximal value. The damage is initiated in the tension zone and then propagates from the lower to the upper fibers of the beam.

Figs. 48, 49, and 50 show that the cracks initiation in concrete around the steel reinforcement is characterized by:

* a of steel stress concentration at the level of the cracked concrete element;

* a discontinuity of the slip value that links the slips the ends of the cracked element.

The same observations concerning the local characterization of the concrete cracking are captured using the different mesh sizes.

## 4. Conclusive remarks and perspectives

The present paper proposes a multiscale steel–concrete bond model within the context of the finite element method. The proposed model consists in defining a macro-element formulation defined at two complementary scales: a global scale and a local one. At the global scale (the scale of the studied reinforced concrete structure), the macro-element is a four-node finite element. At its local scale, the macro-element is an assembly of a set of biphasic elements. Each biphasic element has two parallel three-node bar elements: one bar element
<figure><img src="image_56.png" /><figcaption>Fig. 45. Boundary conditions applied to concrete (a), steel elements and the macro-elements (b).</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>Fig. 46. Reaction curves with perfect and imperfect bond conditions compared with the experimental curve.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Fig. 47. Reaction curves associated to different mesh sizes.</figcaption></figure>

represents a concrete interface zone and the other one represents the steel. Moreover, bond laws connect the nodes of the two bar elements.

In addition to the global iterative resolution at the global scale of the studied structure, an inner resolution is done at the level of each macro-element. The internal degrees of freedom are condensed on the four macro-element external nodes. It has been shown in this paper that the implementation of this macro-element formulation has no impact on the main architecture of a classical finite element code.

Numerical examples of 1D and 3D pull-out models are performed using the macro-element formulation. These models validate the macro-element formulation that is able to reproduce the interface behavior for monotonic and cyclic load configurations.

In a second example, macro-elements are used to model a tie-rod experiment. To demonstrate the numerical robustness of the macro-element model, the number of iterations required for convergence is compared to the number of iterations for the same simulation with other interface models of the literature. It has been shown how well the macro-element formulation is able to reduce the number of iterations at the main Newton-Raphson algorithm. Taking into account the interface behavior in the simulation leads to a better representation of the experimental cracking pattern. By comparing the numerical and the experimental curves of the steel stress with respect to concrete strain, one can deduce that the macro-element can reproduce the tie-rod experimental behavior.

The interface model is then used to reproduce the experimental behavior of a beam test. Indeed, beam tests represent a better structural evaluation of the interface behavior with respect to classical pullout tests. The bond law parameters are identified based on empirical expressions proposed in the literature. For this test, an experimental splitting failure mode of the interface is reproduced with the macro-element model.

Some suggestions to upgrade the current version of the macro-element formulation are the following:

* The consideration of normal stresses between the two zones is not possible by using bar elements in the macro element. However, it can be done by replacing the three-node bar elements with beam elements.

* In the current work, a linear behavior is assigned to the interface zone of the macro-element. However, non-linear behavior can also be assigned to this zone which is interesting to investigate the effect of this type of non-linearity on the macro-element behavior.

* In the current work, since the steel reinforcement is represented by one-dimensional truss finite elements, bending and shear effects are not considered. Considering beam steel elements in future works makes it possible to take into account these effects. In addition, simple linear bar geometries are considered. Inspired by the work of [58], the macro-element philosophy can be applied in future works to model complex steel reinforcement shapes. In
<figure><img src="image_59.png" /><figcaption>Fig. 48. Local observations of the slip distribution and the steel stress along the longitudinal position due to concrete damage: coarse mesh.</figcaption></figure>

this case, two rod elements can be used to model steel and inter-
face zones with bond stresses in between, within the framework
of macro-elements.

Finally, additional numerical tests using the macro-element can be
interesting to perform. These applications can give a better under-
standing on the effect of the interface zone thickness on the structural
response of reinforced concrete elements. Applications for large-scaled
structures subjected to both static and dynamic loads will be addressed
in a future paper.

CRediT authorship contribution statement

**Maryam Trad:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. **Ibrahim Bitar:** Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation. **Stéphane Grange:** Conceptualization, Data
<figure><img src="image_60.png" /><figcaption>Fig. 49. Local observations of the slip distribution and the steel stress along the longitudinal position due to concrete damage: intermediate mesh.</figcaption></figure>

curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, 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.

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

## Appendix A. Double Lagrange multipliers methodology

The double Lagrange multipliers methodology introduces two Lagrange multipliers for each Dirichlet condition, and Eqs. (11) and (13) are combined as follows:

$$ \begin{bmatrix} F_{int}(U) + L^T \lambda^1 + L^T \lambda^2 \\ LU - \alpha \lambda^1 + \alpha \lambda^2 \\ LU + \alpha \lambda^1 - \alpha \lambda^2 \end{bmatrix} = \begin{bmatrix} F_{ext} \\ U_d \\ U_d \end{bmatrix} \qquad (A.1) $$

Eq. (A.1) can be written in a more condensed form as:

$$ F(U_{tot}) = F_{tot} \qquad (A.2) $$

Where:

$$ U_{tot} = [ U \quad \lambda^1 \quad \lambda^2 ]^T, F_{tot} = [ F_{ext} \quad U_d \quad U_d ]^T \qquad (A.3) $$

The linearization of the term $F(U_{tot})$ of Eq. (A.2) gives:

$$ dF(U_{tot}) = k_{tot}dU_{tot} \qquad (A.4) $$
<figure><img src="image_61.png" /><figcaption>Fig. 50. Local observations of the slip distribution and the steel stress along the longitudinal position due to concrete damage: fine mesh.</figcaption></figure>

Where:

$$k_{tot} = \begin{bmatrix} k & L^T & L^T \\ L & -\alpha I & \alpha I \\ L & \alpha I & -\alpha I \end{bmatrix} \qquad (A.5)$$

I is the identity matrix. The value of α is chosen to ensure a good conditioning of the matrix $k_{tot}$.

## Appendix B. Cyclic bond law

The bond law of [40] is here described. For the monotonic initial envelope curve of the law, the bond stress is assumed to be the sum of two types of stresses: the friction stress $τ_f$ and the bearing stress $τ_b$. The maximum stress value $τ_1$ is estimated by considering a contribution of

25% of the friction stress and 75% of the bearing stress. This envelope curve is reduced due to the unloading/reloading cycles (see Fig. B.1).

Eq. (B.1) details the constitutive equations of the initial monotonic envelope curve.

$$ \left\{ \begin{array}{ll} \tau = 4 \frac{\tau_1}{g} g & \text{for } g \le 0.1 g_1 \\ \tau = \tau_1 \left[ 1 - 0.6 \left( \frac{g-g_1}{0.9 g_1} \right)^4 \right] & \text{for } 0.1 g_1 < g \le g_1 \\ \tau = \tau_1 & \text{for } g_1 < g \le g_2 \\ \tau = \tau_1 \left[ 1 - 0.75 \frac{g-1.1 g_1}{g_3-g_1} \right] & \text{for } g_2 < g \le g_3 \\ \tau = 0.25 \tau_1 & \text{for } g_3 < g \end{array} \right. \qquad (B.1) $$

For cyclic loads, two damage variables are defined: $d_b$ and $d_f$. The two variables are respectively associated to the stresses $τ_b$ and $τ_f$. Due to the loading/unloading cycles, reduced curves $τ_{b,r}$ and $τ_{f,r}$ are defined
<figure><img src="image_62.png" /><figcaption>Fig. B.1. Bond law proposed in [40]: monotonic envelope curve (a); cyclic response (b).</figcaption></figure>

<table><caption>Table C.1
Concrete parameters for the beam-end model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Section 3.2</th><th>Section 3.4</th><th>Section 3.5</th><th>Section 3.7</th><th>Units</th></tr></thead><tbody><tr><td>E<sub>c</sub></td><td>Young's modulus</td><td>28</td><td>30.4</td><td>28</td><td>28</td><td>GPa</td></tr><tr><td>ν<sub>c</sub></td><td>Poisson's ratio</td><td>0.22</td><td>0.22</td><td>0.2</td><td>0.22</td><td>-</td></tr><tr><td>f<sub>t</sub></td><td>Tensile strength</td><td>3.12</td><td>2.6</td><td>2.565</td><td>2.7</td><td>MPa</td></tr><tr><td>ε<sub>d0</sub></td><td>Damage threshold</td><td>L<sub>e</sub>/ε<sub>t</sub> = 1.1143 × 10<sup>-4</sup></td><td>L<sub>e</sub>/ε<sub>t</sub> = 8.5526 × 10<sup>-5</sup></td><td>L<sub>e</sub>/ε<sub>t</sub> = 9.1607 × 10<sup>-5</sup></td><td>L<sub>e</sub>/ε<sub>t</sub> = 9.6429 × 10<sup>-5</sup></td><td>-</td></tr><tr><td>A<sub>t</sub></td><td>Local Mazars' model input</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>B<sub>t</sub></td><td>Mazars' model input</td><td>Regularized</td><td>Regularized</td><td>Regularized</td><td>-</td><td>-</td></tr><tr><td>A<sub>s</sub></td><td>Mazars' model input</td><td>1.1</td><td>1.2</td><td>1.3</td><td>1.2</td><td>-</td></tr><tr><td>B<sub>s</sub></td><td>Mazars' model input</td><td>700</td><td>700</td><td>650</td><td>700</td><td>-</td></tr><tr><td>β</td><td>Mazars' model input</td><td>1.06</td><td>1.06</td><td>1.06</td><td>1.06</td><td>-</td></tr><tr><td>G<sub>f</sub></td><td>Fracture energy</td><td>150</td><td>150</td><td>150</td><td>150</td><td>N/m</td></tr></tbody></table>

as:

$$ \begin{cases} \tau_{b,r} = (1 - d_b) \tau_b \\ \tau_{f,r} = (1 - d_f) \tau_f \end{cases} \qquad (B.2) $$

Where:

$$ \left\{ \begin{aligned} d_b &= 1 - e^{2.5 \left( \frac{g_{max}}{g_3} \right)^{0.8}} \\ d_f &= \frac{g_{max}^{+} + g_{max}^{-}}{g_3} \left( 1 - e^{-0.4 \left( \frac{g_{max}}{g_3} \right)^{0.75}} \right) \end{aligned} \right. \qquad (B.3) $$

The term $g_{max}$ is a combination of the maximum positive and negative slips $g_{max}^{+}$ and $g_{max}^{-}$.

$$ g_{max} = 0.75 \max(g_{max}^{+}, g_{max}^{-}) + 0.125(g_{max}^{+} + g_{max}^{-}) \qquad (B.4) $$

The stress value $\tau_{rev}$ indicated in Fig. B.1 is calculated as follows:

$$ \tau_{rev} = \frac{\max(g_{max}^{+}, g_{max}^{-})}{g_1} \tau_{f,r} \quad \text{with} \quad \frac{\max(g_{max}^{+}, g_{max}^{-})}{g_1} \le 1 \qquad (B.5) $$

Unloading and reloading stiffnesses are equal to the initial stiffness of the law: $k_{in} = 4 \frac{\tau_{f,r}}{g_1}$.

## Appendix C. Regularized Mazars’ constitutive law parameters

The local Mazars’ model input $A_t$ is not used since the law is regularized.

## References

[1] Phan TS, Rossi P, Tailhan JL. Numerical modelling of the concrete/rebar bond. Cem Concr Compos 2015;59:1–9. http://dx.doi.org/10.1016/j.cemcompcon. 2015.02.003.

[2] Aguilera MH. Un modèle global homogénéisé pour la fissuration des voiles en béton armé sous chargements sismiques. École Centrale de Nantes; 2016.

[3] Ngo D, Scordelis AC. Finite element analysis of reinforced concrete beams. ACI J 1967;6(3):152–63. http://dx.doi.org/10.14359/7551.

[4] Nilson AH. Nonlinear analysis of reinforced concrete by the finite element method. J Proc 1968;65(9):757–66. http://dx.doi.org/10.14359/7510.

[5] Reinhardt HW, Blaauwendraad J, Vos E. Prediction of bond between steel and concrete by numerical analysis. Matériaux Constr 1984;17(4):311–20. http://dx.doi.org/10.1007/BF02479089.

[6] Giry C. Modélisation objective de la localisation des déformations et de la fissuration dans les structures en béton armé. Université de Grenoble; 2011.

[7] Ramirez ND. Etude de la liaison acier-béton: de la modélisation du phénomène à la formulation d'un élément fini enrichi 'béton armé'. École normale supérieure de Cachan-ENS Cachan; 2005.

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

[9] Abbas M, Bary B, Jason L. 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.

[10] Mang C. Modélisation de la liaison acier-béton dans le calcul de structures en béton armé. Université Paris Ouest; 2015.

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

[12] Raous M, Karray MA. Model coupling friction and adhesion for steel-concrete interfaces. Int J Comput Appl Technol 2009;34(1):42–51. http://dx.doi.org/10. 1504/IJCAT.2009.022701.

[13] Casanova A. Prise en compte de la liaison acier-béton pour le calcul de structures industrielles. École normale supérieure de Cachan-ENS Cachan; 2012.

[14] Gutiérrez R, Stempniewski L, Fleming W. Modelling steel-concrete interaction using the extended finite element method. Obras y Proyectos 2018;24:6–12. http://dx.doi.org/10.4067/s0718-28132018000200006.

[15] Sellier A, Millard A. 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 2019;213:182–96. http://dx.doi.org/10.1016/j.engfracmech.2019.04.008.

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

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

[18] Khebizi M, Guenfoud H, Guenfoud M. Modélisation des poutres en béton armé par des éléments Multicouches. 2014.

[19] Christelle C. Formulation d'un modèle homogénéisé de plaque en béton armé pour des applications sismiques. Université Pierre et Marie Curie-Paris VI; 2013.

[20] Capdevielle S. Introduction du gauchissement dans les éléments finis multifibres pour la modélisation non linéaire des structures en béton armé. Université Grenoble Alpes; 2016.

[21] Yu JC, Wang JT, Pan JW, Guo N, Zhang CH. A dynamic FEM-DEM multiscale modeling approach for concrete structures. Eng Fract Mech 2023;278:109031. http://dx.doi.org/10.1016/j.engfracmech.2022.109031.

[22] Petrolo M, Carrera E, Cinefra M, Zappino E. Finite element analysis of structures through unified formulation. John Wiley & Sons; 2014.

[23] Shen J, Pagani A, Arruda MRT, Carrera E. Exact component-wise solutions for 3D free vibration and stress analysis of hybrid steel-concrete composite beams. Thin-Walled Struct 2022;174:109094. http://dx.doi.org/10.1016/j.tws. 2022.109094.
[24] Shen J, Arruda MT, Pagani A. 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.

[25] Sahyouni A, Grange S, Briançon L, Burtin P, Racinais J, Prunier F, et al. Développement d'un macro-élément pour modéliser un massif renforcé par des inclusions rigides -étape 1, essai de chargement. 11èmes journées nationales de géotechnique et de géologie de l'ingénieur. 2022.

[26] Jr Bitencourt, Luís AG, Manzoli OL, Trindade YT, Rodrigues EA, Dias-da Costa D. Modeling reinforced concrete structures using coupling finite elements for discrete representation of reinforcements. Finite Elem Anal Des 2018;149:32-44. http://dx.doi.org/10.1016/j.finel.2018.06.004.

[27] Verpeaux P, Charras T. Multiplicateur de Lagrange, Condensation Statique et Conditions Unilatérales. 10ème colloque national en calcul des structures. 2011.

[28] Eligehausen R, Popov EP, Bertero VV. Local bond stress-slip relationships of deformed bars under generalized excitations 1982. http://dx.doi.org/10.18419/opus-415.

[29] Harajli MH. Development/splice strength of reinforcing bars embedded in plain and fiber reinforced concrete. ACI Mater J 1994;91(5):511-20. http://dx.doi.org/10.14359/4163.

[30] Desnerck P, De Schutter G, Taerwe L. 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.

[31] Yankelevsky DZ. A two-phase one dimensional model for steel-concrete interaction. Comput Struct 1997;65(6):781-94. http://dx.doi.org/10.1016/S0045-7949(97)00075-8.

[32] Yankelevsky DZ, Adin MA, Farhey DN. Mathematical model for bond-slip behavior under cyclic loading. Struct J 1992;89(6):692-8. http://dx.doi.org/10.14359/4143.

[33] Kwak HG, Kim SP. Bond-slip behavior under monotonic uniaxial loads. Eng Struct 2001;23(3):298-309. http://dx.doi.org/10.1016/S0141-0296(00)00008-0.

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

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

[36] Haskett M, Oehlers DJ, Mohamed Ali MS. Local and global bond characteristics of steel reinforcing bars. Eng Struct 2008;30(2):376-83. http://dx.doi.org/10.1016/j.engstruct.2007.04.007.

[37] Verderame GM, De Carlo G, Ricci P, Fabbrocino G. Cyclic bond behaviour of plain bars. Part II: Analytical investigation. Constr Build Mater 2009;23(12):3512-22. http://dx.doi.org/10.1016/j.conbuildmat.2009.07.001.

[38] Richard B, Rastiello G, Giry C, Riccardi F, Paredes R, Zafati E, et al. CastLab: an object-oriented finite element toolbox within the Matlab environment for educational and research purposes in computational solid mechanics. Adv Eng Softw 2019;128:136-51. http://dx.doi.org/10.1016/j.advengsoft.2018.08.016.

[39] RILEM. Essais portant sur l'adhérence des armatures du béton :2 : Essai par traction. Mater Struct 1970;3:175-8. http://dx.doi.org/10.1007/BF02478968.

[40] Murcia-Delso J, Stavridis A, Shing B. 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.

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

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

[43] Farra B. Influence de la résistance du béton et de son adhérence avec l'armature sur la fissuration. École polytechnique fédérale de Laussanne; 1995.

[44] Mantoglou A, Wilson JL. The turning bands method for simulation of random fields using line generation by a spectral method. Water Resour Res 1982;18(5):1379-94. http://dx.doi.org/10.1029/WR018i005p01379.

[45] Cairns J, Plizzari GA. Towards a harmonised European bond test. Mater Struct 2003;36:498-506. http://dx.doi.org/10.1007/BF02480826.

[46] ASTM A944-10. Standard test method for comparing bond strength of steel reinforcing bars to concrete using beam-end specimens. ASTM International; 2010, http://dx.doi.org/10.1520/A0944-10R15.

[47] Trujillo PB, Jolin M, Massicotte B, Bissonnette B. Bond strength of reinforcing bars with varying encapsulation qualities. ACI Struct J 2018;115(6):1707-17. http://dx.doi.org/10.14359/51702415.

[48] Rex J, Sharma A, Hofmann J. A new test specimen to investigate the bond behavior of post-installed reinforcing bars. Struct Concr 2019;20(2):583-96. http://dx.doi.org/10.1002/suco.201800265.

[49] Gao X, Ren X, Li J, Zhang Y. Bond behavior between steel reinforcing bars and concrete under dynamic loads. Struct Concr 2018;19(6):1806-17. http://dx.doi.org/10.1002/suco.201700205.

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

[51] Harajli MH. Comparison of bond strength of steel bars in normal-and high-strength concrete. J Mater Civ Eng 2004;16(4):365-74. http://dx.doi.org/10.1061/(ASCE)0899-1561(2004)16:4(365).

[52] fib Model Code for concrete structures 2010. http://dx.doi.org/10.35789/fib.BULL.0055.

[53] Harajli MH, Hout M, Jalkh W. Local bond stress-slip behavior of reinforcing bars embedded in plain and fiber concrete. ACI Mater J 1995;92(4):343-53, http://dx.doi.org/10.14359/999.

[54] Lemnitzer L, Schröder S, Lindorf A, Curbach M. Bond behaviour between reinforcing steel and concrete under multiaxial loading conditions in concrete containments. In: 20th international conference on structural mechanics in reactor technology. 2009.

[55] Gilbert RI, Nejadi S. 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.

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

[57] Ragueneau F. Fonctionnement dynamique des structures en b'eton : influence des comportements hystériques locaux. École normale supérieure de Cachan; 1999.

[58] Dur R, Farias MM, Pedroso DM, Meschke G. Reinforcing bars modelling using a rod-solid interface element without the need for mesh compatibility. Finite Elem Anal Des 2021;197:103634. http://dx.doi.org/10.1016/j.finel.2021.103634.