Experimental and numerical evaluations of the bond behaviour between ribbed steel rebars and concrete

The study of interfacial behaviour between ribbed steel rebars and concrete is a subject that has been widely studied. However, the definition of the bond stress distribution throughout the embedded length of the steel rebar is still controversial due to the difficulty of experimentally obtaining such distribution for a fixed load magnitude. It is also undeniable its relevancy for the better understanding and model reinforced concrete (RC) structures. So, the definition of the local behaviour between the ribbed steel rebar and concrete is critical to correctly simulate the adherence between both materials. In this matter, the local bond-slip models recommended in codes seem to satisfy some researchers while others suggest prudence in using them. Therefore, only choosing the correct bond-slip relationship may lead to exact interpretations and conclusions of the structural behaviour of a concrete structure but with the existing different bond-slip types, researchers can be misled inadvertently. This work aims to clarify some of these aspects by numerically simulating several pull-out tests under different conditions and checking their influence (or not) on real-scale specimens. After the validation of the numerical model through a proposed new bond-slip relationship, other parameters were studied also. Although the type of the bond-slip relationship influences the detachment of the steel rebar from the concrete, the yielding of the former material was found to be the main parameter that masks the differences in the behaviour of real-scale RC structures when different types of bond-slip relationships were considered in the numerical simulations.


