Development of Finite Element Analysis for Intermediate Length Coupling Beams Considering Bond-Slip Interface

Finite element analysis is performed on four reinforced concrete coupling beams of intermediate length using 2-D plane stress elements, under monotonic load up to failure. The model is verified using the results from (Nabilah and Koh in KSCE J Civil Eng 21:2807–2813, 2017). The bond-slip interface for the longitudinal reinforcement is modeled in the finite element, as it is found that it better predicts the load-deformation behavior compared to perfect bond. The comparison between finite element analysis and the experiment found that the model is able to predict the overall behavior of the structure, especially the maximum load capacity. The maximum deformation and the shear deformation from the finite element analysis are found to be underestimated, due to the inability of the model to predict shear deformation accurately. Flexural deformation (due to flexure and slip) is found to be well predicted, as the bond-slip behavior is modeled in the analysis. Generally, the shear deformation and slip are found to be significant in the intermediate length coupling beam and should not be ignored in the analysis. Finally, the effective stiffness prediction using finite element analysis is found to be overestimated and should be determined instead using existing equations.


Introduction
The behavior of beams is largely influenced by its span to depth ratio (or aspect ratio, L/d). For a slender beam with L/d greater than 4, the plane section theory could effectively be used for the determination of its flexural behavior. For a deep beam with L/d less than 2, shear is transferred by a single diagonal strut across the beam and it tends to fail in diagonal tension when shear reinforcement is insufficient. The behavior of intermediate beam with L/d between 2 and 4 is by far the least studied (Wallace 2012). Hence, an investigation on the behavior of intermediate length beams is warranted to determine its capacity in terms of load and deformation, to ensure it is not overestimated in the design.
Most of the analytical work on the coupling beam involves the development of a strut-and-tie model for the strength calculation (Shuraim and El-Sayed 2016;Zhao et al. 2018;Nabilah et al. 2019). Finite element modeling, considered as micro-modeling technique, is difficult and require specialized knowledge to obtain a meaningful result. However, if properly conducted, finite element analysis could give insightful information on the behavior of the structural member (Zhao et al. 2004;Mohr et al. 2007;Bower and Rassati 2008;Brena et al. 2009;Shastri 2010). Zhao et al. (2004) developed a 2-D finite element analysis (FEA) by modeling the dowel action and concrete confinement. Smeared crack model is used, and the beam is loaded monotonically up to failure. It is found that the model could predict the crack behaviour with reasonable Page 2 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 accuracy, based on experimental work conducted. Shastri (2010) developed a finite element model using ABAQUS, using the damaged plasticity model. The coupling beams are based on experimental work by Galano and Vignoli (2000) with a beam aspect ratio of 1.5, and by Bristowe (2000) with an aspect ratio of 3.6. Both beams undergo cyclic load up to failure. It is generally found that the finite element model could match the experimental result up to yielding of the steel reinforcement, while the postyield behavior could be further improved. Numerous models (micro and macro modeling, strut and tie model) have been developed to determine the shear capacity and ductility of coupling beam undergoing significant lateral load. Most of the models developed could capture the flexural and shear mechanisms rather accurately. Another significant factor affecting the deformation capacity (and ultimately initial stiffness and ductility) of the coupling beam is the slip of reinforcing bars (Toprak et al. 2015;Ding et al. 2017). This parameter should also be explicitly modeled to obtain an accurate representation of the overall behavior of the member, however, often ignored in the analysis. Zhao et al. (2004), in their research, found that the inability to model the slip resulted in stiffer behavior after initial cracking as the loading progresses.
In seismic design, the effective stiffness (I e /I g ) is an important parameter to assess the response of structural components under the applied load. Based on the design codes, the value recommended is in the ranges of 0.3 to 0.5. However, these values are considered to be an overestimation based on experimental work on coupling beam. Several researchers have proposed equations to determine the effective stiffness (Paulay and Priestley 1992;Taranath 1998), however, its determination involves large uncertainties due to many factors affecting it (Vu et al. 2014). In addition, existing equations on the effective stiffness is based on coupling beams with an aspect ratio of less than 2.0 (deep beam), and very minimal study is dedicated to intermediate length beams (L/d between 2.0 and 4.0).
In this study, a new finite element model is developed for intermediate length coupling beam to study its overall behavior in terms of maximum shear load and deformation capacities, cracking and failure mode, effective stiffness as well as to investigate the components contributing to the overall deformation of the coupling beams. In the development of the finite element model, emphasis is given on the material constitutive models, bond-slip interactions and element formulation.

