Bond-slip analysis via a cohesive-zone model simulating damage , friction and interlocking

A recently proposed cohesive-zone model which effectively combines damage, friction and mechanical interlocking has been revisited and further validated by numerically simulating the pull-out test, from a concrete block, of a ribbed steel bar in the post-yield deformation range. The simulated response is in good agreement with experimental measurements of the bond slip characteristics in the post-yield range of deformed bars reported in the literature. This study highlights the main features of the model: with physically justified and relatively simple arguments, and within the sound framework of thermodynamics with internal variables, the model effectively separates the three main sources of energy dissipation, i.e. loss of adhesion, friction along flat interfaces and mechanical interlocking. This study provides further evidence that the proposed approach allows easier and physically clearer procedures for the determination of the model parameters of such three elementary mechanical behaviours, and makes possible their interpretation and measurement as separate material property, as a viable alternative to lumping these parameters into single values of the fracture energy. In particular, the proposed approach allows to consider a single value of the adhesion energy for modes I and II.


INTRODUCTION
ohesive-zone models (CZMs) are widely used to simulate initiation and propagation of cracks along structural interfaces.They represent an effective alternative approach to fracture-mechanics-based methods for a wide variety of problems at very different scales, such as crack growth in dams, mortar-joint failure in brick masonry, bond-slip response of reinforcing bars in concrete, debonding of adhesive joints, delamination or fibre-matrix debonding in composites, among many others.Many of such problems entail combination of de-cohesion and frictional sliding, C which is often accompanied by dilatancy, in turn associated with the interlocking effect created by the asperities of the fracture surface.Therefore, in order to formulate models that properly account for the underlying physics of the problem it is essential to capture the distinct types of dissipation due to fracture and friction, the influence of the geometry of the asperities on the interlocking effect.Several interface models accounting for damage-friction coupling have been proposed in literature, see e.g.Del Piero and Raous [1] and references therein.Some of them are based on nonassociative softening plasticity, as for the multidissipative interface model proposed by Cocchetti et al. [2] and the contributions given by Bolzon and Cocchetti [3] and by Červenka et al. [4] in the field of concrete dams analysis, and by Giambanco et al. [5].A different strategy was followed by Alfano and Sacco [6], Alfano et al. [7] and, more recently, Sacco and Toti [8], where interface damage and friction have been combined in a cohesive zone model based on a simplified micromechanical formulation.The main idea was to consider a representative area at a micromechanical scale, which is assumed to be additively decomposed into an undamaged and a fully damaged part; moreover, it is supposed that friction occurs only on the latter.The evolution of damage is assumed to depend on the elastic energy in the undamaged part while the frictional behaviour is governed by a Coulomb law.To simulate dilatancy and interlocking this approach was adopted by Serpieri and Alfano [9], within a multi-scale framework in which, at a small scale, the asperities of the interface are represented in the form of a periodic arrangement of distinct inclined planes, denominated Representative Interface Area (RIA).On each of these planes the interaction is governed by the formulation proposed in Alfano and Sacco (2006).The above formulation by Serpieri and Alfano was recently revisited by Serpieri et al. [10], where it was shown that use of a single damage variable, combined with the choice of having a threshold damage function only depending on the damage variable itself and an equivalent displacement norm, requires coincidence of fracture energies in modes I and II to preserve thermodynamic consistency.Furthermore, it was shown that the enhancement of the model to account for friction and interlocking, based on the formulation proposed by Serpieri and Alfano [9], results in retrieving the experimental evidence that the measured fracture energy in mode II is typically quite higher than in mode I.For mixedmode loading, with positive (opening) mode I, the proposed model predicts an increase of the total fracture energy with increasing mode II loading component.Both facts are well supported by good agreement between experimental and numerical results.In this paper the model formulated in Refs.[9,10] is revisited and applied to the study of the pull-put test of a ribbed steel bar from concrete specimen, in which the interface model is used to simulate the interaction between steel and concrete within a two-dimensional axisymmetric simulation.The 'bond-slip' interaction between steel and concrete has been widely investigated in the last decades, both theoretically and experimentally [11][12][13][14].Many studies have revealed the significant role played by complex mechanisms including the progressive development of primary transversal and secondary inclined cracks [15], the influence of stress triaxiality on the formation of inclined struts in the concrete surrounding the steel bar [16], the fact that with increasing confinement the failure mode changes from being more brittle (splitting) to more ductile (pull out) [17].All these studies also revealed that, for ribbed bars (the only ones nowadays used), interlocking is the predominant mechanism of stress transfer between concrete and steel with respect to adhesion and friction (the latter intended as the theoretical friction that would be obtained if the fracture surfaces were flat).All these complex mechanisms are influenced by aspects that are not limited to the point-wise interaction between concrete and steel at the interface but include other structural details.For this reason, all the models developed somehow depend on details such as the diameter of the bars, the type of confinement, the cyclic nature of loading, among many others, and typically include a number of corrective coefficients that need to be empirically determined and are of difficult, if not impossible, physical interpretation for the practicing engineer.Hence, the aim of this paper is to show that very accurate prediction of the bond-slip interaction between concrete and steel bars can be achieved through a physically well justified model, developed within a general and thermodynamically consistent framework in which each contribution to the overall energy dissipation on the interface is accounted for with a simplified, yet effective multiscale description of the microscale.To show the effectiveness of the formulation, the model is validated by simulating the pull-out test of a steel bar from a concrete specimen and comparing the numerical results with the experimental ones reported by Shima [18].These have been chosen as benchmark test data because the relatively large dimensions of the concrete specimen, as well as the insertion of a clay sleeve with very small adherence in an initial part at the loaded end of the steel bar, prevents the formation of splitting cracks and results in rather negligible concrete damage or plasticity, except in a thin region immediately adjacent to the steel bar and its ribs, whose progressive damage is indeed incorporated in the cohesive-zone model.Therefore, the nonlinear response up to failure during the test is predominantly the result of the interface interaction, which is the focus of this investigation and is simulated by the proposed cohesive-zone model.This makes calibration of the model easier and gives more confidence in the robustness of the validation.

