Numerical assessment of bond-slip relationships for naturally corroded plain reinforcement bars in concrete beams

Abstract Reinforced Concrete (RC) heritage structures are often affected by corrosion. Consequently, knowledge about the effect of corrosion on the bond between reinforcing bars and surrounding concrete is critical when assessing the structural performance of these structures. In earlier work, structural tests were carried out on segments of edge beams taken from a decommissioned RC bridge. The specimens had naturally corroded plain reinforcement bars and three-point bending tests were conducted, to investigate their anchorage capacity. In this study, non-linear finite element analyses (NLFEA) were carried out to gain further insight into the bond behaviour of the tested specimens, including the effect of corrosion on the bond-slip relationship. Two different, one-dimensional (1D), bond-slip relationships were calibrated for each tested bar, to account for loss of bond upon yielding. The calibration process was based on a comparison between significant numerical and experimental results, including load–deflection curve, crack pattern and asymmetrical distribution of the yield penetration along the length of the bar. Good agreement between the FE analyses and experimental tests was observed. Finally, the calibrated bond-slip relationships for nine beams with different corrosion levels, casting positions, and visual damage are presented and discussed. The loss of bond at yielding and yield penetration asymmetry are both shown to be crucial factors for adequately describing structural behaviour.


Introduction
Corrosion of reinforcing bars is a major cause of deterioration in reinforced concrete (RC) structures. Corrosion reduces the crosssectional area of the bar, leads to concrete cracking and cover spalling, and affects the bond between bar and concrete [1]. As corrosion propagates, the serviceability and safety of RC structures are seriously affected. Since many existing structures are deteriorated due to reinforcement corrosion [2], reliable methods of assessing the performance of these structures are needed urgently.
Numerous studies have been conducted on the effect of corrosion on the structural performance of RC beams. Based on experimental results [3][4][5][6], by decreasing both ductility and reinforcement area [7], corrosion may lead to alterations in the failure mode, (from flexure to shear, shear-compression or anchorage, for example) and reduce the loadcarrying capacity of the structure. However, it should be noted that the majority of studies focused on ribbed reinforcing bars, as these had become the norm in construction from around the mid-1950s in America and Canada and around the mid-1960s in Europe [8]. Even so, plain bars are still widely present in historical RC structures which, in many cases, now require assessment and rehabilitation [8]. Only a few studies of beams which include corroded plain bars are available in the current state-of-the-art [9][10][11] and the behaviour they describe differs significantly, as compared to corroded RC members with ribbed bars.
One important underlying aspect is the difference in bond mechanisms between ribbed and plain bars. This means there is an anticipated difference in the impact of corrosion on their bond mechanisms. The main difference between plain and deformed bars is the lack of ribs. The mechanical interlock between ribs and concrete is the major component in the bond of deformed bars, while the bond in plain bars fundamentally depends on the adhesion and friction of the bar-concrete interface [12]. The different bond mechanisms indicate that the results obtained from deformed bars cannot be directly extended to plain bars [8,13].
A limited number of bond-slip relationship for plain bars can be found in codes and literature. In Fig. 1, three examples are presented. The bond-slip relationship for plain bars in Model Code 2010 [34] is calculated as a function of the concrete compressive strength expressed in MPa. Good and poor bond conditions (depending on the placement and confinement of the bar) are distinguished, and a reduction factor of 0.5 is suggested for poor conditions. The bond-slip relationship is described by an ascending non-linear curve up to the maximum bond stress, followed by a constant branch. The slip corresponding to the maximum bond stress is constant and equal to 0.1 mm. Other bond models can be found in the literature for uncorroded plain bars, most commonly based on the results of pull-out tests. Verderame et al. [38] proposed a bond-slip model with as well an initial non-linear brunch, but with a larger slip value at peak; 0.23 mm. The model introduces a linear descending branch for the after-peak behavior, followed by a plateau. The compressive strength of the concrete remains the main parameter for the definition of the maximum bond strength. A third bond-slip constitutive model worth mentioning is by Feldman and Barlett [10]. Their empirical model differs from the others in the use of a non-linear descending branch and in including parameters in addition to the compressive strength for the calculation of the bond strength: the roughness of the bar (R y ) and the embedded length (l d ).
However, there is a dearth of available research on the bond behaviour of corroded plain bars [13][14][15][16]. Such studies indicate relevant parameters (such as corrosion level, presence of confinement, cover thickness, bond at yielding and the initial condition of the bars) that influence the bond behaviour. Also, the casting position (top or bottomcast) was deemed a highly significant factor due to the varying density of the concrete surrounding the bar and, thus, the possibility of corrosion products accumulating in the cement pores [11,13]. Additionally, there are concerns about using results obtained from artificially corroded specimens to predict the behaviour of naturally corroded bars [17][18][19]39].
Modelling, specifically nonlinear finite-element analyses (NLFEA), is a widely used tool [20,21] to improve our understanding of how the bond between concrete and naturally corroded plain bars is affected by corrosion products. Detailed bond and corrosion models are often implemented at interface elements between reinforcement and concrete, in three-dimensional (3D) NLFEA [22,23]. However, considering the computational cost and amount of data needed, such an approach may be impractical for assessing entire structural systems. A pre-defined, one-dimensional (1D) bond-slip relationship in 3D NLFEA has already been proved as a valid alternative; it provides a reliable and timeefficient description of structural behaviour [21]. This method was chosen for this study, to examine the structural behaviour of beams with naturally corroded plain bars.
In reality, the beams were cut from the edge beams of a decommissioned bridge in Sweden, that had been previously tested by the authors [11]. For most of the beams, yielding of the reinforcement bars limited the load-carrying capacity. Yielding was followed by anchorage failure and, thus, the deformation capacity of the beams was limited. The loss of bond after yielding led to a change in the load-carrying mechanism, from beam to arch action.
This study aims to improve understanding of the structural behaviour of the tested beams; more specifically, how the bond-slip relationship of plain reinforcement bars is affected by yielding and corrosion damage and how this, in turn, affects beam behaviour. Consequently, 3D NLFEA was conducted on the tested beams, including individually calibrated bond-slip curves for each reinforcing bar. To account for the effect of yielding on bond behaviour, the tensile bars are divided into yielded and unyielded zones in the FE model, with distinct bond-slip relationships assigned. By comparing simulation results with experimental measurements, the NLFE models (and, thereby, the 1D bond-slip relationships) are calibrated and validated for each bar in the nine tested beams. Finally, the overall numerical results and calibrated bond-slip relationships are discussed.