Research Significance
Currently, there is a limited study on the development of a suitable finite element model for the coupling beam, where most of the theoretical model is based on the strut-and-tie concept. Moreover, the existing study did not emphasize the need to model the bond-slip interface, where it is found that it contributes to significant deformation and cannot be ignored. This study also concentrates in intermediate-length beams, whereas most of the existing works were done for short and deep coupling beams.

Experimental Work
To assess the behaviour of intermediate length coupling beams, 4 beams were constructed with a cross-section of 330 mm by 180 mm with 15 mm cover, and L/d of 3.4 and 2.7. The beams have varying longitudinal and shear reinforcements; however, the mode of failure is ensured to be shear after yielding of reinforcement (ductile shear). To ensure proper bonding between rebar and surrounding concrete, adequate anchorage is provided through bar bending and development length as accordance to ACI 318 (1995). The beams were tested upright as shown in Fig. 1 and subjected to monotonic load until failure. The detailing of tested specimens is summarized in Table 1, and testing procedure and instrumentations can be obtained from Nabilah and Koh (2017).

Material Constitutive Model
As this structure is expected to undergo large nonlinear deformation, care must be taken in modeling its tensile behavior. In this research, the rotating smeared crack approach is used, which assumes that the crack orientation changes continuously. To describe the cracking behavior, the exponential tension softening curve is used, which is a function of the mode-I fracture energy, G f and the crack bandwidth, h. The fracture energy is approximated using Eq. (1) based on fracture mechanics, while the crack bandwidth is taken as the square root of the total area of the element for higher-order 2D element. The tensile strength, f t is taken as 0.9 of the split tensile strength and the compressive strength, f' c is in MPa. For compression, parabolic compression curve by Feenstra (1993) is used, where the compressive fracture energy, G c is given in Eq. (2). The confined compressive strength is as given by Mander et al. (1988).
The steel is modeled as bilinear, with a modulus of elasticity (E s ) up to yield, and ultimate stress at a maximum strain of 0.1. All material properties are obtained from Nabilah and Koh (2017). The constitutive models used in the finite element analysis are given in Fig. 2. (1)

Bond-Slip Interface
When a reinforced concrete structure is loaded, multiple cracks will form around the bars (bond cracking) and within the concrete to the surface. This behavior is termed tension stiffening, which is the combination of the tension softening effect in micro-level and the bond across the reinforcing bar. Depending on the modeling scale in finite element, a perfect bond between the rebar and surrounding concrete may be assumed, provided proper tension softening/stiffening model is included in the analysis (fib 2008). However, for a medium-scale model such as in this research, the bond-slip behavior

