Prediction of the IC debonding failure of FRP-strengthened RC beams based on the cohesive zone model

Intermediate crack (IC) debonding failure is one of the common bending failure forms of fiber-reinforced polymer (FRP)-strengthened reinforced concrete (RC) beams. In this paper, a new prediction model for IC debonding in FRP-strengthened RC beams is proposed based on fracture mechanics and cohesive zone model (CZM), which takes into account the coupling effect of many parameters and has the advantages of high precision and simple expression. The nonlinear behavior of FRP-strengthened RC beams and the influence of flexural cracks are reasonably considered in this model, whereas all existing analytical models based on the CZM neglect this effects. To verify the accuracy of this model, we established a database containing 248 test data from the existing literature. By comparing the differences between the predicted and experimental results, we analyzed the causes of the error and established a semiempirical model. To test the reliability of the model, it is evaluated using the database constructed in this paper together with four representative strength models. The results show that the semiempirical model has a high accuracy.


INTRODUCTION
In recent years, externally bonded fiber-reinforced polymer (FRP) has become a common method for strengthening reinforced concrete (RC) structures (Teng et al., 2002). The load acting on the structure transfers the effective stress to the FRP through the interface adhesive layer, thereby enhancing the structural load carrying capacity. The debonding of the FRP and concrete interface is a common form of bending failure of FRP-strengthened RC beams, which may occur at the end of the FRP plate or at the intermediate crack (IC) (Fu et al., 2017). In the past 20 years, scholars from all over the world have performed much research on the prediction model of the IC debonding failure of FRP-strengthened RC beams (JSCE, 2001;ACI, 2008;Lopez-Gonzalez et al., 2016;Hoque et al., 2017;Li and Wu, 2018), but there is not a generally accepted model (Elsanadedy et al., 2014;Lopez-Gonzalez et al., 2016).
At present, the most recognized model in engineering design, namely the strength model, predicts the ultimate bearing capacity of FRP-strengthened RC beams through axial force prediction of FRP when IC debonding failure occurs. This kind of model has been adopted by many standards (ACI 2008;CECS146, 2003;JSCE 2001) and is still being improved (Lopez-Gonzalez et al., 2016;Fu et al., 2018;Li and Wu 2018). However, the evaluation based on experimental results show that there are few models that can accurately predict the capacity of IC debonding (Said and Wu, 2008;Elsanadedy et al., 2014;Lopez-Gonzalez et al., 2016). The strength model is generally obtained by directly or simplifying the stress state at the crack and then calibrating the database. Few parameters are considered in such strength model. Taking the Said and Wu (2008) 's model as an example, only three parameters of concrete compressive strength, FRP elastic modulus and thickness are considered in the model, and related research (Elsanadedy et al., 2014;Leung et al., 2006;Obaidat et al., 2013;Zhang et al., 2016) has confirmed that other parameters (such as height of the RC beam, tensile strain of tensile steel bars, and elastic modulus of FRP-concrete interface adhesive layer) will also affect the strain of the FRP when debonding occurs, thus affecting the debonding capacity.
A class of numerical models, such as the global energy balance model (Achintha and Burgoyne 2009;Hoque et al., 2017) , the non-linear local deformation model (Aiello andOmbres 2004, Ombres 2010) and the finite element model (Lu et al., 2007;Li and Wu 2018) , can directly analyze the stress state of the cracks in a FRP-strengthened RC beam, which considers the nonlinearity of material constitutive relations and more influence parameters, thus can accurately reflect the nonlinear behavior of strengthened beams during the loading process. The nonlinear behavior of the strengthened beam during the loading process is analyzed by examining the force state of the strengthened beam under each load, and using the corresponding criteria to determine whether debonding has occurred. However, in the analysis process of this kind of model, it is often necessary to use specific software to assist with the analysis, which is complicated compared to the strength model (Hoque et al., 2017;Shukri et al., 2018).
The expression for the IC debonding prediction model of the FRP-strengthened linear elastic precracked plain concrete beam based on the cohesive zone model (CZM) proposed by Wang (Wang, 2006a;Wang, 2006b;Wang and Zhang, 2008), which can directly analyze the stress state of the strengthened beam, is very concise. In this type of model, the crack is replaced by a rotating spring that does not consider the geometric size, so that the effect of the crack on debonding can be analyzed on a two-dimensional plane. After determining the relationship between the spring rotation stiffness and the interface bond-slip relationship, the debonding between FRP and concrete can be analyzed by nonlinear fracture mechanics method, which takes into account the coupling effect of several parameters. Although such models have been continuously improved in the past 10 years (Bennegadi et al., 2016;Chen and Qiao, 2009;Hadjazi et al., 2012;Houachine et al., 2013;Hadjazi et al., 2016), these analytical model is always created to quantitatively analyze the relationship between interface shear stress and the debonding failure of the FRP-strengthened linear elastic pre-cracked plain concrete beams. In such studies, the flexural cracks are pre-set before loading begins, and the effect of steel reinforcement and the nonlinear behavior of the strengthened beam are not considered. Therefore, such model has limited practical value and is not suitable for predicting IC debonding failure of FRP-strengthened RC beams.
To overcome the shortcomings of the above-mentioned prediction models, a new concise prediction model of IC debonding failure of FRP-strengthened RC beams based on the CZM is proposed by analyzing the nonlinear behavior of strengthened beams and the influence of flexural cracks during loading. This work is arranged as follows: first, we analyze and discuss the simplification idea of nonlinear behavior and multi-crack effects of strengthened beams in Section 2. The interface shear stress distribution equation for critical debonding is established according to the general derivation idea of the analytical model based on CZM. Then, in Section 3, the theoretical analytical solution of debonding capacity is determined after the boundary conditions at the time of debonding failure according to the interface shear stress distribution equation. We compared and analyzed the experimental and predicted results of the IC debonding failure capacity of 248 beam specimens subjecting three-or four-point loads. A more accurate semi-empirical model was proposed, whose result together with the prediction results of the four current and highly recognized strength models were compared. Finaly, we analyzed how to apply this model in practical situations in Section 4.
1. Note an assumption used commonly in the literature, both the FRP and RC beams are regarded as Euler-Bernoulli beams (Faella et al., 2008;Narayanamurthy et al., 2012;Razaqpur et al., 2020); 2. The RC beam did not experience prior loading after being strengthened with FRP and before being tested statically to debonding failure; 3. IC debonding occurs at a major flexural crack (Yao et al., 2005), which can be replaced by a rotating spring (Rabinovitch and Frostig, 2001;Rabinovitch, 2008) (Figure 2(a), Figure 2(b)). Under such loading conditions, the major flexural crack is located directly below the loading point (Lu et al., 2005); 4. After the IC debonding begins, the debonding will extends to the FRP end nearest to the major flexural crack rapidly, resulting in debonding failure (Lu et al., 2005). Thus only part of the strengthened beam between the spring and the support nearest to the spring of the strengthened beam needs to be analyzed (Figure 2(b)); 5. The bending moment of the interface adhesive layer is not considered.