Introduction
The success of reinforced concrete (RC) structures lies in the good bond quality between ribbed steel rebars and concrete. Therefore, it is not surprising that this topic has motivated several experimental (e.g. [1][2][3][4][5]) and numerical (e.g. [6][7][8]) studies where the pull-out test is often used [9] rather than the tests cited in ACI 408R-03 [10] such as the beam end, beam anchorage, and beam splice tests. Among other subjects, the experimentally based analyses cover the influence of high-performance concrete [1], lightweight aggregate concrete [4], steel corrosion [2,5], or bond conditions of the steel rebars [3] on the bond-slip relationship between the steel rebar and the concrete. However, these studies generally miss showing the ability of the derived bond-slip relationships to predict the load vs. displacement response of the steel rebar-to-concrete joint with precision. Thus, to mitigate that gap, some original numerical implementations of the bond-slip relationships with the Finite Element Method (FEM) can be found in the literature [11][12][13], which allows us to simulate, as close to the reality as possible, different type of RC structures.
Despite the existence of codes, such as Model Code 2010 [14], recommending the use of a bond-slip relationship for steel rebars embedded into concrete, some authors [15,16] are not comfortable with its use. Although based on several experimental data, two main criticisms are pointed out in the recommendation in Model Code 2010 [14]. One is the inconsistency provided by the model when the crack width in concrete is predicted. To explain this, Tan et al. [15] argued that the definition of the anchorage length may provide good results in small specimens but when it comes to real-scale RC structures, some contradictory results can be found. Consequently, the generalized applicability of the Model Code 2010 [14] model could be quite limited and prevent, therefore, its use in more generalized situations. Beconcini et al. [16] had difficulties obtaining precise results when the bond-slip relationship recommended in Model Code 2010 [14] was used to simulate the bending rotation response of RC beams to column joints. So, they decided to look closely at the bond-slip relationship and studied its influence on the numerical results by comparing it with the experimental data. One conclusion was then obtained, the use of a bond-slip relationship with lower slips associated with the bond stresses, whether in the elastic stage (pre-peak stage) or softening stage (post-peak stress), led to more precise numerical results when compared with the experiments than those obtained from the use of the bond-slip recommended in Mode Code 2010 [14]. To better consider the real conditions that the steel rebar is subjected to, Beconcini et al. [16] suggested that the pull-out test could be changed. The embedded length corresponding to five diameters of the ribbed steel rebar (ϕ s ) is another aspect suggested by these authors [16] that could be changed, i.e. the increase of the embedded length value was recommended.
Nowadays, the experimental determination of the bondslip relationship is still widely accepted and, for its determination, the recommendations made by RILEM [17] and Model Code 2010 [14] are commonly followed, e.g. [18][19][20][21]. In the case of RILEM [20], an embedded length of 4.2ϕ s is recommended to be used in a pull-out test whereas Model Code [14] recommends the use of an embedded of 5ϕ s . Behind these recommendations of using short embedded lengths, two main reasons can be pointed out. The first one is to avoid any other failure modes than the interfacial detachment of the ribbed steel rebar from concrete. The other one aims to facilitate the calculation of the bond stresses developed within the steel rebar-to-concrete interface by assuming, based on the loads transmitted to the steel rebar, an average bond stress distribution throughout the embedded length. The main problem here is when the loads are low and the embedded length is not fully mobilized and some regions close to the steel rebar unpulled end are still undeformed. Consequently, the elastic (E) stage of the bond-slip relationship, i.e. before the maximum bond stress is reached, will be mistakenly predicted. Despite this E stage representing a quite short part of the complete bond-slip relationship, its influence is reflected in the embedded length that can be mobilized when the steel rebar is being pulled out from the concrete. In Model Code 2010 [14], the E stage is defined by a power function where α is a dimensionless exponential parameter which, by recommendation, it is set equal to 0.4. Considering also that to obtain a constant bond-slip relationship it only needs to take α = 0, the corresponding load vs. displacement response can be defined by a square root equation [22] suggesting, therefore, that a misunderstanding between the definition of the local and global behaviours exists here.
Furthermore, the use of a short embedded length in a pullout test such as 4.2ϕ s or 5ϕ s , may lead to the determination of higher displacements [22], which agrees with the observations made by Beconcini [16] as mentioned before. So, the discussion on the most appropriate embedded length that should be used in a pull-out test is actual and some researchers are suggesting an embedded length of 10ϕ s instead [19]. However, according to Belarbi et al. [23], such embedded length may cause excessive bond stresses during the pullout test. So, to obtain a compromised solution between the standard 5ϕ s and 10ϕ s recommended by RILEM [17] for, respectively, a pull-out and beam test, they advised the use of an embedded length of 7.5ϕ s as a compromise solution between these two tests [23]. Nevertheless, none of these studies has proved yet the consistency between the local and the global behaviours independently of the embedded length used in the standard pull-out test. However, in a previous study carried out by the author [22] that consistency was pursued by adopting a simple bond-slip relationship with two constant plateaux. For analytical convenience, the E stage was not considered and the first plateau corresponded to the maximum bond stress developed within the steel rebar and the concrete. To calculate this maximum bond stress, an original method was proposed in [22], which consisted in approximating, through the lowest square minimization process, the experimental strains in the steel rebar and the strains calculated from a formula based on a constant bondslip relationship. In this minimization process, the strains should be representative of the same displacement measured at the steel rebar pulled end and only the specimens where the steel rebars have ruptured were considered, i.e. the specimens with the highest embedded lengths. The specimens with the lowest embedded lengths were useful to calculate the residual bond stresses developed in the second plateau of the bond-slip relationship, i.e. corresponding to the friction (F) stage. Again, for analytical convenience, a discontinuous transition between these two plateaux was considered, which is a relevant limitation of the model since the transition between both stages experimentally observed could not be predicted with high precision [22].
In the present work, an improved nonlinear bond-slip relationship is proposed in which, the discontinuity between the maximum and residual bond stresses is solved through a nonlinear formula that smoothed that transition. With no discontinuities, the proposed new bond-slip relationship avoids however, the analytical definition of the complete detachment of the ribbed steel rebars from the concrete since the resulting 2nd order differential equation that governs the full detachment process of the steel rebar from the concrete has no known analytical solution. That is why the finite element method (FEM) is used to predict a series of experiments that after being validated allowed us to carry out an extensive parametric study where different parameters were analyzed, such as: (i) the local bond-slip relationship and its influence on the surrounding concrete; (ii) embedded length; (iii) concrete strength; (iv) dimensions of the concrete block; (v) yielding stress of the ribbed steel rebars; and (vi) diameter of the ribbed steel rebar. This parametric study aims to identify which parameter (or parameters) may have critical interference with the detachment process between a steel rebar and a concrete block. In the end, the load vs. mid-span displacement curves obtained from three flexurally tested real-scale RC beams available in the literature [24][25][26][27][28] were also modelled and the influence of the bond-slip relationship type was analyzed. These tests were monotonically loaded until the yielding of the steel reinforcements. All results were quite close to the experimental data independently of the bond-slip relationship used with no meaningful differences between all simulations as well. The main reason is the very small slips developed between the steel rebars and the concrete before the yielding of the steel reinforcements. However, due to the different nature of the chosen bond-slip relationships (proposed new model, Model Code 2010 [14], and Bigaj's model [29]) considered in this work, different crack patterns were observed from all the simulations carried out in this study.

Mechanical properties of the materials
The mechanical properties of the materials experimentally defined in another work of the author [22] were used here as reference. Briefly, the mean strengths of the concrete are f cm = 37.2 MPa and f cm,cube = 42.3 MPa after 28 days of age whilst, at the time of the experiments, i.e. after 83 days, the mean strengths of the concrete are f cm = 40.8 MPa and f cm,cube = 46.3 MPa. Therefore, the latest values were used here. Furthermore, three different diameters (ϕ6, ϕ8 and ϕ12) of the ribbed steel rebars were also considered. Based on four tensile tests per diameter, the mean yielding stress of the 6-mm diameter steel rebar is f smy = 560.7 MPa whereas the 8 and 12-mm diameters steel rebars are f smy = 530.5 MPa and f smy = 523.8 MPa, respectively. At that point, the mean strain of the 6 and 8-mm diameter steel rebar is ε smy = 0.28% whilst the 12-mm diameter steel rebar is ε smy = 0.27%.

Geometry of the specimens
A total of 33 specimens were considered in the experimental programme where the diameter of the steel rebars and the embedded lengths varied. As mentioned before, and due to some limitations of the available testing equipment, the ribbed steel rebars varied between 6 to 12 mm, and the embedded lengths depended on the diameter of the steel rebars. These variations aim to produce specimens where the embedded lengths were too short and sufficiently long enough so it could be obtained different failure modes. For instance, with the too-short embedded lengths, it would be expected to obtain an interfacial failure mode between the steel rebar and the concrete block whilst the rupture of the steel rebar would be expected to occur when longer embedded lengths were considered. The value of the anchorage lengths was calculated according to Model Code 2010 [14] and whose results were considered as reference values: 140, 180 and 270 mm for the specimens with 6, 8 and 12-mm diameter steel rebars, respectively. So, the specimens with the shortest embedded lengths had less than 100 mm from the reference embedded length whereas the specimens with a longer embedded length had an extra length of 100 mm. Table 1 summarizes this classification given to the embedded lengths considered in this study.
The specimens consisted in embedding the ribbed steel rebar in the centre of a concrete block with a cross-sectional area of 200 × 200 mm 2 as recommended in EN 10038 [30]. Depending on the diameter of the steel rebar, the length of the concrete block varied. To ensure that the embedded length was preserved after the casting of the concrete blocks, a PVC tube was placed on the steel rebars at the bottom of the concrete block creating, therefore, an unbonded length there. Figure 1 shows the geometry and dimensions of all specimens briefly described here.

Brief description of the test procedure
The test setup consisted of pulling out the ribbed steel rebar from a concrete block. To that end, a steel frame (no. 14 in Fig. 2) was fixed to the strong concrete floor of the laboratory (no. 19 in Fig. 2) through two Dywidag steel rods (no. 1 in Fig. 2). Before initiating the tests, the specimens were placed first at the centre of the steel frame and on top of two C-shape steel profiles (no. 5 in Fig. 2) that were supported by two hollowed square steel profiles (no. 6 in Fig. 2) at its centre. The purpose of using these two hollowed steel profiles was to allow the installation of a displacement transducer (no. 8 in Fig. 2) that measured the displacements of the unloaded end of the ribbed steel rebars. A 300-mm of side steel plate with 20 mm thickness and with a hollow at its centre with 30 mm of diameter (no. 7 in Fig. 2) was placed at the top of the concrete block and it was used as a reaction steel plate in the test. Then, six threaded rebars with a diameter of 12 mm each (no. 10 in Fig. 2) were used to fix the concrete block to the steel double C-shape profile (no. 5 in Fig. 2). Since the displacements at the steel rebars were very important to measure and define the maximum bond stress developed within the interface between the steel rebar and the concrete block, the recommendation made in [17] of letting an initial unbonded length of 5ϕ s was not followed. Otherwise, it would have made it impossible to measure the displacements at the pulled end of the steel rebar. The influence of the concentration of normal stresses at the steel rebar-to-concrete interface was mitigated by adopting a small hollow in the reaction steel plate. However, the size of this hollow was sufficient enough to introduce a steel nut (no. 15 in Fig. 2) with three fixing points (through three small steel bolts). To this steel nut, three small steel cantilevers were welded which allowed us to measure the displacements of the pulled end of the steel rebar. Three displacement transducers were placed on the reaction steel plate and each one measured the displacements of each three steel small steel cantilevers. The average of these three displacements of the steel rebar pulled end was considered as the displacement of the steel rebar pulled end.
The loads applied to the specimens were ensured by a hydraulic jack (no. 3 in Fig. 2) that was placed at the top of the steel frame. A 12-mm diameter threaded rebar (no. 10 in Fig. 2) passed through the hydraulic jack and a grip (no. 9 in Fig. 2) was attached to its bottom end. The top of this threaded rebar was fixed with a steel nut. Between this fixing point and the hydraulic jack (no. 3 in Fig. 2), a pressure cell of 200 kN (no. 2 in Fig. 2) was placed which allowed the measurements of the loads transmitted to the ribbed steel rebar during the pull-out test. Finally, by switching on the hydraulic pump (no. 16 in Fig. 2), the test of the specimen began. The displacements and loads measurements were saved with the help of a data logger (no. 18 in Fig. 2) that sends the signal to a desktop computer (no. 17 in Fig. 2) allowing us to save the data after the pull-out test ended for later analysis.

ID of the specimens and summary of the experimental results
To easily identify the specimens and their geometrical characteristics, the following nomenclature is adopted: "a"-L "b"-f "c"-Y "d"-D "e"-ϕ "f"-"g", where "a" stands for the type of specimen which, in the case of the experimental specimens, "Sp" is considered; "b" correspond to the embedded length; "c" stands for the mean compressive strength of concrete obtained from uniaxial tests of concrete cylinders (f cm ); "d" corresponds to the yielding stress of the ribbed steel rebar; "e" denotes for the side of the concrete block; "f" is the diameter of the ribbed steel rebar; and "g" corresponds to the number of repetitious tests carried out under the same conditions. For instance, the third specimen tested with an embedded length of 80 mm and with an 8-mm diameter steel rebar is specimen Sp-L80-f37.2-Y530.5-D200-ϕ8-03 which the concrete strength is 37.2 MPa, the yielding stress of the steel rebar is 530.5 MPa and the cross-sectional area of the concrete block is 200 × 200 mm 2 . Table 2 shows the ID of the specimens and summarizes also the main results obtained from each test. The results show that the maximum loads tend to increase with the embedded length. However, for the specimens with a 12-mm diameter steel rebar such a trend did not occur because the rupture of the steel rebar was always reached in all these specimens and therefore, the maximum loads shown in Table 2 correspond, approximately, to the rupture load of those steel rebars.
In terms of the maximum bond stress developed within the steel rebar-to-concrete interface, it can be noticed that only using long embedded lengths it is possible to determine this parameter. However, to define the F stage of the bondslip relationships, specimens with short bonded lengths should be analyzed, i.e. specimens whose failure modes have an interfacial source between the steel rebar and the concrete. In these cases, as can be seen from Table 2, the maximum bond stress is quite low. To bypass this drawback, the average bond stress is still commonly calculated instead, e.g. [31][32][33][34][35].

The finite element model
To make the parametric study and highlight the influence that each one has on the bond response of the ribbed steel rebars from concrete, several numerical models were simulated according to the Finite Element Method (FEM). The latest 5th version (v5.9.1) of the commercial software ATENA 2D was used [36], which is known to be quite suitable for modelling RC structures. Aiming to significantly reduce the number of unknowns and non-linear equations to be solved through the Newton-Raphson method, the 2D version of ATENA software was used rather than its 3D version. At the same time, the computational effort to solve the numerical simulations was reduced without losing the precision of the final results. All numerical models of the pull-out tests were discretized through a mesh with 5 mm quadrilateral finite elements as shown in Fig. 3.
The loads applied to the specimens were simulated with a monotonic displacement control rate of 0.005 mm per step. A total of 3000 steps were simulated. Therefore, unless the failure of the model occurred at an earlier displacement, a maximum displacement of 15 mm at the top of the steel rebars was assumed. At each step, the corresponding loads transmitted to the steel rebars were controlled through two monitoring points (no. 6 in Fig. 3) located at two pinned supports spaced at the same distance from the centre of the models. Although the displacements were known at each loading step, two monitoring points (no. 5 in Fig. 3) located at the bottom corners of the steel plate (no. 3 in Fig. 3) were used to check the numerical simulations process and identify any possible errors during the numerical simulations.
The ribbed steel reinforcements were modelled through a truss (or tie) element which was defined through a multi-linear stress-strain relationship. Therefore, the real constitutive behaviour of the ribbed steel rebars was approximated with several finite lines. Similarly, and despite the nonlinearities of the local bond-slip relationships herein considered, multilinear lines were also considered to approximately define them. Nevertheless, two other and well-known local bondslip relationships are available in the ATENA 2D database. The first one is recommended in Model Code 2010 [14] whilst the second one is the bond-slip relationship proposed by Bigaj [29]. It should be mentioned that in these two bondslip relationships a good bond quality between the ribbed steel rebars and the concrete element was always assumed.
At the top of the truss element, the relative displacements between the ribbed steel rebar and the concrete block were prevented and no interfacial slips can be developed at that point. Thereby, no premature detachments of the truss    Fig. 3) from the loaded steel plate (no. 3 in Fig. 3) can occur. On the opposite side, i.e. at the ribbed steel rebar unloaded end, no interfacial slip constraints were considered. Therefore, and throughout the embedded length, the development of the bond stresses is ruled by the adopted local bond-slip relationship, and the detachment of the ribbed steel rebar (no. 4 in Fig. 3) from the concrete block (no. 1 in Fig. 3) is always guaranteed. Nevertheless, it should be bear in mind that such options will not compromise the full failure mode of the specimens whatever it could be, e.g. full interfacial detachment of the ribbed steel rebar from the concrete block, rupture of the steel rebar or splitting of the concrete block.
To conveniently simulate the mechanical properties of the materials, the following models were used in ATENA 2D [36]: (i) the "Constitutive Model SBETA" material to model plain concrete; (ii) the "Reinforcement" material to model the ribbed steel rebars; and (iii) the "Reinforcement Bond" material to model the local bond-slip relationship between the ribbed steel rebar and the concrete element. The behaviour of concrete includes the following effects in all simulated models [36]: (i) nonlinear behaviour in compression including hardening and softening; (ii) fracture of concrete in tension based on the nonlinear fracture mechanics; (iii) biaxial strength failure criterion; (iv) reduction of compressive strength after cracking; (v) tension stiffening effect; (vi) reduction of the shear stiffness after cracking (variable shear retention); and (vii) rotated crack model. It should be mentioned that the model used to simulate the concrete is complex and it is beyond the scope of this work to deeply describe it. Therefore, the readers are strongly advised to carefully look at ATENA's theory manual [36] for better and more detailed information about this model. Nevertheless, this material model has been successfully used to numerically simulate RC structures, e.g. [37][38][39][40], with particularly excellent results when the rotated crack approach is used rather than the fixed crack approach [25,26]. For this reason the rotated crack approach model is exclusively used in this study.