Experimental program and materials
Gullspång Bridge (1935-2016) was a concrete bridge reinforced with plain steel bars and located in the town of Gullspång, Sweden. The Fig. 1. Bond-slip relationships for plain bars embedded in concrete. bridge was exposed to different conditions during its service life, including freeze-thaw cycles and the use of de-icing salts. Heavy corrosion damage was visible in parts of the structure. The edge beams of the bridge were extracted before demolition and divided into smaller members for research purposes. The experimental study comprised two different investigations into the effect of corrosion on the bond of plain reinforcement bars; one focusing on structural behaviour [11] and one focusing on local bond strength [16].
During the structural investigation, twenty 900 mm beams were cut from the edge beams and classified according to the extent of visible corrosion damage within the anchorage zones. Three categories were defined: (1) reference beams (R), meaning beams with no visible cracks, (2) cracked beams (C), meaning beams with visible surface cracks, and (3) severely damage beams (S), meaning beams with corrosion-induced cracks and concrete spalling. The results of structural tests conducted on these members are presented in Robuschi et al. [11]. Nine out of twenty tested beams were selected and studied in this work.

Three-point bending tests
A three-point bending test set-up was used to bring the specimens to failure, see Fig. 2. The test configuration was designed to evaluate the influence of corrosion on the anchorage capacity of the tensile bars [24]. A narrow support (NS, 50 mm × 100 mm), was placed on one end of the beams, to minimise the confinement effect on the tensile reinforcements due to support pressure. A full support (FS, 50 mm × 250 mm), was used at the other end of the beam, to prevent the beam from tipping during loading. To assess the influence of the casting position of the bars (in other words, top or bottom-cast bars), the beams were positioned in two different configurations in the rig; as positioned on the bridge or upside down. For the present study, two beams were chosen as positioned on the bridge and seven upside-down. The reasons for testing more upside down beams were twofold; (1) because top-cast bars are known to have a weaker bond between concrete and steel and were therefore deemed more interesting and (2) because severe corrosion damage occurred, with bars falling off in the bottom inner position. Thus, there were not so many relevant results available for beams tested as positioned on the bridge.

Yield penetration and corrosion level of tensile reinforcement bars
After testing, the tensile reinforcement bars were carefully taken out of the beams. Corrosion products and concrete attached to the bar surface were removed by sandblasting, according to [25]. The bars were then 3D-scanned. This allowed the yielded parts of the bars to be located and their average corrosion level to be estimated, by calculating the average loss of cross-sectional area in the unyielded segments according to [11]. Geometrical information (such as the uncorroded area and perimeter) was also obtained for each bar by 3D scanning.

Results of the bending tests
The tests were characterized by the opening of one or two major bending cracks, that was closely followed by yielding of the tensile reinforcements. In most of the tested beams, their load-carrying capacity was limited by yielding of the bars. After yielding, slipping of the two tensile reinforcement bars occurred at different deflections, as shown in Fig. 3. However, the failure was defined as a bending failure, since slipping of the bars occurred at a point where the bending cracks had already reached an opening width of several centimetres. Of the nine modelled beams, only two failed in anchorage before the tensile reinforcement yielded. This study numerically models the behavior of the beams in 3-point bending and makes use of the experimental measurements both as input for characterizing the bond of the naturally corroded, plain reinforcing bars and for the calibration of the FE models.
Three parameters were of particular importance: mid-span deflection at end-slip, loss of load capacity at end-slip and yield penetration along the length of the bar. The mid-span deflection and loss of load capacity at end-slip were critical for the calibration of bond-slip relationship (as discussed in Section 4). The yield penetration was preassigned to allowed for the use of simple, 1D bond-slip relationships, as explained in the following. The corrosion level of the rebars was as well pre-assigned in the model, while the crack pattern and crack widths were used in the calibration procedure.
The need for pre-assigning the yield penetration in the model was due to the asymmetric yield penetration along the length of the bars observed in the tests. This caused the unyielded length at the side were anchorage failure took place being significantly shorter than at the other side of the beam, see Fig. 4. This behaviour was observed for most of the tested beams, despite the loads being equidistant from the supports. Factors such as heterogeneity of the material properties of the concrete, non-uniform corrosion along the bar, unequal stirrups disposition and use of different support types (full and narrow support), are likely linked to the asymmetrical distribution of the yield penetration. To evaluate this asymmetry, an asymmetry factor (a) was defined as: Fig. 2. Test set-up and view of narrow support in the beam tests [11]. The beam shown was tested as positioned on the bridge. Measurements in mm.
where Δl is the difference between the centre of the bar and the centre of the yielded region, and l is the total bar length, Fig. 4.