Analysis of the influence of multiple cracks
The experimental results show that in the process of bending loading, the number of flexural cracks considered in the model determines the distribution of shear stress on the interface between the FRP and concrete. At present, there are three analysis methods for considering the influence of cracks: (1) Considering the effect of multiple cracks (Lu et al. 2007;Pan et al., 2010;Li and Wu, 2018). From the distribution of flexural crack of RC beam, this modeling method has the highest goodness of fit with the experimental phenomenon. However, the calculation cost of these numerical models is too high, which is not conducive to the promotion of engineering design. The research by Li and Wu (2018) shows that the prediction accuracy of these models may even be lower than that of the traditional strengthened model.
(2) Take no account of the effects of any cracks (Narayanamurthy et al., 2012, Razaqpur et al., 2020. Compared with the model considering the influence of multiple cracks, this kind of numerical model (Razaqpur et al., 2020) can ensure the prediction accuracy of FRP laminate strain to a certain extent, and its calculation cost is greatly reduced at the same time. However, in the analytical model (Narayanamurthy et al., 2012), the prediction results show a large degree of dispersion, and one of the reasons is that the model overly weakens the effect of cracks.
(3) Considering the effect of a single major flexural crack (Wu and Niu, 2000;Wu and Niu, 2007;Shukri et al., 2018), the calculation cost of this method is generally between (1) and (2). Although the interface shear stress distribution equation obtained based on this assumption is somewhat different from the actual situation, Shukri et al., (2018) found that such a model is more accurate in predicting debonding failure than the model considering multiple cracks.
From the perspective of the authors, it is difficult to accurately measure the effect of multiple cracks on IC debonding. The crack spacing will directly affect the interface shear stress distribution, thus affecting the prediction results of debonding failure. However, the dispersion of prediction results of the current FRP-strengthened RC beam crack spacing model are usually relatively large (Ceroni and Pecce 2009), which will significantly increase the prediction results of IC debonding failure. In addition, the calculation of the model considering the influence of multiple flexural cracks is often complicated, and it is difficult to give an analytical solution that is easy to calculate. Without considering the influence of cracks, the distribution equation of interface shear stress is often greatly different from the experimental results. Therefore, this study will try to analyze the prediction method of IC debonding failure under the condition that only the major flexural crack is considered.

Simplification of the nonlinear behavior of strengthened beams
In the process of bending loading, the constitutive relation of concrete and steel reinforcement (Figure 3(a)) under tension and compression has obvious nonlinear characteristics, so the flexural and compressive stiffness of the normal section at a certain position change constantly during loading. The interface shear stress distribution equation is a function related to stiffness, which will directly affect the interface shear stress, thus affecting the debonding bearing capacity. The numerical model (Aiello and Ombres 2004;Achintha and Burgoyne 2009;Ombres 2010;Hoque et al., 2017;Razaqpur et al., 2020) can often simulate the variation law of the stiffness of the normal section during the loading process. Although this method is very rigorous in theory, the calculation cost is obviously too high.

Figure 4
Typical experimental and analytical moment vs. curvature curves A large number of test results have confirmed that the flexural stiffness and compressive stiffness of the mid-span section of FRP-strengthened RC beams commonly used in engineering will undergo obvious three-stage variation with the cracking of concrete in tensile zone and the yield of tensile steel reinforcement (Garden et al., 1998;Gao et al., 2004;Attari et al., 2012). Take the evolution law of flexural stiffness as an example, the analytical moment vs. curvature curve includes three segments with different slopes when the normal section stiffness in each stage is approximately regarded as a constant ( Figure 4) (Attari et al., 2012). Therefore, under such conditions, the distribution of interface shear stress throughout the loading process can be expressed by the three equilibrium equations corresponding to the three groups of stiffness.
Considering that IC debonding failure usually occurs after the yield of tensile steel reinforcement (Said and Wu 2008), once the stiffness of Stage III in Figure 4 is determined, the distribution equation of interface shear stress under this state can be derived to determine the debonding bearing capacity. In such a state, the contribution of tensile steel reinforcement and concrete in the tensile zone to flexural stiffness and compressive stiffness of the normal section can be ignored. At the location of the crack, the RC beam compressive stiffness and the FRP compressive stiffness and , the bending stiffnesses rc D and frp D (relative to each neutral axis) in stage III can be obtained as follows: where c x is the neutral axis depth of the fully cracked plated beam (Figure 3(b)). c E and frp E are the elastic moduli of the concrete and FRP, respectively. This method can approximately reflect the relationship between the load and the deformation of the strengthened beam after the yield of tensile steel reinforcement (Kabir et al., 2018). It is noteworthy that Narayanamurthy et al., (2012) also used similar methods to calculate the flexural stiffness and the compressive stiffness of RC beams when predicting the debonding of the FRP plate ends. However, Narayanamurthy did not consider the three-stage variation of the stiffness, and believed that the section stiffness remained at the stage II ( Figure 4) throughout the loading process, which was one of the reasons for the large degree of dispersion of the predicted results of the model.

Bond-slip law of FRP-concrete interface
The existing IC debonding prediction models of FRP-strengthened plain concrete beams based on the CZM all derive the interface shear stress distribution law by defining a nonlinear interface bond-slip relationship (Wang, 2006a;Wang, 2006b;Chen and Qiao, 2009;Hadjazi et al., 2012;Houachine et al., 2013;Hadjazi et al., 2016;Bennegadi et al., 2016); thus, the bond slip curve (τ δ − curve for short, whereτ is the Interfacial shear stress and δ is the relative size of the slip) plays an important role in the prediction of IC debonding. The constitutive relation of the FRP-concrete interface is described by the bilinear model ( Figure 5), which is widely used to define the interface behavior of FRP-strengthened RC beams due to it being convenient for use and accurate prediction of interface debonding (Liu et al., 2007;Faella et al., 2008;Shukri et al., 2018;Razaqpur et al., 2020).The constitutive equations for the slip law expressed by the following equations: Latin American Journal of Solids and Structures, 2020, 17(6), e296 6/29 where max τ is the shear strength of the interface (the corresponding bond slip reaches 0 δ ), max δ is the bond separation slip when the interfacial shear stress reduces to zero. The area surrounded by the bilinear represents the interface fracture energy f G , which can be calculated as: (7) where ct f is the concrete tensile strength; a G and a h are the shear modulus and thickness of the adhesive, respectively. Figure 6 shows the evolution of the interface shear stress along the X-axis (Figure 2(b)) when the interface constitutive relation between the FRP and the RC beam is defined by the bilinear model and only the influence of a single crack is considered (Wang, 2006a;Hadjazi et al., 2012). E represents the region where the relative slip size δ is smaller Latin American Journal of Solids and Structures, 2020, 17(6), e296 7/29 than 0 δ and the relation between τ and δ satisfies the first equation of Eq. (6); S represents the region where the relative slip amount δ is greater than 0 δ and smaller than max δ , and the relation between τ and δ satisfies the second equation of Eq. (6). During the loading process, when the external load is relatively small, the interface shear stress distribution will reach the linearly elastic stage at first (Figure 6(a)). This stage ends when

Evolution of interfacial shear stresses
With the further increase of the external load, the distribution of interface shear stress will enter the elastic-softening stage ( Figure 6(c)), and a softening region (S region in Figure 6(c)) whose length is a appears near the crack. This stage ends when , and the length of softening region reaches its maximum value u a at this point ( Figure 6(d)).
If the external load continues to increase, the interface debonding will develop rapidly and macroscopic failure will occur (Lu et al., 2007). In order to predict debonding failure, the interface shear stress distribution equation for the elastic-softening stage should be obtained first, and then the analytical solution of the debonding bearing capacity should be determined according to the boundary conditions at the end of the elastic-softening stage (Wang, 2006a;Wang, 2006b;Hadjazi et al., 2012;Bennegadi et al., 2016). IC debonding failure usually occurs after the yield of the tensile steel reinforcement (Said and Wu, 2008) (i.e. Stage III in Figure 4), during which the interface shear stress distribution at the major flexural crack will always remain at the elastic-softening stage. This is because that the width of concrete cracks is generally greater than 0.01mm, while 0 δ ranges from 0.0021 to 0.0066 mm in the database established in Section 3.1 below.
The database was established by collecting 248 samples from 48 references, which covers the parameter variation range of common FRP-strengthened RC beams in engineering. Considering that the crack width is approximately equal to the sum of the slip size on both sides of the crack (Lu et al., 2009), it is obvious that the slip size at the crack after the tensile steel reinforcement yields will be greater than 0 δ , and the interface shear stress will present the distribution as shown in Figure 6(c). Next we will describe how to establish the interface shear stress distribution equation in Stage III.

Interface shear stress equation in Stage III
In this stage, the axial force, the shear force and the bending moment of the RC beam, the FRP and the whole section of the strengthened beam show following relations under any external load in Stage III.: where "T " represent the FRP-strengthened RC beam; are respectively the increments of the axial force, the shear force and the bending moment of the section after the end of Stage II, and their values vary with the change of the external load.
The infinitesimal isolator is selected at the crack after the tensile steel bar yields (Figure 3(b)), the following equilibrium equations are established: where ( ) 0 T N x = . The bending moment equilibrium equation is arranged in the center of the concrete in the compression zone (Figure 3(b)), we have: The constitutive law for RC beam and FRP read: represent the increment of axial and vertical deformation after the end of Stage II, respectively. To simplify the analysis, the curvature of RC beams and the FRP can be usually considered as the same (Smith and Teng 2001;Hadjazi et al., 2012): Substituting Eq. (11), Eq. (13), Eq. (19) and Eq. (20) into Eq. (16), we have When considering the contribution of shear deformation of interface colloid to the relative slip between the FRP and concrete, a equation similar to Chen and Qiao (2009) can be used to define the interface longitudinal displacement compatibility condition： (1) Interface shear stress distribution in region E Substituting Eq. (22) . Differentiating both sides of Eq. (23) with respect to x , we have: Differentiating both sides of Eq. (25) with respect to x and considering equilibrium Eq. (13), Eq. (17) and Eq. (18) gives the governing equation of shear stress along the interface between FRP and concrete: The solution of Eq. (26) can be expressed as .
where (2) Interface shear stress distribution in region S Substituting Eq. (22)  . Differentiating both sides of Eq. (32) with respect to x , we have: The solution of Eq. (35) reads According to the boundary conditions where, sy ε is the yield strain of the tension steel reinforcement; c ε , ' The maximum softening region length u a is a function related to the rotational stiffness r K of the spring and the

Semiempirical model
To evaluate the prediction effect of the model, a test database containing 248 test data is established (as shown in Appendix A table A1). The samples meet the following conditions: (1) The FRP is directly bonded to the RC beam through adhesive layer, and the strengthened beam has no anchorage or is only anchored at one end; (2) All the RC beams are strengthened with constant-thickness carbon, glass or aramid FRP sheets; (3) Failure of the specimens was due to IC debonding.  The average of the predicted to experimental value ratio (AVG), average absolute error (AAE), standard deviation (SD) , coefficient of variation (COV) and the range of predicted to experimental value ratio (Range) were used to describe the statistical characteristics of the predicted results (Table 1).The relationship between the predicted value p P and the experimental value e P obtained by Eq. (38) is shown in Figure 7(a), the results show that the predicted result of the theoretical model is generally smaller than the experimental result and has high dispersion, which is further confirmed by the statistical parameters shown in Table 1.There are two main reasons for this error: (1). Compared with the analysis method that considers the effects of multiple cracks, ignoring the effects of multiple cracks when the external load on the strengthened beam is the same will cause a greater slippage at the main crack, so the predicted value of the debonding bearing capacity will be lower. (2) Since the attenuation of the elastic modulus of the concrete during loading is not considered, Eq. (2) and Eq. (4) often overestimate the compressive stiffness rc C and bending stiffness rc D of the RC beam.

Figure 8 Calibration process
In order to further improve the accuracy of the model, the theoretical model needs to be calibrated. We found that the traditional polynomial fitting based on the least square method used by Narayanamurthy et al. (2012) could not effectively reduce the degree of dispersion of the predicted results, so the method used by Said and Wu (2008) to calibrate each parameter one by one is adopted in this paper. As mentioned earlier, the theoretical model proposed in this paper simplifies the attenuation of the elastic modulus of the concrete and the effect of multiple cracks on the predicted results, compared to the experimental phenomena, so the relevant parameters in the model need to be calibrated. However, we found that the degree of dispersion of the predicted results was not closely related to the change of c E . Considering that u a is related to the number of flexural cracks and will significantly affect the debonding bearing capacity (Chen and Qiao, 2009), the results of parameter analysis show that multiplying the calibration coefficient ϕ on the basis of u a can improve the prediction accuracy and reduce the degree of dispersion of the model. As we can see in (47) The prediction results of Eq. (47) are summarized in Figure 7(b) and Appendix A.

Comparison of the prediction results between the semiempirical model and the existing strength model
We selected four representative and highly recognized strength models to compare with the accuracy of the prediction results of the SEM. These strength models are the Said and Wu model (2008) Through the method of section analysis (Lopez-Gonzalez et al., 2016) (where the material strength reduction coefficient is taken as 1 and the material constitutive relation is consistent with the above analysis), the four strength models are used to predict the debonding failure of the 248 samples selected in this paper. The predicted results are shown in Figure 9 and Appendix A, statistical parameters are shown in Table 1. Except the Elsanadedy's model is slightly conservative, the AVG values of the other strength models are all around 1, so the degree of dispersion and the range of the predicted results of the model will become an important criterion for evaluating the model. It can be clearly seen from Figure 9 that the Said model, the ACI model and the Li model have a wide range of predictions due to the great errors of some samples; although the dispersion degree of Elsanadedy model looks slightly higher, the predicted range of this model is the lowest among all the strength models, which is further supported by the statistical parameters shown in Table 1.
In general, the prediction results of this model are slightly better than those of the above-mentioned strength models. On the one hand, it can be seen from Table 1 that the statistical parameters characterizing the degree of dispersion and prediction range in this model, such as SD, COV and Range, are superior to those of the above-mentioned strength model. Even when some samples (such as the No. 6 sample in Appendix A) that lead to a great error in the prediction results of the strength model are excluded, the coefficient of variation (COV) in this model is still slightly lower than the Said model with the lowest COV in the above-mentioned strength model. On the other hand, Figure 10 shows the distribution range of the absolute values of the relative errors of the samples. Orange, green, and purple represent the percentages of the samples with the absolute values of their relative errors ranging from 0 to 0.1, 0.1 to 0.2, and greater than 0.2, respectively. Ideally, the orange part should have maximum height, while the purple part should have minimum height. As we can see from Figure 10, the accuracy of SEM is obviously higher than that of the traditional strength model.

RECOMMENDATIONS FOR ANALYSING UNIFORMLY LOADED BEAMS
In the practical situations, the FRP-strengthened RC beams are generally subject to the uniform load. We used the test data of Pan et al. (2009) andMazzotti andSavoia (2009), to analyze the prediction results of the four strength models under uniformly distributed loads. The prediction results are shown in Appendix Table A2 and Figure 11. It should be pointed out that the specimens in groups D2-P2-L2, D2-P4-L2 and D2-P8-L2 were identical, but were subject to four-point, six-point and eight-point bending loads respectively. Pan's test result showed that with increase of the load points, the load condition became closer to ideal uniformly distributed loads and the IC debonding capacity of the strengthened beam also became greater (Pan et al., 2009). However, the same was not true for the four strength models, and the predicted values of the ultimate bending moments of these models were independent of the loading forms (Figure 11(a)). Although the stress increment model can reflect the influence of load distribution degree on the ultimate load to a certain extent (Pan et al., 2009), the prediction results of such models tend to have relatively high dispersion (Lopez-Gonzalez et al., 2016). Therefore, the prediction effect of such models was not analyzed in this paper.
Considering that there are relatively few studies on IC debonding under uniformly distributed loads, we suggest that under such loads, the debonding bending moment value should be the same as that of the strengthened beam under three-point loads, which is relatively conservative in theory. As shown in Appendix B and Figure 11(b), the prediction results show that this method is reliable.

CONCLUSION
In this study, a new IC debonding prediction model for FRP-strengthened RC beam is proposed. On the basis of results of the analysis described in the paper, the following conclusions are reached: 1. For the first time, CZM and fracture mechanics for predicting IC debonding failure of the FRP-strengthened RC beam is proposed based on the theoretical derivations and the available experiments, which can be used in engineering design conveniently because of its conciseness of calculation. Some parameters missing in available models have been incorporated into the new model proposed in the paper.
2. The nonlinear behavior of strengthened beams and the influence of flexural cracks are the main reasons for the complicated calculation of some numerical models. In order to reduce the calculation cost and better predict debonding failure, the analytical solution of IC debonding capacity is derived on the basis of considering a single major flexural crack and simplifying the three-stage variation of the stiffness of the strengthened beam during flexural loading to a trilinear model. Although the prediction results of such theoretical models tend to be conservative and show a high degree of dispersion, the accuracy can be improved and the degree of dispersion of the prediction results can be reduced by calibrating the length of the maximum softening region u a . This simplification and calibration method provides a more concise way to predict the debonding failure of similar nonlinear bi-material beams.
3. Compared with existing strength models, this model has less dependence on the empirical formula. We evaluated the predicted results of this model together with four highly recognized strength models by establishing a large experimental database containing 248 samples. The results show that, the semi-empirical model proposed in this paper has the highest accuracy and the lowest dispersion.