Table 1 Detailing of tested specimen. (Modified from Nabilah and Koh 2017).
Beam  Page 4 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 between the rebar and the concrete must be modeled to obtain a better prediction of the result in terms of the deformation profile. Bond behavior affects the structure in both serviceability (SLS) and ultimate limit states (ULS). In SLS, bond influences the width and spacing of transverse crack, tension stiffening and curvature, while at ULS, effects the rotation of the plastic hinge region (fib 2013). As coupling beams undergo large inelasticity and significant deformation, the bond-slip interface must be modeled to obtain an accurate estimation of shear-deformation behavior. A large slip occurs due to significant transverse and longitudinal cracks near the reinforcement, allowing it to move relative to the concrete. The bond between the rebar and surrounding concrete is modeled through an interface element with zero thickness. This interface element assumes the relationship between traction and relative displacement normal to the reinforcement direction (t n and s n respectively) as linear elastic (Fig. 3a), while in the shear direction (t t and s t in Fig. 3b) as nonlinear. Consequently, the bond stress-slip relationship can be derived for shear direction and are developed and calibrated in many literatures (Lundgren and Gylltoft 2000;Lowes et al. 2004).
For monotonic loading of ribbed reinforcement embedded in structural concrete, the non-linear bondslip model by fib Model Code for Concrete Structures (fib 2013) is found to be adequate and simple for modeling in finite element analysis. The bond stress-slip behavior is represented in Eq. (3) and is translated in Fig. 4.
Assuming the pull-out failure mode, τ b represent the bond stress, τ bmax is the maximum shear stress, taken as 2.5 f ′ c (units in MPa), and s is the slip of the reinforcing bar. The values s 1 , s 2 and s 3 are the values based on for s 2 ≤ s ≤ s 3 τ bf = 0.4τ bmax for s 3 < s experiment and semi-empirical expressions, namely 1 mm, 2 mm and clear distance between the ribs (taken as 6 mm based on measurement) respectively. These values are generally similar (and slightly lower) than other models as summarized by Murcia-Delso and Shing (2015). The bond stress-slip relationship is described in four distinct phases as given in Eq. (3). The ascending branch describes the bonding between the ribbed reinforcement within the concrete matrix, with the inclusion of local crushing and micro-cracking effects. The region of plateau exists for confined concrete, where the bond is sustained after the maximum bond strength is reached. The descending branch represents the decrease of the bond strength due to the shearing of the concrete between the ribs until the residual bond stress of 0.4 of the maximum shear stress is reached.

Finite Element Model
The structure is modeled in 2-D using plane stress elements, as the out-of-plane stresses could be neglected. Hence, the concrete elements are modeled using 8-noded quadrilateral isoparametric plane stress elements. The ratio of the longer to the shorter length of each element is ensured to be not more than 1.5. As

Cross-sectional view (a)
Longitudinal view (b) t n , s n t t , s t  Page 5 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 steel reinforcement does not provide significant bending stiffness, 2-noded linear 2-D truss elements are used to model the bars. The steel reinforcements (both longitudinal and transverse reinforcements) were modeled as a discrete element with a bond-slip interface within the concrete. The element sizes are approximately 13.5 mm. The coupling beam is modeled upright, similar to the testing configuration. However, the finite element model shown in Fig. 5 is horizontal, to reflect the actual nature of the beams. The supports are modeled as elastic element since sufficient reinforcements are provided in the test sample to ensure rigid supports. The bottom support (or left support in Fig. 5) is restrained in the X and Y directions, while the top support (or right support in Fig. 5) is restrained only in rotation to ensure the equal end moments can be achieved. The loading is applied at the surface of the rotational restrained (top support) as deformation controlled, and all elements are assigned with own self-weight. The finite element model for beam B2.7-2 is shown in Fig. 5, as an example. Figure 6 compares the result obtained from the experimental work with the result from the numerical analysis with and without considering the bond-slip behavior of the longitudinal reinforcement for beam B3.4-2. As shown, the load-deformation curve closely matches the experimental result when a proper bond-slip model is considered in the analysis.