Proposed new bond-slip relationship
The simplified local bond-slip relationship proposed by the author in [22], and considered in this study as a reference, has two constant stages. The first one has maximum bond stress of τ bmax = 11.75 MPa if the 6-mm diameter steel rebar is considered whereas τ bmax = 12.64 MPa if 8 and 12-mm diameter steel rebars are assumed. The second stage corresponds to the residual stress where τ bf = 1.95 MPa for the 6-mm diameter steel rebar and τ bf = 3.99 MPa in the other two cases with diameters of 8 mm and 12 mm. The transition between these two stages was made at an interfacial slip (s 0 ) of 1.601 mm in the lowest diameter and 2.446 mm in the two largest diameters. However, in the present proposed new bond-slip relationship, s 0 will be recalculated since the transition between τ bmax and τ bf is now smoothed and its value will change.
As mentioned, the bond-slip relationship proposed in [22] was a simplified relationship due to analytical convenience that allowed the derivation of closed-form solutions to calculate slips, strains, and bond stresses distributions throughout the embedded length as well as the estimation of the load-displacement responses of the pull-out tests of ribbed steel embedded in a concrete block. That simplified bond-slip relationship had two limiting simplifications. One is the assumption of no elastic (E) stage due to its very short domain when compared with the other stages. The other one is the abrupt and discontinuous transition between the constant (C) stage and the friction (F) stage which, in reality, is much smoother and continuous. Despite the proposed new bond-slip relationship also neglects the E stage, the smooth and continuous transition between the C stage and the F stage is now considered. The proposed new bond-slip relationship (τ-s) is mathematically defined according to: where τ bmax and τ bf are the maximum bond stress and the residual bond stress, respectively; a is a parameter to be calibrated with the experimental results; and s 0 is the mid slip of the transition between the maximum and residual bond stresses (see Fig. 4).
The proposed new bond-slip relationship defined in (1) is adjusted to the experimental data though the lowest square minimization process between the predicted and experimental bond stresses is used, in which parameters a and s 0 are determined according to: where τ b,pred (s) and τ b,exp (s) are the predicted bond stresses obtained from Eq. (1) and the experimental data, respectively. Figure 5 compares the predicted bond-slip curves with the experimental results after carrying out the minimization process defined in (2). The results conducted to a = 0.8838 mm −1 and s 0 = 3.436 mm for the specimens with a 6 mm steel rebar diameter whereas, for higher diameters, i.e. 8 and 12 mm, the results led to a = 0.7074 mm −1 and s 0 = 5.727 mm. These results suggest that the increase of the bar diameter makes the transition between τ bmax and τ bf smoother since parameter a decreased with the increase of the diameter of the steel rebar.
The second-order differential equation that governs the detachment of the steel rebar from the concrete is defined according to [22,41]: where where E s and ϕ s are the elastic modulus and the diameter of the steel rebar, respectively; E c and A c are the elastic and cross-sectional area of the concrete block, respectively. When Eq. (1) is introduced into Eq. (3), no analytical solution can be derived and only through the implementation of a numerical strategy will be possible to obtain a solution. Alternatively, the bond-slip relationship can be implemented in a commercial finite element code which was the solution followed in this work.