Properties of concrete used in the analyses
Material tests were conducted on concrete cylinders drilled from the edge beams. A total of eight uniaxial compression tests were conducted to determine the compressive strength (f c ) [11]. The modulus of elasticity (E cm ) was obtained from three uniaxial compression tests. Additionally, splitting tensile tests were conducted on four additional cylindrical samples, according to EN 12390-6:2009 [26]. The fracture energy (G f ) was calculated from the compressive strength according to CEB-FIP Mode Code 90 [27], assuming the maximum aggregate size to be 16 mm. The concrete material parameters used in the analyses are provided in Table 1.

Tensile properties of reinforcement bars
Tensile tests were carried out on bar segments cut from the unyielded zones of each tensile reinforcement. Detailed information on the tensile tests is provided in [11]. The average corroded cross-sectional area of the corresponding segment was used to calculate the tensile stress. The stress-strain curve obtained from the tensile test of each bar was used as input to the corresponding reinforcement element in the model. For five of the bars, no information was available; due to damages resulting from the extraction process, they could not be tested. However, all the bars in the beams had similar properties, and similar load-displacement curves [16]. Therefore, the yield stress in these bars was assumed equal to the average yield stress of uncorroded tested bars (see Fig. 5). This is based on the assumption that the yield stress does not change with corrosion at low corrosion levels; this is consistent with observations in several studies [36,37]. The main properties of each curve, such as yield stress (f y ), yield strain (ε y ), ultimate stress (f u ) and ultimate strain (ε u ) are listed in Table 2. The bars are labelled according to their position in the crosssection of the beam; TI (top-cast, inner side of the bridge), TO (top-cast, outer side of the bridge), BO (bottom-cast, outer side of the bridge) and BI (bottom-cast, inner side of the bridge) (see Fig. 2).

Numerical modelling
This section introduces the 3D NLFE analyses of the presented beams. The analyses were used to simulate their structural performance, with 1D bond-slip relationships used to represent the bond between the naturally corroded plain bars and the concrete.

FE model description
The finite element software DIANA 10.2 [28] was used for the analyses. The geometry of each model was based on the specific measurements of the corresponding specimen. An overview of a 3D model is shown in Fig. 6. Spalling of concrete cover was included in the model by removing the concrete surrounding the bar at the spalled locations, as shown in Fig. 6b. Corrosion of the reinforcement was implemented as an average cross-sectional reduction of reinforcement area in the unyielded segments. Corrosion was not considered in the cross-section of the yielded segment of the reinforcement bar; due to the deformation induced by yielding, it was not possible to calculate the average corrosion level of the segments with 3D scanning data. The influence of corrosion damages on the yielding behaviour was however observed to    be negligible; likely due to the relatively small corrosion levels, 1-2%. 3D tetrahedral elements (TE12L) were used for the concrete. The longitudinal tensile reinforcement bars were modelled using embedded reinforcement with 1D bond-slip relationships. The bond-slip relationships were described by multi-linear curves. By contrast, the compressive reinforcement bars and stirrups were modelled as fully bonded and using embedded reinforcement elements, see Fig. 6c. The loading plate and supports consisted of a steel plate and a wood one (as in the tests) and were modelled using triangular-prism elements (TP18L). The geometrical dimensions of support and loading plates were also the same as in the test set-up. The full support (FS) was restrained in the vertical direction, while the narrow support (NS) was restrained in both vertical and longitudinal directions. The load was applied via displacement control on a master node at the top of the steel loading plate. A line of nodes was tied to the master node; this allowed rotation of the plate around the y-axis but imposed the same displacement in z-direction.
An incremental static analysis was conducted, using specified displacement increments. These increments were set to a vertical displacement (z-axis) of 0.01 mm until a 1 mm mid-span deflection was reached. The increments were then increased to 0.2 mm to save computation time. The analyses were conducted using a Quasi-Newton (BFGS) scheme, based on force and energy convergence criteria, with respective tolerances of 0.01 and 0.0001.

Material models for concrete and steel reinforcement
For concrete, a total strain-based rotating crack model (using nonlinear fracture mechanics) was chosen to describe cracking. The tensile behaviour was modelled using a stress-strain curve proposed by Hordijk [29]. The crack band width (h b ) was assumed as equal to the average element size of 0.02 m. This assumed size of the localization zone was later verified in the analysis results. The compressive behaviour of concrete was described by an ideal elasto-plastic relationship. The reduction of compressive strength due to lateral cracking was taken into account using a model developed by Vecchio and Collins [30]. The stress confinement was also taken into account, as suggested by Selby and Vecchio [31]. The concrete properties used in the analyses are listed in Table 1.
The steel reinforcement was modelled using a Von Mises plasticity model, including strain hardening. A simplified multilinear curve obtained from tensile testing of each individual bar was used as input. The tensile properties of the reinforcement bars are provided in Table 2.