Load-Deformation Analysis
When perfect bond (no slip) is assumed, the loaddeformation curve is found to be stiffer prior to yielding, and the sudden drop of load capacity is observed at two points prior to failure. In addition, the maximum deformation is found to be lower compared to the model with the proper bond-slip interface. At the peak of the load as marked by 'X' in Fig. 6, it is observed that the load carrying capacity reduces slightly after yielding of the reinforcement before continuing up to failure. This is because as the reinforcements are modeled to be perfectly bonded to concrete, the sudden increase in the axial strain of longitudinal reinforcement after yield causes an increase in local tensile strain in concrete. This resulted in the reduction of the load carrying capacity before the stresses are redistributed in the elements and stabilize the whole system. The sudden high tensile strain is not observed when the bond-slip is modeled, resulting in a more gradual stress redistribution hence eliminating the drop in the load capacity. This phenomenon is also observed in all other beam samples modeled in finite element analysis. Hence, it can be concluded that bond-slip has to be modeled in finite element analysis to obtain a more accurate result, and the model by fib Model Code for Concrete Structures (fib 2013) is found to be adequate for the coupling beams.
In subsequent analysis, the bond-slip behavior is explicitly modeled in the finite element analysis unless otherwise indicated.
The summary of maximum load (P max ) and maximum deformation (∆ max ) for both experimental work (Exp) and finite element analysis (FEM) are shown in Table 2. Overall, the maximum load capacity is found to be in good agreement with the experimental result, however, the maximum deformation is found to be underestimated for most of the beams. Figure 7 shows the load-deformation curves for all tested beams compared to the finite element analysis results. In all of the beams, the stiffness before cracking of concrete obtained from the finite element analysis are  Page 6 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 found to be almost similar to the experimental result. However, most of the models tend to be stiffer until the yielding of reinforcement, possibly due to the inability of the model to accurately determine the shear deformation. The finite element model is found to underestimate the maximum deformation, particularly in beams B3.4-3 and B2.7-3. These two beams have a higher amount of longitudinal reinforcement compared to its counterpart, hence having much larger yield and ultimate load capacity. A closer inspection of the ductility (ratio of the maximum to the yield deformations) obtained from the experiment and finite element model for each beam shows similar result. This is due to the slightly stiffer pre-yield load-displacement curve and smaller maximum deformation in the finite element analysis. Hence, it can be concluded that the finite element model could effectively estimate the ductility of the beams. It should also be noted that during testing of beam B3.4-3, excessive support movements and rotations are observed due to inadequate boundary condition. That is possibly also the reason for large underestimation of maximum deformation compared to the experiment for the beam. Generally, the finite element models are able to capture the overall behaviour of the structure and predict its point of failure with reasonable accuracy.  Page 7 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 Figure 7 also shows the components that contribute to the total deformation of the beams (refer to Fig. 7a for explanation). In this analysis, three major deformation components were identified, namely shear, flexure and slip. In general, shear contributes to very small deformation, and the value increases as the load increases beyond concrete cracking. Although small, the shear deformation is found to be quite significant, accounting for 20 to 27 percent of the total deformation before yielding, with a larger contribution in the shorter beam. Flexure contributes to large deformation even before concrete cracks, however, as deformation increases, the slip of the longitudinal bars become significant. This shows that in intermediate length coupling beams, shear deformation and slip is significant and cannot be neglected in determining the deformation capacity (in other words ductility) of the member.

Deformation Components
The deformation of the beam member is divided into its components, namely flexure and shear, and they are measured appropriately in the experimental work conducted. However, during the experiment, the deformation due to slip is not specifically measured (due to limitations in instrumentation), hence the flexural deformation obtained is inclusive of both flexure and slip of the longitudinal reinforcement. Figure 8 shows the components of shear and flexural deformations for each beam based on the experimental work and finite element analysis.
In general, the finite element analysis is able to predict the flexural deformation (flexure and slip) of the beams rather accurately, except for B3.4-2 (Fig. 8a) where the analysis overestimated the flexural deformation after cracking of concrete. For all of the beams, the deformation due to shear is found to be underestimated in the finite element analysis, compared to the experimental result. The differences in the values become profound after the concrete cracks and increases with the applied load. It is also noted that this error is more noticeable in beams B3.4-2 and B3.4-3 (Fig. 8b, c). A possible reason is that the shear demand (and subsequently shear deformation) is higher in the shorter beams and cannot be captured accurately in the model.