Multiscale formulation
n this section the fundamental ideas behind the frictional cohesive-zone model proposed in Refs.[9,10] are reported and discussed.For all the details of the derivation of the model and its implementation within an incremental nonlinear finite-element formulation, the reader is referred to the above references.An interface  is considered between two parts 1  and 2  f a body  .In this paper we will focus on two-dimensional problems whereby  is represented by a line.The model considers two scales: a macroscale where the interface is assumed to be smooth (Fig. 1(a)) and a microscale where the asperities of the fracture surface are considered (Fig. 1(b)).The terms microscale and macroscale are here used with aim to distinguish the smaller from the larger scale.In the context of the finite-element method, which is considered here, this implies the use of smooth interface elements to model the macroscale, whereby the details of the asperities do not have to be captured by the spatial discretization.Instead, at each integration point of each interface element, the cohesive law, which relates the relative displacement vector s to the interface stress σ , is determined by resolving a microscale problem for a Representative Interface Area (RIA).The RIA is assumed to be large enough to capture the details of the asperities.If the actual microscale is periodic, the repeating unit of the interface defines the RIA.In case it is not periodic or not perfectly periodic, the RIA can be taken large enough and with a geometry rich enough so that the response of the RIA can be considered statistically representative of the response of an 'infinitesimally small' area at the macroscale.In turn, this implies the assumption of 'separation of scales', which is typically required in multiscale homogenization [19].While the approach could be theoretically implemented by representing the asperities through a curved profile, it is assumed that such profile is defined by a finite number  N of microplanes.A further assumption that is made in the model is that the two sides of the interface profile within the RIA are infinitely rigid.In other words, the deformation of the asperities, as well as their damage, wear and fracture, are not considered.This means that a single vector s characterizes the relative displacement between the two sides of the RIA.On the k th microplane a local reference system