Parametric analyses
In this section, the influences of different parameters on the detachment of ribbed steel rebars from concrete are presented and discussed, such as the influence of the: (i) test set-up; (ii) bond-slip relationship; (iii) embedded length; (iv) concrete strength; (v) dimensions of the concrete block; (vi) yielding of the ribbed steel rebar; (vii) diameter of the ribbed steel rebar; and (viii) bond-slip relationship on the surrounding concrete. It should be pointed out also that the comparisons with the experimental data are made, unless mentioned otherwise, with the 8-mm diameter steel rebar specimens with an embedded length of 80 mm.

Influence of the bond-slip relationship
Four local bond-slip relationships were considered. The two first ones are the simplified bond-slip relationship proposed by the author in [22] and the improved and new bond-slip relationship defined in (1). The third bond-slip relationship is the one recommended in Model Code 2010 [14] which is defined as follows: where s 1 and s 2 are the slips at the beginning and end of the plateau under maximum bond stress; s 3 defines the initiation of the F stage; and α is a dimensionless parameter that despite it could take a value between 0 and 1, is recommended to use α = 0.4. For good bond conditions between the steel rebar and concrete, Mode Code 2010 [14] recommends using also s 1 = 1.0 mm, and s 2 = 2.0 mm, while for the quantification of s 3 , it should be used the clear distance between ribs which, in the present case, is 5.0 mm. The fourth bond-slip relationship herein considered is the one proposed by Bigaj [29]. This bond-slip relationship does not consider the F stage and it has a bi-linear softening stage. It assumes also, after a final slip, no interaction between materials, i.e. τ = 0 MPa after the final slip is exceeded which means that, after this point, the bond stresses are no longer transferred between both materials. However, the maximum bond stresses are close to those experimentally determined. Bigaj's model [29] also considers the influence of the ribbed steel rebar diameter by affecting only the slips and making no changes to the values of the bond stresses.
The numerical results of the load-displacement curves (P-δ 0 ) are confronted against the homologous curves experimentally obtained (see Fig. 6). The load-displacement responses of the shortest embedded lengths of L = 40 mm (with ϕ s = 6 mm) and L = 80 mm (with ϕ s = 8 mm) show that the use of the proposed new bond-slip relationship leads to very close results with the experiments. However, for higher embedded lengths that lead to the rupture of the ribbed steel Fig. 6 Influence of the local bond-slip relationship on the load-displacement response of the specimens with: a a 6-mm diameter steel rebar; b an 8-mm diameter steel rebar; and c a 12-mm diameter steel rebar rebar, no meaningful differences between the simplified bond-slip relationship and the new one can be observed. Another interesting aspect found from the results shown in Fig. 6 is the P-δ 0 curves obtained from the use of Bigaj's model [29]. This model is the unique one that predicts the detachment of the ribbed steel rebar from the concrete block for the two shortest embedded lengths in each diameter of the steel rebar. This reveals that Bigaj's model [29] leads to the highest anchorage lengths. Moreover, this model [29] was the one that showed the highest differences from the experimental results, particularly from those specimens with the shortest embedded lengths.
Comparing with the experimental data obtained from specimens Sp-L40-f37.2-Y560.7-D200-ϕ6-01|3, Sp-L90-f37.2-Y560.7-D200-ϕ6-01|2 and Sp-L80-f37.2-Y530.5-D200-ϕ8-01|3, the numerical results obtained from the use of the bond-slip relationship recommended in Model Code 2010 [14] were not very precise as well. However, its precision seems to increase in the specimens where the rupture of the ribbed steel rebars was experimentally observed, i.e. in the specimens with the longest embedded lengths. When the embedded length increases, the precision of all bond-slip relationships is almost the same, which means that the shape of the bond-slip relationship does not influence the P-δ 0 responses in those cases. Still, it should be noticed that the bond-slip relationships recommended in Model Code 2010 [14] and proposed by Bigaj [29] predict the rupture of the steel rebars at a higher displacement than that predicted by the proposed new bond-slip relationship.

Influence of the embedded length
The influence of the embedded length on the P-δ 0 curve was analysed by changing the embedded length in the numerical models. To obtain an interfacial failure and thereby obtain the corresponding P-δ 0 curve, the yielding of the steel rebars was ignored in all simulations. The following embedded lengths were considered: 30 mm (which corresponds to 5ϕ s of ϕ s = 6 mm), 40 mm (5ϕ s of ϕ s = 8 mm), 60 mm (5ϕ s of ϕ s = 12 mm), 80, 150, 180 and 280 mm. The results are shown in Fig. 7. Although the graphs were not limited by the yielding load of the steel rebars, it should be noticed that beyond this load the results were put behind a shadow that becomes darker after the corresponding rupture load of the steel rebar is exceeded.
Based on the results shown in Fig. 7, it is possible to observe that the increase in the embedded length increases the strength of the specimens, i.e. corresponding to the detachment of the ribbed steel rebars from the concrete block. Ignoring the yielding of the steel rebar, before its detachment from the concrete block, residual bond stresses will develop within the interface between the steel rebar and concrete which, by assuming a constant distribution, explains this increase in the load capacity of the specimens. However, as the steel rebar is pulled out from the concrete block, the bonded area is reduced which, from the point of view of equilibrium, will make the load decrease with the increase of the displacements even after the model has reached its maximum load capacity. This leads to a plateau which is due to the existence of friction between the ribbed steel rebar and the surrounding concrete and since the bond stresses at this F stage are constant, the value of the corresponding load can be calculated according to [22]: where L is the embedded length of the ribbed steel rebar. Figure 8 shows the maximum load vs. embedded length obtained from the use of the proposed new bond-slip relationships, recommended in [14] and proposed by Bigaj [29]. This graph also reveals that the increase of the embedded length increases the maximum load transmitted to the steel rebar. However, it should be noted that with Bigaj's model [29], the maximum load of the specimens with high Fig. 7 Influence of the embedded length on the load-displacement curve when ϕ s = 8 mm is assumed embedded lengths tend to a plateau, which means that, if the steel rebar did not yield, it would be possible to define an embedded length beyond which the loads would not continue to increase further. Unlike the other two models where residual stresses are considered (F stage), Bigaj's model [29] assumes that after a final slip, no further bond stresses are transmitted between the ribbed steel rebar and concrete, which explains this plateau at maximum load for higher embedded lengths.

Influence of the concrete strength
To analyse the influence of the concrete strength on the P-δ 0 response, the numerical models of the specimens with a steel rebar with 8 mm of diameter and an embedded length of 80 mm were used. The concrete strength in these models was changed and the following concrete strengths were considered: C20/25, C25/30, C30/37 and C35/45. Unlike Model Code 2010 [14] and once the experimental programme did not consider different concrete strengths, the proposed new bond-slip relationship has no parameters that may consider such influence. Therefore, to make as many valid comparisons as possible between bond-slip relationships, the values of the maximum bond stress (τ bmax ) and residual bond stress (τ bf ) were modified according to the same proportion stated in Model Code 2010 [14]. To avoid that the yielding of the steel rebar could be limiting the evaluation of the final results, the yielding of the steel rebars was ignored once again. By assuming this, the failure modes of the models were conditioned by ensuring always that the interfacial detachment of the steel rebar from the concrete block would occur.
A good bond condition between the ribbed steel rebar and the concrete was assumed which, according to the Model Code 2010 [14] recommendation, the maximum bond stress should be determined according to: where f cm is the mean strength of concrete. Thus, for the additional different concrete strengths herein considered, the maximum bond stresses adopted in the proposed new bond-slip relationship are: 9.98 MPa (ϕ s = 6 mm) and In the case of Bigaj's model [29], for a ribbed steel rebar with a diameter of 8 mm, the C20/25, C25/30, C30/37 and C35/45 concrete lead to a maximum bond stress of 10.28, 11.03, 12.00, and 13.02 MPa, respectively. Figure 9 shows that all the bond-slip relationships led to higher load capacities with the increase of the concrete strength. This can be explained by the increase of the maximum bond stresses in all bond-slip relationships with the increase of the concrete strength. The results also show that such an increase in the load capacity of the models is followed by a slight increase of the initial stiffnesses in the P-δ 0 curves of the simulations that used the bond-slip relationships in [14] and [29].
Moreover, the use of Bigaj's model [29] is the unique bond-slip relationship that leads to the interfacial detachment of the ribbed steel rebar from the concrete independently of its strength (see Fig. 9c). Except the model with a C35/45 concrete where the yielding of the steel rebar is exceeded, the use of the proposed new bond-slip relationship led to the same failure modes in all the other cases, i.e. interfacial failure mode (see Fig. 9a). However, by not reaching the rupture load of the steel rebar in this case, the interfacial detachment of the steel rebar from concrete will be ensured. Bigaj's model [29] Proposed new model Finally, the use of the bond-slip relationship recommended in Model Code 2010 [14] is the unique one that led to the rupture of the steel rebar when C30/37 and C35/45 concretes were used. Therefore, to prevent the rupture of the 8-mm diameter steel rebar, the embedded length of 80 mm should be decreased.