1D bond-slip relations
As shown in different studies, bond is strongly affected by bar yielding [32,33]. This phenomenon is even more relevant with plain bars than ribbed ones, as a major part of the bond capacity comes from friction. Since the experimental tests were significantly affected by yielding of the tensile reinforcement bars, the loss of bond at yielding needed to be included in the modelling. Thus, the bars were divided into different zones according to test measurements, as shown in Fig. 6c. Also, two different bond-slip relationships were defined for each bar in the analyses; one for the yielded zone and another for the unyielded zones.
The so-called yielded zone was defined over the length of the bar affected by yield penetration, as measured after the tests. However, it should be noted that assigning a single bond-slip relationship to all sections of the bars that yielded during testing is a simplification. Thus, the bond-slip relationship given to the yielded zone does not correspond to the bond-slip curve of a yielded bar. Rather, the bond-slip curve of the yielded zone is designed to describe the bond at the beginning of the test (when the bar was unaffected by yielding) and the subsequent loss of bond associated with yielding. Note that the tensile properties was the same in the yielded and unyielded zones. Thus, the bond-slip relationship was the only parameter changed. It can also be noted that the same bond-slip relationship can be expected in both regions before the reinforcement yields; the chosen division of different bond-slip relationships in yielded and unyielded zones was made to simplify the analyses. Fig. 7 describes the bond-slip relationships used for the yielded and unyielded zones. Such relationships were defined by nine different parameters; five for the unyielded bond-slip and four for the yielded one. The nine parameters were divided into three sets (highlighted in different colours in Fig. 7), due to different roles in the calibration procedure (details in Section 4). The bond-slip relationship for the unyielded zones resembles the one proposed in CEB-FIP Mode Code 2010 [34], which is characterised by four branches; ascending, plateau, descending and residual. The following parameters define the transition between branches: τ max and τ f are, respectively, maximum bond stress and residual bond stress for unyielded zones; s 1 and s 2 are the slips at the start and end of the maximum bond stress plateau. Finally, s 3 corresponds to the slip when residual bond stress is reached. However, the bond-slip relationship for the yielded zone consists of only three branches; ascending, descending and residual. Thus, no plateau is included for the yielded zone. The main parameters defining the transition between branches are τ ymax and τ yf (representing maximum and residual bond stresses respectively) and s y1 and s y2 (corresponding to the slip at maximum bond stress and slip at the beginning of the residual bond stress branch respectively). In addition, at the steel concrete interface the normal stiffness modulus was set to 2 × 10 14 N/m 3 , while the shear stiffness modulus was calculated as τ max /s 1 , τ ymax /s y1 .

Calibration procedure of bond-slip parameters
As explained in the previous section, nine unknown parameters characterise the bond-slip relationships in a single bar. This section presents the calibration process for these parameters. Assumptions and interpretations about the influence of each parameter on the structural behaviour are described, plus the calibration methodology used. As an example of the parameter calibration process, a thorough description is given of the entire process for a single tensile reinforcement bar; namely TO in beam 17F. Note that, once the manually calculated, fixed parameters are defined for both bars, the calibration for a single bar in the beam becomes essentially independent of the other bar. To conclude, the model and, consequently, the calibration of the bond relationships, is validated against the experimental results by thoroughly comparing various aspects; load-deflection, crack pattern, crack width, yield penetration and yield penetration asymmetry.
As mentioned earlier, the nine parameters representing the bond-slip relationship of each yielded tensile bar were grouped into three sets and highlighted in different colours, Fig. 7. The first set contains parameters chosen to have the same values for all bars. This includes τ ymax , τ yf , s y1 and s 1 . The second set contains two parameters manually calculated from test results, τ max and s 3 . The third set consists of three parameters, τ f , s y2 and s 2 ; these were calibrated by comparing the results of the FE analyses with the experimental results for each beam individually. The bars that slipped before yielding only required determination of five parameters, since only the bond-slip relationship for the unyielded zone was used. A parametric study was undertaken, to understand the influence of each parameter on the structural behaviour of the modelled beam. An overview of the results appears in Table 3.

Constant parameters
The parameters τ ymax , τ yf , s y1 and s 1 were chosen to have constant values for all bars, based on similarities in their influence on the anchorage behaviour of the various beams. This is described and discussed in more detail for each parameter below.

Maximum bond stress for yielded zone, τ ymax
The maximum bond stress for the yielded zone, τ ymax , describes the maximum bond reached before yielding took place in the predefined yielded zone. In reality, this value likely varies along the bar in the yielded zone. In areas close to the crack, where the bar first started yielding, the maximum bond stress was likely lower than in areas further from the crack. However, by way of simplification, this study chose a single value for the entire yielded zone. Upper and lower bounds were used to give trial values. Experimental results (in the form of a known force transfer along a known length of the bar when yielded), gave an estimate of 1 MPa [11]. This was deemed a lower bound, as it was affected by yielding. The upper limit was evaluated based on the bar reaching the yield force in the cracked section when the bar first started yielding. Thus, the average bond stress (acting along the distance between the first bending crack and the end of the beam) had to carry the yield force. This corresponded to a value of between 2.5 and 3 MPa in the different tested beams.
Accordingly, three τ ymax values between these boundaries were investigated in the parametric study; 1, 2 and 3 MPa. It was observed that τ ymax largely affected crack width and mid-span deflection (Δ m ) at end-slip. Fig. 8 shows the crack pattern and crack width in the various analyses. As may be observed in Fig. 8d, the results for τ ymax equal to 2 or 3 MPa showed good agreement for crack width, while a value of 1 MPa underestimated it. Further, in Table 3, it may be seen that τ ymax equal to 1 or 2 MPa gave reasonable values for the mid-span deflection at endslip, while τ ymax equal to 3 MPa did not. Accordingly, τ ymax at 2 MPa was chosen, as this value gave the best agreement for both crack width and mid-span deflection at end-slip. As most of the beams showed similar crack patterns (one major crack  accompanied by a minor crack) and because the main focus was on reproducing the anchorage behaviour, this calibration was deemed sufficient. Thus, the same τ ymax value of 2 MPa was used in all the models.