Thermodynamical formulation of the model
On each microplane, the stress is obtained implementing the model developed by Alfano and Sacco [7].The main idea behind such model is that an infinitesimal area dA of microplane can be decomposed into an undamaged part, whose I area is   1 k D dA , and a damage part, whose area is k D dA , k D being a damage parameter ranging from zero, in the case of no damage, to unity in the case of complete damage (Fig. 2).Assuming a linear elastic relationship for the undamaged part and a Coulomb-type frictional-contact relationship on the damaged part, the following expression of the free energy is obtained: where  k s the ratio between the area of the k th microplane and the sum of all such areas while, denoting with x  the negative part of x i.e.
in (local) modes I and II and f k s denotes the frictional slip on the k th microplane.Note that all microplanes are characterized by the same stiffness.The interface stress is obtained as the derivative of the free energy with respect to the relative displacement s: The frictional slip on each microplane is obtained through a non-associate Coulomb-type relationship defined by: Where  k represents the slip-onset function and  denotes the friction coefficient.
The evolution of damage variable k D , for each microplane, is defined to recover the two classical bilinear cohesive-laws in pure modes I and II depicted in Fig. 3, for the relevant local stress-displacement relation: In such relationships 0I s and 0II s represent the values of the relative displacements corresponding to the onset of damage in modes I and II, respectively, while cI s and cII s represent the 'critical values' corresponding to the attainment of complete damage.Furthermore, 0  I and 0  II denote the strengths in modes I and II, respectively, and it is easy to  According to [6], the evolution of the damage variable k D is given by: where The above evolution law for the damage variable k D could allow to consider different values of the fracture energies [11]  in modes I and II, cI G and cII G .However, it is shown in [10] that thermodynamically consistency is more 'neatly' achieved if these two values are taken to be equal.In fact, there is an ongoing debate in the scientific community concerning the underlying physical justification for a mode-mixity dependent fracture energy, that for zero or positive mode-I opening typically increases from the lowest value in pure mode I to the highest value in pure mode II, to the point that many authors have questioned whether the mode-II fracture energy can be considered as a true material property at all [20].In this aspect lies one of the primary strengths of the proposed approach, because the different contributions to the energy dissipation, i.e. to the measured fracture energy, are accounted for separately in the model and with a simple, yet physically sound, approach.In particular, the fracture energies cI G and cII G only represent here the energy dissipation due to the rupture of bonds; regardless of any theoretical reason on how thermodynamic consistency can be ensured, physical arguments suggest that such energy should be equal.

Experimental set up and finite element model
he model briefly described in the previous section has been implemented in a finite-step algorithm, providing the solution of the interface law at the integration points of interface elements, and coded in a user-subroutine (UMAT) for the finite-element code ABAQUS.To validate the efficiency of the model in capturing the concretesteel interface interaction against experimental results published in the literature, the pull-out test of s SD30 steel bar from a cylindrical block of concrete, studied by Shima et al. [18] has been numerically simulated.The experimental apparatus devised by Shima et al. is reported in Fig. 4(a).The diameter D of the steel bar is 19.5mm, and the set-up and specimen geometry are such that an initial 10 195  D mm long unbonded region in the vicinity of the loaded end of the bar is created.These features have been purposefully designed by Shima and coworkers in order to move the stress concentration well inside the concrete block and therefore avoid concrete damage and splitting at the top end of the concrete block.To this end, an initial 195mm long clay sleeve surrounding the bar in this region was inserted.