Influence of the dimensions of the concrete block
The dimensions of the concrete block to adopt in a pullout test are another important aspect that deserves to be analysed. For instance, if the cross-sectional area is quite limited then, it would be expected that the maximum load could be affected. Consequently, the precise definition of the bond-slip relationship will be compromised. The failure mode of the specimens may be also shifted, e.g. from interfacial between the steel rebar and the concrete (or from rupture of the steel rebar) to splitting of the concrete block. For this motive, the following cross-sectional areas for the concrete block were considered: 50 × 50 mm 2 , 100 × 100 mm 2 , 150 × 150 mm 2 and 200 × 200 mm 2 . To understand how these dimensions may be influenced by the embedded length and by the diameter of the steel rebars, two embedded lengths (L = 80 mm and 280 mm) and two ribbed steel rebars (ϕ s = 8 mm and ϕ s = 20 mm) were considered. It should be kept in mind that the concrete block cross-sectional area of 200 × 200 mm 2 was the dimension adopted in the experimental programme which corresponds to the minimum dimension recommended in EN 10080 [30]. The three previously described bond-slip relationships were considered as well and the results are shown in Fig. 10.
All the results with an embedded length of 80 mm seemed to agree on one point: the cross-sectional area of 50 × 50 mm 2 is quite small to carry out a pull-out test since a splitting failure mode in the concrete block is predicted. As result, the load capacity of the specimens is reduced significantly when compared with the other simulations with the highest cross-sectional areas of concrete. Nevertheless, Bigaj's model [29] led to the lowest load reduction of the load capacity of the remaining simulations, approximately 15.2% of the average maximum load capacities predicted in the other simulations. On the other hand, the highest reduction occurred in the simulations where the bond-slip relationship recommended in Model Code 2010 [14] was used, i.e. a reduction of approximately 45.3% of the average maximum load capacities predicted in the other simulations. So, with a steel rebar with a diameter of 8 mm and an embedded length of 80 mm, a cross-sectional area of 100 × 100 mm 2 should be sufficient enough to obtain an interfacial failure mode (see top Fig. 10a-c).
However, to obtain an interfacial failure mode with higher steel rebar diameters, it is expected that the cross-sectional area of the concrete block needs to be increased. Thus, when the steel rebar diameter increased to 20 mm, the use of the proposed new bond-slip relationship as well as the use of the one recommended in Model Code 2010 [14] showed that only with a cross-sectional area of 200 × 200 mm 2 it is possible to obtain an interfacial failure mode between the steel rebar and the concrete block. However, according to the results that used Bigaj's model [29] a cross-sectional area of 150 × 150 mm 2 seems to be sufficient enough to obtain the same interfacial failure mode of the specimen (see bottom of Fig. 10c). Furthermore, by increasing the cross-sectional area of the concrete block from 100 × 100 mm 2 to 150 × 150 mm 2 , a short load capacity increment of approximately only 4.2% was predicted, i.e. from 162.7 kN to 169.9 kN.

Influence of the yielding stress of the ribbed steel rebar
To analyse the influence of the yielding stress of the ribbed steel rebars in the P-δ 0 curve obtained from a pull-out test, new four numerical models, based on the previous ones, were developed where the following yielding stresses were assumed: 500, 750, 1000 MPa and assuming no yielding at all. It is important to bear in mind that the perfect elastic-plastic constitutive relationship of the ribbed steel rebars was considered in these simulations. The results are shown in Fig. 11. It can be seen in this figure that the lowest initial stiffnesses were obtained from the bond-slip relationship recommended in Model Code 2010 [14]. When the yielding stress of 500 MPa is adopted, the simulations carried out with the proposed new model and with the one recommended in Model Code 2010 [14] predicted both the yielding of the steel rebar. Since a perfect elastic-plastic model was used to simulate the stress-strain relationship of the steel rebar, Fig. 11a, b show a plateau at a maximum load which corresponds to the yielding load of the steel rebar with f ym = 500 MPa. Moreover, Fig. 11b shows another P-δ 0 curve that is different from the others, which is the curve when f ym = 530.5 MPa is assumed, i.e. corresponding to the experimental test. Here, the steel rebar after yielding, the load continued to increase, although at a lower rate, due to the hardening of the steel rebar. Before the steel rebar could reach its rupture load, the interfacial detachment of the steel rebar from the concrete block is predicted which delayed the specimen reaching its final state, i.e. the F stage. As could be expected, the models where the yielding stress of 500 MPa was used, led to similar results to those experimentally obtained since the yielding stress is approximately the same (530.5 MPa). In the case of the steel rebars with yielding stresses higher than 530.5 MPa, the proposed new bond-slip relationship and the one recommended in Model Code 2010 [14] led to results, despite with different initial stiffnesses, that ended in a plateau with the same load in each case. The differences between the maximum loads reached in each case are due, of course, to the maximum bond stress (τ bmax ) adopted in each bond-slip relationship.

Influence of the diameter of the ribbed steel rebar
Like in previous subsections, to prevent any disturbance of the final results and lead to the interfacial failure mode of the models between the steel rebar and the concrete block, the yielding of the steel was ignored in the following analyses. Also, the original numerical models were changed by varying the diameter of the steel rebar as follows: ϕ6 mm, ϕ8 mm, ϕ10 mm, ϕ12 mm, ϕ16 mm and ϕ20 mm. Based on the results shown in Figs. 12a, b, it is possible to asseverate that a dependency between the maximum load and the diameter exists. Bearing in mind the equilibrium conditions of the steel rebar, with a higher diameter of the steel rebar the bond stresses will be distributed by a larger bonded area, which explains the load capacity increase of the specimen. After the detachment of the steel rebar, a plateau at residual load can be seen. The bond-slip relationship recommended in Model Code 2010 [14] led to higher maximum loads than that obtained from the proposed bond-slip relationship. This can be explained by the difference in the maximum bond stress used in those models where, for instance, in the model recommended in Model Code 2010 [14] τ bmax = 15 MPa whereas in the proposed model we have τ bmax = 12.46 MPa.
On the other hand, since the definition of the slips in Bigaj's model [29] increases with the diameter of the steel rebar, the increase in the loads observed in Fig. 12c can be justified by the embedded length adopted in these simulations. In other words, an embedded length of 80 mm is sufficient enough to produce a bond stress distribution throughout the embedded length that makes the load transmitted to the steel rebar increase. The increase of the slips with the diameter of the steel rebar whether at a maximum load magnitude or at the second post-peak load decay are both quite clear in Fig. 12c.