Residual bond stress for yielded zone, τ yf
The residual bond strength of a yielded bar, τ yf , is expected to be very small [11,32,33,35]. τ yf values varying from 0.2 to 0.6 MPa were investigated. From the FE results, it was found that τ yf equal to 0.2 or 0.4 MPa resulted in the anticipated failure mode (end-slip of the bar after yielding). The failure mode turned to concrete crushing when τ yf, was set to 0.6 MPa, without bar slippage. In this case, for comparison, the yield penetration and asymmetry values were measured just before concrete crushing.
The tensile stress distribution along the bar was also affected by the chosen τ yf value. Fig. 9 describes the stress distributions of the bar at a mid-span deflection (Δ m ) of 19.2 mm (when the bar started slipping on the side with the narrow support) for different τ yf values. From the results, it was concluded that τ yf equal to 0.2 MPa produced better agreement in both yield penetration and stress distribution. High τ yf values could transfer higher bond stresses, which shortened the yield penetration length and decreased the arch action in the mid-span. Thus, even though the bar was divided into yielded and unyielded zones from the beginning, the length and asymmetry of the yield penetration could only be described correctly when a small τ yf value was applied. Given that a similar behaviour was observed in all the tested beams where yielding of the tensile reinforcements took place, τ yf at 0.2 MPa was adopted for all the modelled beams.

Slip corresponding to maximum bond stress for yielded zone, s y1
The slip corresponding to maximum bond stress for the yielded zone was expected to be comparatively small, since the bars start yielding right after the opening of the bending crack, when the slip between bar and surrounding concrete is still limited. s y1 values ranging between 0.005 and 0.1 mm were investigated. The results in Fig. 10 and Table 3 show that s y1 mainly affected the crack width. Large s y1 values influenced the mid-span deflection at end-slip but to a lesser extent. The crack width was the main indicator of correct calibration of the parameter. As seen in Fig. 10, the crack width decreased as the s y1 value increased and the best agreement with experimental measurements was found for an s y1 value of 0.01 mm. Since the manner of yielding was similar for all the tested bars, s y1 = 0.01 mm was chosen as constant in this study.

Slip at the start of maximum bond stress for unyielded zones, s 1
According to [34], the slip corresponding to maximum bond stress for the unyielded zone, s 1 , is suggested to be around 0.1 mm. However, the code proposes an exponential curve for the initial branch of the bond-slip relationship. Since this work used a linear relationship instead, a lesser s 1 value was deemed a reasonable compensation. Therefore, the investigated s 1 values ranged between 0.01 and 0.1 mm.
As may be seen in Table 3, the numerical results for s 1 varied between   0.01 and 0.1 mm and showed almost no influence on the main structural parameters. In other words, mid-span deflection at end-slip, residual load capacity and yield penetration and asymmetry remained virtually constant, regardless of the value chosen. There was some influence on the crack width, see Fig. 11. With the increase in s 1 , the bond stiffness decreased and thus the crack width increased. s 1 was adopted at 0.01 mm; the same value as s y1 in this study. However, for the reasons outlined at the beginning of this section, it is noteworthy that the peak real slip of a plain reinforcement bar is probably slightly greater than this value [16].

Manually calculated parameters
The parameters τ max and s 3 were chosen to be manually calculated from test results, by assuming uniformly distributed bond stress in the unyielded zone on the end-slip side. This is described more in detail below. By calculating these from experimental results, it is possible to include the varying damage levels (corrosion and cracking) of the beams.

Maximum bond stress for unyielded zones, τ max
The assumption of uniformly distributed bond stress was used in the previous study to calculate the average bond strength in the unyielded zone at the end-slip side [11]. Although such assumptions are commonly used when calculating bond stress in short pull-out tests, the actual bond stress distribution before anchorage failure is unknown and may significantly deviate from uniformity. Fig. 12 depicts three examples of simplified bond distributions along the unyielded zone. All three have the same area under the curve, but they yield to different maximum bond stresses. For the sake of safety and simplification, this study adopted τ max as the average bond strength, according to the short pullout hypothesis. This corresponded to the lower bound of possible maximum bond stresses.

Slip corresponding to residual bond stress for unyielded zones, s 3
The slip at the start of residual bond stress for unyielded zones (s 3 and s 2 ) defines the descending slope of bond-slip curves for unyielded zones. It may be noted that, in the case of a large difference between s 3 and s 2 , the bond distribution along the unyielded zone is closer to the piecewise linear bond stress shown in Fig. 12. In keeping with the assumption of uniformly distributed bond stress in the unyielded zone at the end-slip side, a small difference between s 3 and s 2 was chosen (a very stiff descending branch). Accordingly, the initial choice was to have the same stiffness for the descending slope as for the ascending slope, but with the opposite sign. This choice correlated with an s 3 value of 0.104 mm. As shown in Table 3, the influence on the results was small when s 3 varied between 0.102 mm and 0.108 mm. (Note that although the variation in s 3 may be deemed small, the variation in stiffness of the descending slope is large. This varies by a factor of four, between the largest and smallest chosen inputs.) Accordingly, s 3 was calculated to give the same stiffness for the descending slope as for the ascending slope, but with the opposite sign.

Individually calibrated parameters
The parameters τ f , s y2 and s 2 were calibrated for each beam. The most significant results of the FE analyses for each are described below.