T
The load has been applied by prescribing the displacement of the top end of the bar up to a maximum value of 5.2 mm, with an applied strain rate of 10 -4 min -1 , the strain being measured at the top end of the bar and reaching a maximum value of 0.032.The strain along the bar has been measured using pairs of 5mm strain gauges, each pair placed on the bar diametrically opposite to each other, and at a distance of 49mm from the next pairs along the bar.A two-dimensional axisymmetric finite-element model has been used for the simulation employing a structured mesh.Both the steel bar and the concrete block have been discretized with 4-node fully integrated axisymmetric elements (named CAX4 in ABAQUS), with a size of 5mm in the concrete and 2.5mm in the steel bar.On the interface, 4-node axisymmetric interface elements (named COHAX4 in ABAQUS) have been used.Details of the geometry, mesh and constraints employed in the FE model are reported in Fig. 4(b).On the region unbonded by the clay sleeve no interface elements have been inserted on account of the prevented adhesion, and also because no data on the clay material properties were provided.The overall length of the bar and of the concrete cylinder is equal to 975mm.The concrete and steel material properties used in the simulation are reported in Tab. 1.For the concrete, a linear elastic material model was used due to the negligible damage and plasticity evidenced in the experimental results.Since the average cylindrical strength f c reported in [18] Figure 4: Pull-out test: (a) experimental set up [18] and (b) details of the geometry and the mesh of the finite element model.
The Poisson's ratio of concrete  c has been taken in accordance with values typically reported in the literature in the case of linear elastic behaviour.For the steel, the values of the Young modulus s E , the Poisson ratio  s , the first yield strength y f and the total strain  h corresponding to the onset of hardening after the initial plateau are those reported in [18].For the steel a J2 small-strain von Mises elasto-plastic model with nonlinear isotropic hardening has been used.The isotropic hardening curve in Tab. 2 relating the hardened yield strength  y to the equivalent plastic strain  p eq has been given in accordance with the uni-axial stress-strain curve reported in [18].

 
   deg.This value has been determined through a sensitivity analysis reported in Fig. 6.The properties of the interface have been determined through calibration and are reported in Tab. 3. The strengths in modes I and II have been taken equal to those determined by Serpieri and Alfano [9] for a similar problem, as they are known to have a limited effect on results.The same values of fracture energies have been taken for modes I and II in accordance with the observations made at the end of the previous section.Their value has been calibrated through a sensitivity analysis, whose results are provided in Fig. 7 and falls well within the range typically found for the adhesion between steel and concrete.The ductility parameter  has been taken very close to 1 because this results in values of the stiffness parameters  I I I K K high enough to well simulate the initial response when the interface is undamaged, yet not too high to avoid ill conditioning.In this range of values the sensitivity of the analysis to  is known to be negligible [22].The value selected for  has been taken as a mid-value in the range 0.7 and 1.4 recommended by the ACI code of practice [23] for the friction coefficient of concrete surfaces.

Numerical results and comparison with experimental data
Fig. 6 shows the interface tangential stress against the distance from the loaded end for varying microplane inclination , compared to the curve experimentally determined.The latter have been obtained in [16] based on the strains measured in the bars.It can be noted that with increasing  the peak value of the tangential stress in the vicinity of the loaded end increases, which is expected because of the increase in the mechanical interlocking with increasing microplane inclination.Fig. 7 shows the influence of the fracture energy and reveals that this parameter has a limited influence on the tangential stress profile, whose values tend to decrease at the loaded end of the bar and to increase at the free end with increasing c G .Fig. 8 reports the axial stress obtained in the FE simulation at the centre of the bar cross section against the distance from the loaded end, compared with the experimentally determined curve.In the FE model the axial force in the bar is completely withstood by the bar in the first 195mm, because of no adherence between steel and concrete in this part.Instead, the experimental results reveal that some adherence is actually present.This explains the initial, relatively small discrepancy in the curve, after which the agreement between numerical and experimental results is good.Similar considerations can be made for Fig. 9, where the numerical and experimental axial strain profile along the same path have been reported.Both curves report a steep change at the point separating the two parts of the bar in the plastic and elastic ranges.