Influence of parameter α
Parameter α is a dimensionless parameter that is used to define the bond-slip relationship proposed in Model Code 2010 [14] and, despite recommended to use α = 0.4, it is allowed to use a value that varies between 0 and 1. Also, no reasons are presented in [14] that can guide engineers to choose another value within that interval. However, if it is assumed that α = 0, then the local bond-slip relationship has an initial and constant maximum bond stress which is similar to the simplified bond-slip model proposed by the author in [25]. On the opposite, i.e. when α = 1 is considered, the local bond-slip model has an initial linear shape. The Model Code 2010 [14] recommends α = 0.4 whether the bond conditions between the steel rebars and concrete are good or poor or if the failure mode is interfacial between the steel rebar and the concrete block or if the failure mode is due to the splitting of concrete. As could be expected, the definition of α will affect the P-δ 0 response of the specimens. Therefore, to evaluate Fig. 12 Influence of the diameter of the steel rebars on the load-displacement curve for an embedded length of 80 mm how this dimensionless parameter α influences the P-δ 0 curve of short and long embedded lengths, new numerical models were created in which the embedded lengths of 40 and 280 mm were considered and where parameter α was set equal to 0, 0.4 and 1. Like in the previous analysis, to highlight the full differences between local bond-slip models, the yielding of the steel reinforcement was ignored.
Based on the P-δ 0 results shown in Fig. 13 it can be asseverated that when α = 0 is considered, the results obtained from the Model Code 2010 [14] are, until the corresponding value of the yielding of the steel rebar, almost identical from the ones obtained from the use of the proposed new bondslip relationship. In this case, these results are also quite close to those obtained from the experiments. The results that were far away from the experimental ones correspond to the numerical models with α = 1. Therefore, it can be stated that the bond-slip relationship recommended in Model Code 2010 [14] provides close results to the proposed new bondslip relationship when an initial rigid-plastic bond-slip relationship is assumed, i.e. when α = 0, or, at least, when α is much smaller than the recommended value of 0.4.

Influence of the bond-slip relationship type in the surrounding concrete
To analyse the influence that each bond-slip relationship has on the concrete, Figs. 14 and 15 show the principal maximum strains developed in the concrete block for 8-mm diameter specimens with embedded lengths of 80 and 280 mm, respectively. In Fig. 14, three different instants of the pullout test are shown and each one corresponds to the maximum load transmitted to the ribbed steel rebar. Since in each case this occurs at different displacements of the simulations, the other two results allow us to make comparisons with another case that shows the principal maximum strains in the concrete at maximum load. The results show that independently of the bond-slip relationship used, the full embedded length is mobilized, i.e. with no undeformed regions. Therefore, the interfacial detachment between the steel rebar and the concrete is always guaranteed (see top Fig. 14b). However, the use of the proposed new bond-slip relationship mobilized the full embedded length at a smaller displacement (δ 0 = 0.140 mm). On the opposite, the use of the bond-slip relationship recommended in Model Code 2010 [14] took the model to the highest displacement (δ 0 = 0.865 mm) at maximum load. In all cases, the principal main stresses at maximum load seemed to be concentrated at the steel rebar unpulled end. It should be noted also that the principal maximum strains obtained by the proposed new model in the three stages of the pull-out test are quite alike due to the smooth decay of the bond stresses with the slip increase. Moreover, the results obtained from Bigaj's model [29] (see Fig. 14c) allow us to observe that as the loads increase until their maximum value, the maximum principal strains increase and then, i.e. after exceeding this load point, tend to decrease. The principal maximum strains reach in each model were similar being equal to 0.0037, 0.0038 and 0.0031% in the proposed new model, Model Code 2010 [14], and Bigaj's model [29], respectively.
The simulations carried out of the specimens with an embedded length of 280 mm show that the use of the proposed new model concentrates the bond stresses in a smaller embedded length (see Fig. 15a). This may be justified by the fact that no E stage is considered in this bond-slip relationship. For this same reason, the maximum load is reached first in this case for a displacement of δ 0 = 0.190 mm. Furthermore, at this instant, the yielding stress of the steel rebar was reached which means that the maximum load is 26.7 kN, i.e. corresponding to the yielding load. Thus, this predicted failure mode is consistent with the one observed in the experiments, i.e. yielding of the steel rebar followed by its rupture. Figure 15b shows the results obtained from the use of the bond-slip relationship recommended in Model Code 2010 [14]. Here, it is possible to observe that the highest bond stresses developed within the steel rebar and the concrete interface are located at the steel rebar pulled end. The same is also predicted in the simulation that used  Bigaj's model [29] as can be seen from Fig. 15c. In these three simulations, the yielding of the steel rebar occurs at a higher displacement of δ 0 = 0.315 mm when the bond-slip relationship in Model Code 2010 [14] is used. In terms of principal strains predicted in these three simulations, the proposed new bond-slip relationship, Model Code 2010 [14] and Bigaj's model [29], led to the following maximum values of 0.0041, 0.0030, and 0.0033%, respectively.
Two other important aspects can be derived from Fig. 15 also. First, is the corresponding stress bulb predicted in the numerical simulations, which facilitates the understanding of using a concrete block with a sufficiently large crosssectional area to avoid the splitting of the concrete as previously discussed. The second one is the influence that such stress bulbs may introduce in concrete resulting, e.g. in an RC beam, column or slab, different crack patterns. Due to its importance, this second topic is discussed in the subsequent section.

Examples with real-scale RC beams
This section aims to emphasise how the bond-slip relationship type may (or may not) affect the load-displacement response of real-scale RC beams. To that end, three different simple supported RC beams tested by the author's group research and whose experimental results are available in the literature [24][25][26][27][28] are considered. In all cases, the RC beams have a cross-sectional area with a T shape and were all flexurally tested under a 4-point bending condition according to the system shown in Fig. 16. All the tests were conducted under a monotonic displacement control until failure, i.e. until the yielding of the tensioned steel rebars. Despite a brief description of these tests being made next, the readers are strongly advised to look at those published works [24][25][26][27][28] for further experimental details about them.
The RC beams had 3.30 mm long which allowed us to have a beam with 3.00 m of span. The loads were applied at one-third and two-thirds of the beam's span. The crosssectional areas of the RC beams were almost similar in all these studies [24][25][26][27][28] and the main differences lie in the mechanical properties of the concrete as well as the ribbed steel rebars used in each beam. All were designed to prevent a premature shear failure and therefore, steel stirrups Ø6//15 mm were assumed. To generate an RC beam under a realistic situation, the flexural reinforcement under tension, i.e. at the bottom of the beam, was set equal to 3Ø12. Figure 17 illustrates the transversal sections of all beams.
The mechanical properties of the concrete used in these RC beams were obtained from the uniaxial compression of 3 cubes with 150 mm of side. All tests were conducted according to Portuguese national standard EN 12390-3 [42]. The mechanical properties of the steel reinforcements were also defined through the uniaxial tensile testing of, at least, three samples and where standard EN 10002-1 [43] was followed. The summary of these mechanical properties can be found in Table 3. In the following subsections, the previous three bondslip relationships are considered also. The assumption of a rigid connection between the steel reinforcements and the concrete is also implemented whose results can be used as reference values. The numerical results obtained from these distinct situations are all confronted against the experimental ones. To identify which bond-slip relationship can provide better results, the precision of the numerical simulations is evaluated through the calculation of the Integral Absolute Error (IAE). This parameter is known to be sensitive to the deviation of a theoretical result and is commonly used for model assessment [1,44,45]: where P i num and P i exp correspond to the loads numerically predicted and those obtained experimentally at the mid-span displacement i of the test, respectively; and n corresponds to the number of steps considered in the numerical simulations until failure.
For the same motives presented before for the simulations of the pull-out tests, the 2D version of the commercial software ATENA v.5.9.1 [36] was used once again. Furthermore, to reduce the computation effort and time consumption of the simulations even more, only one-half of the RC beams were modelled due to symmetry conditions of the 4-point bending tests. All the models were discretized through a 20 mm quadrilateral mesh and, after applying the dead load of the beam and the corresponding weight of the test apparatus over the beam (approximately a total of 7.0 kN), the subsequent loads were simulated through a monotonic displacement control rate of 0.2 mm per step. The loads transmitted to the specimens were controlled through two monitoring points (no. 3 and 4 in Fig. 18) located at the roller support and the loading point whilst the mid-span displacements were monitored by a unique monitoring point located at the top of the RC beam (no. 5 in Fig. 18) and on the symmetry axis of the numerical model (no. 9 in Fig. 18). Like the pull-out tests, the steel reinforcements were simulated with truss elements and "Constitutive Model SBETA" was used to model concrete. The rotated crack approach in concrete was also considered. Fig. 17 Cross-sectional areas of the RC beams flexurally tested by: a Biscaia et al. [25,26]; b Carvalho et al. [24]; and c Franco et al. [27,28]. (Dimensions in mm)