Crack Pattern and Failure Mode
The crack pattern and failure of the beams are observed from finite element analysis (Fig. 9) and compared with the experimental results. In general, the crack at failure in finite element analysis is very similar to the experimental work. For beam B3.4-2, large cracks are observed diagonally as well as at the beam-support interface.
The failure is shown to be due to flexure, however, the observed failure in the experiment is diagonal tension. For beams B3.4-3, B2.7-2 and B2.7-3, the largest cracking strain is observed at the interface between beam and support, indicating that the failure is due to flexure. This failure is also observed during the experiment. In beam B2.7-3, cracks extending to more than half of the beam length is observed before the failure of the structure. It is due to the higher amount of longitudinal and transverse reinforcements, creating confinement which improves the strength of concrete. This generally increases the concrete capacity, while allowing deformation to extend significantly throughout the length. In all beams, the yielding of the bars is concentrated near the beam-support interface. However, the increase of strain in longitudinal reinforcement extended into the significant height of the beam, causing cracking along its length, as shown in Fig. 9. In general, the extent of a diagonal crack in all beams is found to be similar to the experimental work.

Effective Stiffness
Coupling beams under seismic load will undergo significant inelasticity before failure. In design, the effective stiffness has to be properly determined as it will affect the overall behavior and fundamental period of the structure. Vu et al. (2014) proposes an equation for the determination of the effective stiffness (κ Eq ) and is shown in Eq. (4).
where ρ v and ρ l are the transverse and longitudinal reinforcement ratio respectively. Hence, the comparison of the effective stiffness obtained from the experiment, finite element analysis, and theoretical equation is shown in Table 3.
In general, both estimations based on finite element analysis (κ FEA ) and equation (κ Eq ) overestimated the effective stiffness except for beam B3.4-2. As shown in Table 3, Eq. (4) gives slight overestimation of the effective stiffness, which is in-line with the findings by Vu et al. (2014). However, this prediction is still more accurate than the value provided in the design codes. This is because the equation by Vu et al. (2014) has been calibrated using numerous coupling beam experiments, causing the values to be closely estimated with considerably low error. In general, the use of finite element will result in an overestimation in the effective stiffness and is not suitable for its determination. The possible reason is due to the overestimation in the flexural deformation of the beam, resulting to lower effective stiffness in the finite (4) κ Eq = I e I g × 100% = 0.67 1.8 Page 8 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 element analysis. However, the finite element analysis is still very useful to determine the maximum load and deformation, cracking pattern and deformation components as discussed.

Conclusions
Finite element model has been developed for intermediate length coupling beams with L/d of 3.4 and 2.7, based on the paper by Nabilah and Koh (2017). The beams are designed to fail in shear after yielding of reinforcement. The finite element model is developed using plane stress element and the bond-slip interface is modeled using fib Model Code for Concrete Structures (fib 2013). The following are the findings.
• Modeling the bond-slip is necessary to better simulate the overall behavior of the beams. When a perfect bond is assumed, it is observed that there is a sudden drop in the maximum capacity due to stress redistribution in the elements. This phenomenon is eliminated when the slip of reinforcement is modeled, with improvement in the prediction of stiffness. • The model developed is able to predict the behavior of the beams accurately, especially the maximum load capacity. The maximum deformation was found to be underestimated, due to the inability of the model to accurately determine the shear deformation. • The deformation of the beam was divided into two different components, namely shear and flexural deformations. The shear deformation was found to increase beyond cracking, becomes significant before yielding (20 to 27 percent of the total deformation), and is found to be larger for a shorter beam. The shear deformation was underestimated for all beams, signifying the inability of the finite element model to accurately estimate the shear deformation. The flexural deformation (including slip) in the finite element model was well deter-

Shear deformation
Flexure and slip def. Shear deformation Flexure and slip def. Page 9 of 10 Nabilah et al. Int J Concr Struct Mater (2020) 14:33 mined, where the slip becomes significant as the concrete cracks. This shows that the slip is very important in coupling beams and should be modeled. • The finite element model developed was found to be able to predict the failure mode, where most of the beams fail in flexure. • The effective stiffness from the experiment is compared to the one obtained from the finite element and an analytical equation. It is found that the equation estimated the effective stiffness better than the finite element result. Hence, the finite element might not be suitable for its determination, as it generally gives an over prediction of the effective stiffness.

Funding
The authors declare that no funding was received 2 ) c (