CONCLUSIONS
he cohesive-zone model proposed in [9,10], which combines elastic damage, friction and dilatancy due to mechanical interlocking, has been revisited in this paper and further validated by simulating the pull-out test conducted by Shima et el.and reported in [16], obtaining good agreement between experimental and numerical results.The relatively large dimensions of the concrete block and the clay sleeve inserted at the top of the steel-concrete interface, which moves the stress concentration well inside the block, result in negligible damage and plasticity in the concrete block.This gives more significance and robustness to the presented validation because the only significant source of nonlinearity are the concrete-steel interface interaction and the plasticity in the mild-steel bar, with values of the plastic strains in the bar exceeding the end of the plateau of the stress-strain curve of the steel.The latter aspect is also significant because it has been observed that the deformation of steel bars affects the experimentally measured bond-slip law, particularly when reaching the plastic range [16].The main strength of the proposed model is that it effectively separates the total dissipated energy per unit of new crack surface, i.e., the measured mode-mixity dependent fracture energy, into the sum of three contributions.A first one is the rupture of bonds, which is mode-mixity independent and is therefore given by one single value of the fracture energy cI cII G G  nd, for each microplane, this translates in a single variable k D easuring the evolution of damage.A second contribution is given by the frictional dissipation that would be present even if, at the microscale, the interface presented no asperities, i.e. if it was completely flat.For the k th microplane this dissipative mechanism is related to the evolution of the inelastic slip s f k Finally, the third contribution, which has a predominant influence on the bond-slip relationship between concrete and ribbed steel bars, is given by the interlocking effect, which is here modelled through the above described multiscale approach and the definition of a pattern of microplanes in the RIA.

T
The remaining small discrepancies between the numerical and experimental results obtained in the presented validation can be partly attributed to the need of including some adherence between steel and concrete in correspondence of the clay sleeve.On the other hand, these can also be due to some aspects of the underlying physics that have not yet been included in the model, namely the finite size of the asperities in the RIA at the microscale, their deformation, wear, damage and failure and a more representative statistical distribution of their dimensions and geometry.These aspects will be investigated in the continuation of this research.

Figure 1 :
Figure 1: Multi-scale approach linking a smooth interface at the (a) macro-scale level to the (b) RIA characterized by a periodic arrangement of microplanes.

Figure 2 :
Figure 2: Decomposition of each infinitesimal area of microplane into damaged and undamaged parts.
, the areas under the bilinear curves in each mode represent the fracture energies in pure modes I and II and are denoted by cI G and cII G , respectively.

Figure 3 :s
Figure 3: Decomposition of each infinitesimal area of microplanes into damaged and undamaged parts.The assumption is then made that 0 is 19.6 MPa, the value of the Young's modulus c E in Tab. 1 was obtained using the following correlation suggested by the Italian code of practice [21]:

Figure 5 :
RIA geometry.The geometry of the RIA has been taken as shown in Fig.5.The three microplanes are all equal in length, resulting in a value them has zero inclination and the remaining two are symmetric and have an inclination angle equal to 1

Figure 6 :
Figure 6: interface tangential stress vs. the distance from the loaded end for varying microplane inclination 1 3     nd comparison with experimental results in[18].

Figure 7 :
Figure 7: interface tangential stress vs. the distance from the loaded end for varying fracture energy  cI cII G G nd comparison with experimental results in [18].

Figure 8 :
Figure 8: vertical direct stress in the bar vs the distance from the loaded end: comparison of the simulation results with the experimental results in [18].

Figure 9 :
Figure9: vertical direct (total) strain in the bar vs the distance from the loaded end: comparison of the simulation results with the experimental results in[18].
  t k representing the local normal and tangential directions, respectively.Hence, the relative displacement vector is decomposed on the k th microplane into its local mode-I component,

Table 1 :
Input parameters for concrete and steel.

Table 2 :
Hardening curve for the SD30 steel.

Table 3 :
Input parameters for the interface model.