RC beam flexurally tested by Biscaia et al. [25, 26]
The load vs. mid-span displacement obtained from the numerical simulations as well as the IAE values calculated for each model are shown in Fig. 19a until a maximum mid-span displacement of 30 mm. The results show no meaningful differences between using any of the bond-slip relationships herein considered. In an overall overview of the results, the precisions with the experimental data of all numerical simulations are very good being the maximum IAE value at a mid-span displacement of 30 mm equal to 2.63% when Bigaj's model is used and a minimum value of 2.35% when a rigid contact between the steel reinforcements and the concrete is assumed. Although the local bond adherence type does not influence the load vs. mid-span displacement results, the four principal maximum strain patterns of the RC beam are shown in Figs. 19b to 19e at the yielding of the steel reinforcements reveal some differences. The principal maximum strain patterns allow the identification of the crack pattern of the RC beam. For instance, using Bigaj's model [29] would lead to a crack pattern with more generalized cracks along the span of the RC beam, especially at the centre region where the bending moment is constant. In this region, five main cracks can be observed (see Fig. 19e). In this case, the predicted principal maximum strains reached the highest values (approximately 0.76%). Until the yielding instant, the bond stresses have not reached their maximum value and in this case, shown in Fig. 19e, the maximum bond stress developed in the RC beam is 3.66 MPa only. At this instant, the predicted maximum crack opening displacement (COD) is also 0.35 mm.
On the opposite, by assuming a rigid contact between the steel reinforcements and concrete, the numerical results show that the principal strains, although not much different, have the lowest values. At the yielding instant, the highest maximum principal strain (approximately 0.70%) is located between the support and the applied load. Fig. 18 Example of the 2D numerical model used to simulate the RC beam tested by Biscaia et al. [25,26] P/2  [14], and the Bigaj's model [29] are used, respectively A similar principal maximum strain pattern can be found between the numerical simulations with the proposed new bond-slip relationship and with the one recommended in Model Code 2010 [14]. The main difference between these two cases is the bond stress developed at this yielding instant. For instance, with the use of the proposed new bond-slip relationship, a value of 10.37 MPa was predicted whilst by using the bond-slip relationship recommended in Model Code 2010 [14] a maximum value of only 3.70 MPa was obtained. So, considering that the slips between the steel reinforcements and the concrete are very low and the local bond stresses have not reached their maximum values before the yielding of the steel reinforcements, this may help to explain why the load vs. mid-span displacement curves are quite similar between them.

RC beam flexurally tested by Carvalho et al. [24]
The numerical prediction of the load vs. mid-span displacement curve of the RC beam tested by Carvalho et al. [24] is also quite close to the experiment (see Fig. 20a). At the yielding instant, the IAE values calculated for the four different conditions never exceeded 2.00%. After 30 mm of mid-span displacement, the highest IAE value was associated with the model where the bond-slip relationship recommended by Model Code 2010 [14] was used (almost 4.00%). So, like in the RC beam tested by Biscaia et al. [25,26], similar results were obtained for the RC beam tested by Carvalho et al. [24], i.e. assuming (or not) a bond-slip relationship will produce no meaningful impacts on the load vs. mid-span displacement curve of the beam.
Looking at the principal maximum strain patterns of the four simulated conditions in Figs. 20b-e, it can be noticed that Bigaj's model [29] leads, once again, to a more generalized principal strain distribution along the span of the beam. In this case, the principal maximum strain value is approximately 1.13% which is located at the bottom of the beam and below the applied load (see Fig. 20e). Almost the same value (approximately 1.11%) and the location was predicted by the simulation of the RC beam with rigid contact between the steel reinforcement and concrete (Fig. 20b). In these two cases, the number of main cracks developed in the constant bending region is the same, i.e. five main cracks can be observed. At this yielding instant, the maximum bond stress developed between the steel reinforcement and the concrete is approximately 3.45 MPa and the maximum bond stress defined for the bond-slip relationship of 9.34 MPa was not reached so far.
The RC beams where the proposed new bond-slip relationship and the one recommended in Model Code 2010 [14] were used (see Fig. 20c, d, respectively), show some differences. For instance, the principal maximum strains in the RC beam where the bond-slip relationship recommended by Model Code 2010 [14] (see Fig. 20d) are mainly located at four cracks below the application of the loads and had a predicted value of approximately 0.09%. In the other RC beam where the proposed new bond-slip relationship was implemented, the principal maximum strain was mainly located at a single crack below the applied load in which a value of approximately 0.11% was predicted. The bond stresses developed in each model are also below the maximum values defined in each bond-slip relationship. In the first case, i.e. with the use of the bond-slip relationship recommended in Model Code 2010 [14], the maximum bond stress at the yielding instant of the beam is 3.73 MPa (τ bmax = 12.04 MPa) whilst when the proposed bond-slip relationship was implemented, the estimated maximum bond stress is 11.05 MPa (τ bmax = 11.67 MPa).
Once again, these results may help to explain why the four different numerical simulations of the RC beam tested  [14], and the Bigaj's model [29] are used, respectively by Carvalho et al. [24] are all quite close to each other as well as to the experimental data.

RC beam flexurally tested by Franco et al. [27, 28]
The load vs. mid-span displacement curve obtained from the RC beam experimentally tested by Franco et al. [27,28] is shown in Fig. 21a. The four simulations carried out in this work revealed, once again, that all are quite close to the experimental data and all produced also very similar results between them. In this case, however, the maximum load was under-predicted. Considering the failure of the RC beam at the yielding of the steel reinforcements, the experiment lead to a maximum load of 95.  [14], and 85.77 kN (− 10.0%) in the case of using the Bigaj's model [29].
The IAE values until the yielding of the beam are all quite small (less than 1.0%) but after this yielding point, it begins to increase rapidly due to the difference between the experimental yielding load and the predicted ones. The test conducted by Franco et al. [27,28] was purposely interrupted when the mid-span displacement was approximately 20 mm so it could allow us to use it for a future reinforcement of the beam with prestressed stainless steel and, therefore, after this displacement, no further IAE values were calculated and shown in Fig. 21a. Figure 21b-e show the principal maximum strains predicted by the numerical simulations at the yielding instant of the beam. In this case, the crack pattern predicted in the four numerical simulations is different. The model that used the bond-slip relationship recommended in Model Code 2010 [14] had the highest COD of approximately 0.43 mm while the lowest COD was obtained from the model that used the proposed new bond-slip relationship (approximately 0.39 mm). Three main cracks can be observed at the central region where the bending moment is constant when the bond-slip relationships in [14] and [29] were considered. In the other two cases, i.e. using the proposed new bond-slip relationship and considering a rigid contact between the steel reinforcements and concrete, four main cracks can be seen at that same region with a constant bending moment.
Considering the bond stresses developed between the steel reinforcements and the concrete, the results show that with the proposed new bond-slip relationship the maximum bond stress of 10.64 MPa is reached which is the maximum bond stress defined in the bond-slip relationship for the steel rebars with diameters of 8 and 12 mm. So, in this simulation, only negligible slips have developed between the steel reinforcements and concrete. With Model Code 2010 [14], the maximum bond stress was 3.52 MPa which is quite far away from τ bmax = 10.98 MPa whereas by using Bigaj's model [29] a maximum bond stress of 3.39 MPa is reached in the steel rebar with a diameter of 12 mm, which is 61.4% lower (τ bmax = 8.78 MPa).
So, considering all the aforementioned aspects, it can be noticed that until the yielding of the steel reinforcements, the differences pointed out previously were not sufficient enough to provide meaningful differences in the load vs. mid-span displacement curves of the RC beams. In addition, the predicted slips between the steel reinforcements and the concrete were very low in all three cases where a bond-slip was introduced to reflect the local adherence of the steel rebar-to-concrete interface. After this yielding point, the interfacial detachment of the steel reinforcements was Fig. 21 Numerical results of the RC beam tested by Franco et al. [27,28]: a load vs. mid-span displacement; b-e principal strains at the yielding instant when the rigid contact, the proposed new model, the Model Code 2010 [14], and the Bigaj's model [29] are used, respectively prevented and the loads applied to the RC beams remained similar in all the numerical simulations.