Residual bond stress for unyielded zone, τ f
The residual bond stress for unyielded zones, τ f , determined the remaining load capacity after the end-slip of the corresponding bar. As stated above, in most of the beam tests, slipping of the two tensile bars did not start simultaneously, but rather one after the other. This caused two decreasing stages in the load-deflection curve of the beam: (1) When the first bar started to slip, the load capacity decreased and the residual load capacity basically relied on the residual bond stress of that bar and maximum bond stress of the other bar and (2) after slipping of the second bar, the residual load capacity relied on the residual bond stress of both bars, see Table 3. Therefore, the first bar to slip was calibrated first.

Slip corresponding to residual bond stress for yielded zone, s y2
The slip corresponding to residual bond stress for the yielded zone, s y2 , had a critical influence on the mid-span deflection at end-slip of the bar. Increasing s y2 from 0.25 mm to 0.45 mm delayed slipping of the bar (see Table 3); in other words, the beam deflection at end-slip increased. This is because the greater bond stresses delayed the anchorage failure.

Slip at the end of maximum bond stress for unyielded zones, s 2
As may be seen in Table 3, the slip at the end of maximum bond stress for unyielded zones, s 2 , governed the mid-span deflection of the beam at end-slip significantly. Large s 2 values resulted in a longer plateau in the load-deflection curve. This is reasonable, since increasing s 2 values meant that the bar would maintain high local bond values in relatively large slips.

Verification of calibrated parameters
As an example of what was compared and verified for all analysed beams, the results of the FE analyses of beam 17F (with calibrated parameters) are presented below and compared with the experimental results.

Load versus mid-span deflection and crack pattern
The load versus mid-span deflection from testing and FE analyses of beam 17F is shown in Fig. 13. In both test and analysis, the beam behaviour was governed by the appearance of a bending crack at midspan followed by yielding of the tensile reinforcement bars. Thereafter, the applied load remained approximately constant, up to a midspan deflection of 13.3 mm in testing and 13.6 mm in FE analysis (when the top inner (TI) bar at the full support (FS) side slipped). At a deflection of 19.1 mm (test) and 19.2 mm (FE analysis), the top outer (TO) bar at the narrow support (NS) also slipped. Note that the bars are named based on their original position on the bridge and that beam 17F was tested upside down. As may be seen in Fig. 13, the agreement between the analysis and the experiment is very good. The model was able to capture the overall beam response during all stages, including cracking, maximum load, remaining load and deflection when end-slip of the various bars occurred. Fig. 14 shows the crack pattern in both the outer and inner sides of beam 17F, as recorded at the end of the test (Δ m = 23.3 mm). The crack pattern on the outer side of the beam was captured using digital image correlation (DIC), while that on the inner side was measured after testing. On the outer side, one major bending crack was located at approximately mid-span, while a minor crack may be observed to have developed on the left side. Only one bending crack formed on the inner side. Good agreement was shown between the FE results and the observed crack pattern. The crack width was also in good agreement, as shown earlier in Figs. 8, 10 and 11. Note that the crack width compared in these figures corresponded to the crack width of the major bending crack, as measured by DIC on the outer side of the beam at reinforcement level. The values measured for the crack width were 31.5 mm (outer side) and 30.7 mm (inner side). The FE analysis results were similar to the experimental ones, at 31.2 mm (outer side) and 28.3 mm (inner side).

Bar stress and bond behaviour
Tensile stress, bond stress and slip along the length of the studied bar at three different loading stages (cracking, bar yielding and failure) are shown in Figs. 15 and 16. The yield penetration in the FE analyses increased from the onset of yielding until slip occurred. At this point, the yield penetration in the analyses was in good agreement with the experimental measurements of length and position.
The changes in bond behaviour along the length of the bar aligned well with the changes in its tensile stress distribution at the various loading points. At cracking, good bond conditions allowed for transfer of stresses between the reinforcement and its surrounding concrete. Thus, the stress in the bar was symmetrical, decreasing with the distance from the centre of the bar (Fig. 15) and following the bending moment distribution. The bond stress did not reach the maximum bond stress in either the yielded, τ ymax , or unyielded, τ max , zones. The opening of the bending crack is clearly visible in Fig. 16a, as characterised by an initial bar slip at the crack location.
Under increased loading, the main crack opened, leading to increased slip at the crack location. Simultaneously, the bond stress reached the residual values, τ yf and τ f in some parts. Note that abrupt changes in bond stress distribution occurred at the transition between yielded and unyielded zones, as these regions were assigned different bond-slip relationships in the model, see Fig. 16b. These abrupt changes may not be fully realistic as, at this point of loading, the bar was not yet yielding along the entire length assigned to the yielded zone. Thus, this behaviour is a direct consequence of the modelling choice of a predefined yielded zone. Avoiding this behaviour would require a far more complex model. Fig. 16c shows slip and bond stress distribution in the load step just before end-slip occurred. At this load step, it may be observed that the bond along the whole yielded zone had reached the residual value, τ yf and thar the bar is fully yielded in the predefined region (Fig. 15). The loss of bond in the mid-region of the beam, the yielded part, indicated a transition in structural behaviour from beam to arch action. This may also be seen in Fig. 15, where the bar stress in the yielded zone was almost constant.
The transition from beam to arch action was previously observed in the experiments, [11], and has been reported by other researchers for similar cases. Feldman and Bartlett [10] recorded the variation in bar tension force and bond stress at different load levels, for a beam reinforced with plain bars. They observed near-complete bond loss in the central region of the beam, corresponding to the yield penetration. With increasing applied load, increased arch action within the mid-span brought the location of the calculated peak bond stress towards the support reactions, ultimately causing anchorage failure. Furthermore, in Coronelli and Gambarova's [20] modelling of simply supported corroded reinforced concrete beams, there was a complete bond loss in the midspan and stress redistribution at the ultimate limit state; this caused an arch action in the beam.
To conclude, the above comparisons showed that the FE analyses of 1D bond-slip relationships were capable of reasonable representation of beam behaviour. In particular, the changed load-carrying mechanism from beam to arch action and the asymmetric development of the bar yield penetration were clearly captured. This showed that the adopted bond-slip relationships were well-calibrated, thus lending confidence to the calibration procedure. However, it should be noted that the division in yielded and unyielded zones that was made during the analyses was based on the experimental results. In other words, the yield penetration was assumed a priori. Such data may not be generally (routinely) available. Nevertheless, it is clear that both the loss of bond at yielding and the asymmetry of yield penetration were very important in correctly  describing the structural behaviour. The asymmetrical yield penetration caused one of the sides to have a shorter available anchorage length and was thus crucial in the anchorage failure. It is common to use the symmetry of the test set-up and only model half the beam; this would not have worked when describing the asymmetric yield penetration for the studied beams. This therefore warrants further attention in future studies.