Conclusions
In the current work, a new nonlinear bond-slip relationship for ribbed steel rebars embedded in concrete was proposed whose calibration must be based on an experimental campaign. The proposed bond-slip relationship was implemented into a commercial finite element software and, for comparison purposes, two more bond-slip relationships were also considered. The influences of several parameters on the bond between ribbed steel rebars and concrete were studied and the differences obtained from the use of these three bond-slip relationships were analysed thoroughly. Based on all results achieved in this work, the following conclusions can be made: • Due to its smooth transition between the maximum and residual bond stresses, the proposed new bond-slip relationship showed to be more precise with the experimental results than any of the other bond-slip relationships. Its versatility was shown by the precision achieved whether the specimens had short or long embedded lengths and consistency between the local and global behaviour was demonstrated; • The simulations showed that the increase of the embedded length increases the load capacity of the specimens. However, due to the yielding and subsequent hardening of the steel rebars, the load capacity of the test is limited by the maximum load stress of the steel rebars. When considering a perfect elastic stress-strain of the steel rebars, Bigaj's model [29] is the unique one to slow down the increase of the maximum load with the embedded length, which the assumption of non-zero bond stress after reaching the final slip can be pointed out for that behaviour; • Since the increase of the concrete strength increases both the maximum and the residual bond stresses, the load capacity of the specimens is also increased. However, it should be noted that such increases may be limited by the yielding of the steel rebar which was more evident when the bond-slip relationship recommended in Model Code [14] was used. For the analyzed cases, using Bigaj's model [29] the yielding of the steel rebars was never reached which may be explained by the lack of a friction stage in this bond-slip relationship; • Considering the influence of the cross-sectional area of the concrete blocks, the use of the three bond-slip relationships led to a consensual conclusion, i.e. at least a 200 × 200mm 2 cross-sectional area of concrete should be used. With such a dimension, the splitting of concrete can be prevented and an interfacial failure mode or the rupture of the steel rebar can be obtained; • The yielding of the steel rebars is a parameter with high influence on the pull-out tests and, consequently, on the definition of the local bond-slip relationship of the steel rebar-to-concrete joint. Moreover, as the embedded length increased, the influence of the shape of the bondslip relationships on the load-displacement response of the specimens becomes almost marginal due, precisely, to the yielding of the steel rebars; • Since the embedded area increases with the increase of the diameter of the steel rebars, it was not surprising to see that the load capacity of the specimens increases with this parameter. However, to determine some parameters of the bond-slip relationship between the steel rebar and the concrete such as the residual bond stress or the transition slip point between the maximum and residual bond stresses, the yielding of the steel rebars should be avoided; • The definition of the exact value of the exponential parameter α in the bond-slip relationship in [14] produces different P-δ 0 . For instance, when a zero value is assumed, the P-δ 0 curves stay close to the results obtained from the proposed new bond-slip relationship. However, if the yielding of the steel rebar is ignored, huge differences between both can be found. With α = 1, the initial branch of the P-δ 0 curve is also linear and for the same bond stress it provides higher displacements; • The influence of each bond-slip relationship on the surrounding concrete estimated in the numerical simulations is similar when a short embedded length is considered while some differences between them can be found when a larger embedded length is assumed. Thus, the corresponding stress bulb shapes developed in the three models with short embedded lengths are very similar to their corresponding maximum load capacities. In these, the embedded length was fully mobilized which means that an interfacial failure mode was always predicted. However, the maximum loads reached in each simulation occurred at different displacements, being in the proposed new bond-slip relationship the one that led to the smallest displacement, against the model with the bond-slip relationship recommended in [14] where the highest displacement was predicted. With a long embedded length, the proposed new bond-slip relationship concentrates the bond stresses at, approximately on the first part of the steel rebar, whilst the other models predicted a wider range of the bond stress distribution throughout the embedded length. Such difference can be explained by the lack of an E stage in the proposed new model, which should be a subject to look at in the future; • With a calculated IAE of less than 5% until a mid-span displacement of 30 mm in the RC beams tested by Bis-caia et al. [25,26] and Carvalho et al. [24] and 20 mm in the RC beam tested by Franco et al. [27,28], the numerical simulations, independently of the bond-slip relationship used, were very close to the experimental load vs. mid-span displacement curve. The results showed, therefore, that the type of bond-slip relationship did not influence the load vs. mid-span displacement curve of the RC beams herein studied. The main influence is the crack pattern of the RC beams where the use of Bigaj's model [29] conducted into a more generalized crack pattern in the beams. Moreover, the RC beams that used Bigaj's model [29] developed always the highest principal maximum strains.
Finally, it should be recognized here that for wider use of the proposed new bond-slip relationship, it would be needed to refine it by considering other material properties such as those parametrically analysed in this work or others, e.g. diameter of the steel rebars, concrete strength, type of concrete (e.g. with fibres or recycled aggregates), dimensions of the transversal steel ribs, etc. For that achievement, the present experimental program will be extended in the future so it could be possible to refine the values of a and s 0 determined here and adjust the proposed new bond-slip relationship to a wider range of applications.

Appendix A
This appendix presents the numerical models that were used to simulate the pull-out test carried out during this work as can be seen in Table 4. To facilitate the identification of each model and its corresponding characteristics, similar nomenclature given to the tested specimens (see Sect. 2.4) was adopted here. However, instead of using "g" as the number of repeated tests, "α", "α0", "α0.4" or "α1" is used instead, making it easier to identify the models where this parameter was not considered, was set equal to 0, 0.4 and 1, respectively. Also, letter the "a" stands for the type of bond-slip relationship used in the simulations, i.e. "PN" (Proposed New bond-slip relationship), "MC" (Model Code 2010 [14]) or "B" (Bigaj's model [29]). As an example, the model simulated with an embedded length of 80 mm, with concrete strength of 37.2 MPa, a cross-sectional area of concrete with 200 × 200 mm 2 , and yielding stress of the steel rebar equal to 530.5 MPa, corresponds to the model PN-L80-f37.2-Y530.5-D200-ϕ8-α, in which parameter α was not used since it is an exclusive parameter of the bond-slip relationship recommended in Model Code 2010 [14]. In total, 139 different numerical models were idealized and simulated with the 2D version of ATENA v. 5.9.1 [36].