Results and discussion
The parameters in the bond-slip relationships were calibrated according to the procedure described in Section 4, for eighteen bars in nine beams. This section firstly presents and discusses the results of the FE analyses of the beams. The calibrated bond-slip parameters are then presented and discussed.

Overall numerical results
Exemplary comparisons between FE analyses with calibrated bond-   slip relationships and experimental results are presented in Fig. 17. Agreement is judged on the basis of load versus mid-span deflection responses and on the yielding penetration along the length of the bar at failure.
The FE analyses results for the nine beams were compared to the experimental results in Table 4 and Fig. 18. Mid-span deflection at endslip, load capacity at end-slip and yield penetration and asymmetry were compared. Overall, the FE analyses agreed well with the test results, indicating that the final calibrated bond-slip curves were reasonable. However, some disagreement may be observed. The load capacity at end-slip of bar 16H TI was slightly too large (19%) but was too low for bar 10I TO (− 20%). However, the load capacity at end-slip depended, not only on the calibrated bond-slip parameters, but also largely on the concrete compressive strength. In the analyses, the compressive strength of the concrete was assumed to be the same for all beams. Nevertheless, in reality, it varied depending on the deteriorated state and the observed high scatter in concrete properties [11]. As the other important parameters (mid-span deflection at end-slip, yield penetration and asymmetry) agreed well, the calibration of bond-slip parameters was deemed acceptable also for the bars in these beams. The maximum bond stress for unyielded zones (τ max ) was defined as the average bond strength calculated in the experimental study [11]. Thus, similar findings to those reported there may also be observed here;

Bond-slip relationships in FE analyses
top-cast bars generally had a lower τ max than bottom-cast bars when uncorroded but gained bond stress when corroded. τ max increased from 9.7 MPa to 13.4 MPa in top-cast bars, as the corrosion level increased from 0.57% to 4.7%. This may likely be attributed to the lower density of the concrete surrounding top-cast bars, which allowed for greater amounts of corrosion products to expand without spalling the surrounding concrete. Bottom-cast bars, on the other hand, lost bond stress when corroded. There was a 48.2% reduction in τ max when the corrosion level increased from 0.73% to 3.71% for bottom-cast bars. This was also confirmed in a later investigation into local bond-slip, which was also based on specimens extracted from Gullspång Bridge [16]. Moreover, the τ max values in the bars which failed in anchorage before yielding were significantly lower than those which reached yielding.
When compared with previous bond-slip constitutive models, one main difference can be observed: the maximum bond stress is considerably higher than in the constitutive models in Fig. 1, which yield a value between 1 and 2 MPa for the reinforcing bars in Gullspång bridge. This difference may be due to the models intending to provide values on the safe side. The calibrated slip at peak agrees with the value of the model by Feldman and Bartlett [10], 0.01 mm, while this value is larger in the other models (0.1 mm in Model Code [34] and 0.23 mm in Verderame et al. [38]); it is noted that they are comparable to the slip at the end of the plateau in the calibrated curves. The calibration of the parameter s 2 was an important part of the calibration process and had large influence on the structural behaviour of the beam. This indicates that the use of fixed values for the slip in bond-slip constitutive models may not be appropriate, especially if the effect of corrosion is to be taken into account.
Further, in Fig. 19, it may be observed that a higher τ max was accompanied by a higher residual bond stress (τ f ) and a shorter plateau (s 2 -s 1 ). This is further investigated in Fig. 20, which shows τ f and s 2 versus τ max . As may be seen in Fig. 20a, the residual bond stress increased with the maximum bond stress, along an almost linear trend. Spalling of concrete cover decreased the residual bond stress, while splitting cracks did not have any significant effect. In Fig. 20b, it may be seen that the slip parameter s 2 decreased with the maximum bond stress. For the bars that slipped without yielding, the slip parameter s 2 was smaller than for the others. This behaviour was anticipated, as the low maximum bond stress and short plateau (small s 2 ) caused the early anchorage failure in bars that did not reach yielding.
Although there were four parameters in the bond-slip relationship for yielded zones, three of them were determined as constant (Section 4.1). The slip parameter s y2 appears in Fig. 21, plotted against parameter s 3 . Both these parameters describe the slip at the beginning of the residual bond stress. The slip parameter s y2 increased with s 3 . In other words, slip at residual bond stress was greater in the yielded zone when it was large in the unyielded zone. The slip parameters, s y2 and s 3 , were closely related to the mid-span deflection at end-slip. Ductile deformation of the beam led to a ductile behaviour in the bond-slip relationship of the bar. This trend was almost unaffected by the corrosion level, damage level or casting position. Since the test results of mid-span deflection at end-slip were highly scattered among the beams, so were the calibrated slip parameters, s y2 and s 3 .

Comparison of calibrated bond-slip curves with pull-out test results
It is worth to note that the described procedure to calibrate bond-slip relations from beam test results is time-consuming; and also the experiments on beams require rather large efforts. An alternative structural test method was developed and used in a previous paper by the authors [16], in the form of pull-out tests of plain bars with specimens taken from the same bridge. Details on the testing procedure can be found in the aforementioned paper. Pull-out tests are much quicker and easier both to perform and to evaluate bond-slip relations from; it is therefore of interest to compare the results from these two methods.
In Fig. 22, the average peak bond strength over a corrosion interval of 0.5% is presented, both for bottom and top-cast bars. Experimental results from pull-out tests and 3-point bending tests are compared. The   Fig. 17. Comparisons between FE analyses and experimental results for two of the analysed beams.  Fig. 22. In Fig. 23, the residual bond strength, as obtained from the calibration of the bond-slip curves in this work, is plotted against the maximum bond strength, re-presenting data from Fig. 20a. Additionally, the residual bond strength as obtained from the pull-out tests (defined as the bond stress recorded at 4.5 mm of slip) is presented against the recorded maximum bond strength. Only results from top-cast bars were chosen,  given the more significant amount of data obtained in the FE analyses. Two linear interpolations are presented for reference bars, one for pullout data and one from FE results. The choice of reference bars is again based on the amount of data points.
Pull-out tests and calibrated FE analysis results of the beams show good agreement, both in magnitude and in trend. This is particularly of interest, since residual bond strength is hard to define experimentally from 3-point bending tests, but easy to obtain with pull-out tests. Nevertheless, the use of FE analyses allowed us to define the residual bond strength in the beam tests, and, thus, to compare it with the pullout tests results. Fig. 23 confirms what already observed in Fig. 22, indicating pull-out tests results as a reasonable choice for characterising the bond of plain bars. The data acquired present, however, a significant scatter, and individual pull-out test results are not always on the safe side with respect to behaviour in beams. This needs to be taken into account if it is chosen to employ data acquired from pull-out tests to define the residual bond of plain bars. Additionally, it is important that such comparison takes place between specimens with similar characteristics, such as the absence of corrosion induced cracks in the concrete cover, as in this case, or similar estimated corrosion levels. The acquisition of additional data, on different existing structures, is however necessary for reaching a better understanding on the usability of pull-out tests for structural assessment.

Conclusions
This work aimed to gain insights into the effect of corrosion and yielding of plain reinforcement bars on the bond and structural behaviour of RC beams. A total of nine previously tested beams were analysed using 3D NLFEA. Each tensile reinforcement bar was divided into yielded and unyielded zones, with different 1D bond-slip relationships assigned. The bond-slip parameters which characterised these curves were carefully calibrated. From the results and the calibration of the bond curves, the following conclusions may be drawn: 1. The FE analyses with calibrated 1D bond-slip relationships might provide a consistent description of the behaviour of the damaged RC beams with plain bars. Load-deflection curve, crack pattern, loadcarrying mechanism and asymmetrical yielding development agreed well with the experimental results. Thus, the shape used for the bond-slip relationship was suitable and the chosen parameters were adequately calibrated to describe the interaction between concrete and corroded plain bars. 2. Nine bond-slip parameters were needed for each bar (Fig. 7). Four parameters (τ ymax , τ yf , s y1 and s 1 ) might be set as constant for all bars, while two other parameters (τ max and s 3 ) might be manually calculated from test results and the other parameters. Only three parameters, τ f , s y2 and s 2 , needed to be individually calibrated for each bar, by comparing FEA and experimental results. Further, τ f and s 2 were correlated to τ max , which simplified the calibration process.
3. For the bond-slip relationships for unyielded zones, a higher maximum bond stress was followed by a shorter plateau and higher residual bond stress. Spalling of concrete cover decreased the residual bond stress, while splitting cracks did not have any significant effect. 4. A main finding of this study was that, to be able to describe the structural behaviour adequately, the loss of bond at yielding and the asymmetry of the yield penetration were both crucial factors. Whilst the former was already reported by other researchers, the latter, to the authors' knowledge, had not been previously observed. The asymmetry of the yield penetration caused one of the sides to have a  shorter available anchorage length. This, therefore, governed the anchorage failure mode.
To conclude, this study presents a novel approach to the assessment of concrete structures with plain bars. The use of FE models to interpret structural tests on naturally corroded structures showed the importance of including the bond behaviour of plain bars for a correct assessment, and specifically, highlighted the effect of loss of bond at yielding and the asymmetric behaviour of the yield penetration. This is particularly of interest since it is common to utilize symmetries in a structure to reduce computational cost. This could lead to a significant overestimation of the anchorage length. In the presented analyses, the position of the yield penetration was assumed as an input parameter based on measurements from experiments. This is an important limitation, as it is an unknown parameter a priori. However, observations on the asymmetry in the distribution of the yield penetration are given and could be used to estimate a possible reduction of the unyielded length. Additionally, comparisons with pull-out test results suggest that they can provide useful information to calibrate the bond-slip curve of plain bars; however, this deserves more attention in future studies.

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.