Fracture-based interface model for NSM FRP systems in concrete

This paper introduces a numerical simulation tool using the Finite Element Method (FEM) for near-surface mounted (NSM) strengthening technique using fibre reinforced polymers (FRP) applied to concrete elements. In order to properly simulate the structural behaviour of NSM FRP systems there are three materials (concrete, FRP and the adhesive that binds them) and two interfaces (FRP/adhesive and adhesive/ concrete) that shall be considered. This work presents the major details of a discontinuous-based constitutive model which aims at simulating NSM FRP interfaces implemented in the FEMIX FEM software. This constitutive model was adapted from one available in the literature, originally employed for fracture simulation in meso-scale analyses of quasi-brittle materials, which is based on the classical Flow Theory of Plasticity combined with fracture mechanics concepts. The most important features of the implemented constitutive model are the consideration of both fracture modes I and II and the possibility of performing 2D and 3D analysis. In the end, results based on FEM simulations are presented with the aim of investigating the soundness and accuracy of the constitutive model to simulate NSM FRP systems’ interfaces. 2016 Elsevier Ltd. All rights reserved.


Introduction
One of the most effective techniques to strengthen and/or repair concrete structures consists on inserting a reinforcing material into a groove opened in the concrete cover of the element to be strengthened. This solution is known in literature as near-surface mounted (NSM) technique. Regarding the reinforcing material, fibre reinforced polymer (FRP) bars with rectangular, square or round cross-section have been widely used due to their several advantages when compared with steel [1]. In terms of the employed adhesives to bind FRP bars to concrete, epoxy adhesives are the most common ones.
Considering the local bond behaviour of concrete elements strengthened with NSM FRP systems, five failure modes can be found: two have cohesive nature and occur either (i) within the adhesive layer binding FRP to concrete or (ii) into the concrete surrounding the groove; other two failure modes have adhesive nature since they occur in the existing interfaces, namely, between (iii) FRP and adhesive or (iv) adhesive and concrete; finally, if none of the previous four had occurred, then the failure will happen by (v) FRP tensile rupture [2].
In the context of the present work, the bond behaviour of concrete elements strengthened with NSM FRP systems is discussed from the standpoint of numerical simulation within the Finite Element Method (FEM). Since the bond behaviour of such strengthening systems is normally studied by conducting bond tests, the focus of this work is more specifically devoted to the FEM simulation of NSM FRP bond tests [2].
In order to properly simulate this bond behaviour, three ''continuum" materials and two ''interfaces" need to be correctly simulated. These simulations include both physical representation and the material modelling of each one of them.
There are already available in literature accurate non-linear constitutive models aimed at simulating the post-elastic and failure behaviour of concrete and adhesive, e.g. [3][4][5]. The FRP can be simply assumed as a linear elastic. In terms of physical representation, line, surface or volume FE elements can be used to simulate all the three materials, depending on the type of FEM simulation to be performed. The simulation of those three materials is already well established, both in terms of elements and constitutive models.
On the contrary, in the case of the interfaces, there is a lack of proper constitutive models. Hence the following paragraphs present a review of the strategies that have been adopted to simulate the existing interfaces, focusing on the main target of this work, i.e. the simulation of NSM FRP bond tests. Isoparametric zerothickness interface elements have been classically used for this purpose. Depending on the type of elements adopted to simulate FRP, adhesive and concrete, the interface FE elements can be lines or surfaces, in order to assure compatibility between FE elements.
From a literature review on experimental programs of pullout tests, described in [2], some included simulation of the interfaces' behaviour. Essentially two types of strategies to simulate the interfaces were found: (i) the first strategy consist on simulating the bond behaviour with a set of closed form analytical expressions which were deduced from the physics of the observed phenomenon. Typically these mathematical expressions translate the different stages of stress transfer during a pullout test [6][7][8][9][10]; (ii) the second strategy consists on the use of advanced numerical tools, namely, the Finite Element Method (FEM) for simulating the interfaces. This later one can in turn be divided into two.
The first group consists in using closed-form analytical expressions as constitutive models of interface elements. The use of such strategy has proved to be very effective in terms of capturing the global behaviour of the entire system. Even though it has been widely used in the past, since the scope of this work is limited to bond tests, only four examples were identified. Three of them use 2D FEM simulations of beam [11] and direct [12,13] pullout tests. In all these three examples, only FRP and concrete were simulated with finite elements; the adhesive was simulated by the interface elements used in between FRP and concrete. Hence, the interface elements were used to simulate the joint behaviour of the adhesive and the two interfaces (FRP/adhesive and adhesive/ concrete). Finally, the fourth example consists on a 3D FEM simulation of beam pullout tests [14] in which the adhesive was simulated with volume finite elements.
The second group using FEM analyses, corresponds to approaches based on discontinuous constitutive models for zerothickness interfaces, which represents the main subject of the present work. Since the bond behaviour in NSM FRP systems has an inherent three-dimensional nature, only four works using this approach in 3D FEM analyses were found in the literature. Three examples consist of direct pullout tests, two with round [15,16] and one with rectangular FRP bars [17]. The fourth example consists of a beam pullout test with square FRP bars [18].
In [15], the adhesive/concrete interface was modelled with a frictional model based on a Coulomb yield surface. The interface FRP/adhesive was modelled using an elasto-plastic interface constitutive model originally developed for internal steel reinforcement. The yield surface of this model was defined by two functions. In tension, a Coulomb yield surface with zero cohesion and non-associated flow rule was adopted. In compression, another surface sets the limit in compression considering an associated flow rule.
In the second example of direct pullout tests with round FRP bars [16], the interface adhesive/concrete was not modelled, thus full bond was assumed between these two materials. The interface FRP/adhesive was modelled using a Mohr-Coulomb yield surface. This surface was limited by a normal stress equal to the tensile strength of the epoxy and by a limit value of tangential stress.
In the last example of direct pullout tests [17] the interfaces were not simulated since the experimental failure mode was not interfacial. Instead, a cohesive failure of the concrete surrounding the bonded length occurred. Once again, full bond between concrete and adhesive and adhesive and FRP strip was assumed.
Similarly, in the fourth example [18] the authors also considered full bond in both interfaces since the failure in their tests was cohesive within adhesive and/or concrete. Hence, no interface constitutive model was used.
Comparing the two strategies using FEM analyses presented above, in practical terms, there is essentially one main difference between them. While the first (using an analytical expression as constitutive model) is generally based on assuming a priori an analytical expression for the interface bond-slip law, the second strategy is completely conceived within the general framework of constitutive theories (e.g., fracture mechanics, plasticity, damage, among others) where the interface bond-slip law is not known a priori.
It is worth mentioning that the first strategy, based on analytical expressions, normally needs a lower number of parameters to be adjusted (depending on the analytical expression adopted), which may explain the higher use of such strategy when compared with the second one.

Interface constitutive model
Regarding the interface's constitutive model, it should have the ability of describing the two possible fracture modes in concrete elements strengthened with NSM FRP systems. Fig. 1 presents an example of a NSM FRP direct pullout test where the FRP can be seen moving simultaneously in x g 1 and x g 2 directions. The sliding movement in x g 1 direction is associated with fracture mode II while the opening movement in x g 2 direction is associated with mode I. As previously referred, most formulations used in literature are based on adopting ''a priori" analytical expressions for describing the interface bond-slip law. Furthermore, these formulations are generally based on assuming a fracture process in pure mode II and neglecting the effect of the interface normal stresses and the occurrence of out-of-plane displacements.
The constitutive model implemented in this work was provided with separate modules which allow performing 2D and 3D analyses considering either only mode II or both modes I and II of fracture simultaneously. This represents one of the key contributions of this work.
All the work presented in this paper was developed in the framework of FEMIX 4.0 [19], which is a freeware FEM software based on the displacements method. Fig. 2 presents the interface elements available in FEMIX on which the constitutive model was implemented. Particularly, it includes two line interface elements, with 4 and 6 nodes, which are schematically presented in Fig. 2a and b, respectively. Even though each of those interfaces can be used in 2D and 3D simulations, in this work only 2D line interface elements are addressed. FEMIX also includes two surface interface elements with 8 and 16 nodes ( Fig. 2c and d, respectively). Table 1 presents the three modules of the implemented constitutive model: the first module is used for 2D and 3D FEM analyses where only fracture mode II is considered (developed by [20]). Hence, the non-linear elasto-plastic behaviour is considered for local direction x l 1 , while the remaining directions behave elastically. a second module was developed for 2D FEM analysis where both fracture modes I and II are available (published in [21]). the third module addresses 3D FEM simulations where all the local directions have an elasto-plastic coupled behaviour (proposed by [22]).
The following section summarises the formulation of the three modules composing the implemented constitutive model. The most relevant expressions of all modules are presented together in the A. Further detailed information regarding each module should be found in [20][21][22].

Formulation
The constitutive model presented in this section is based on the classical Flow Theory of Plasticity. The basic assumption of this theory, in the context of small displacements, is the decomposition of the incremental joint relative displacement (designated as slip from this point onwards) vector, Ds, in an elastic reversible part, Ds e , and a plastic irreversible one, Ds p . The later is defined according to a general flow rule which depends on the plastic multiplier Dk and the plastic flow direction m. Hence, the relationship between slip and stress in the constitutive model is obtained by the following expressions, where Dr and D are the incremental stress vector and the constitutive matrix, respectively.
Dr e ¼ D e Ds e ¼ D e ðDs À Ds p Þ ð 2Þ Assuming that in a generic stress state n À 1, previously converged, the slips and stresses vectors and the hardening parameters are known, all these parameters need to be updated when a new increment of slip vector is added (step n). This update is performed by using the backward Euller method presented in the local return-mapping algorithm flowchart of Fig. 3. Block (2) of this flowchart corresponds to the beginning of the new step where the stress is updated by adding the new increment of slip. Then, if the new stress state lies inside the yield surface (i.e. the third residue f q 3;n is negative -see block (4) in Fig. 3), the actual stress state is in the elastic phase, otherwise it has a plastic component that must be accounted for. This is made using an iterative Newton-Raphson method which requires the estimation of the Jacobian matrix, J, in order to estimate the variations d of stress and state variables (i.e., hardening parameter j and plastic multiplier Dk) in the new iteration. The Jacobian matrix (Eq. (4)) is obtained by deriving the three functions used to estimate the residues necessary to check the stress state (see block (3) in Fig. 3), as shown in Eq. (4). This algorithm is then repeated until convergence is reached, i.e. until all three residues are lower than a predefined tolerance (see block (5) in Fig. 3).
J ¼ @f 1 @r @f 1 @j @f 1 @Dk @f 2 @r @f 2 @j @f 2 @Dk @f 3 @r @f 3 @j @f 3 @Dk 2 6 6 6 4 3 7 þ Dk @m @r Dk @m @j m À @Dj @r 1 À @Dj @j À @Dj @Dk n @f @j If the constitutive relation presented in Eq. (2) is true for elastic increments, it ceases to be when entering into the elasto-plastic regime. Hence, the elastic constitutive matrix shall be replaced by an elasto-plastic one. In this case, the expression of this new matrix can be deduced by imposing the consistency conditions and the Kuhn-Tucker condition presented in Eq. (5). Taking into account that the constitutive model was formulated under the work-hardening hypotheses, this condition can be rearranged to obtain the plastic multiplier (Eq. (6)), where the parameter H is defined according to Eq. (7). Replacing the plastic multiplier in the constitutive relation of the interface model, the elasto-plastic constitutive matrix can be obtained (Eq. (8)). Hence the new relation between slips and stresses is finally defined according Eq. (9).
Dk P 0; f ðr; jÞ 6 0; Dkf ðr; jÞ ¼ 0; Df ðr; jÞ ¼ 0 H ¼ À @f ðr; jÞ @k ð7Þ Appendix A includes all the expressions used in the formulation of the constitutive models. This includes yield function f, hardening variables U, yield surface gradient n, plastic potential g, plastic potential variables W, plastic flow direction m and hardening law Dj. In the following paragraphs few important comments are presented regarding those parameters.
In all the three constitutive models (CM) the hardening parameter is the plastic work, since, as referred before, work hardening was admitted in all formulations. However, the way the hardening parameter affects the yield surface is different in each CM since it depends on different variables. Hence, in CM II there is only a single hardening variable which is the shear strength, c, while in the other two CM, three variables exist: tensile (v) and shear (c) strengths and the friction angle (tan /).
The plastic potential surface of CM II and I/II_2D is not explicitly defined. However, since the formulation only requires the direction of the plastic flow, that is provided instead. The major difference between these two CM is that in CM II the plastic flow is associated while in CM I/II_2D a non-associated flow rule is admitted. Additionally, in CM I/II_2D an additional parameter exists which is the dilation stress, r dil . This stress corresponds to the normal stress at which the dilatancy vanishes when compression and shear stresses occur at the same time. The plastic potential of CM I/II_3D is described through a hyperbola identical to the yield surface but with different shear strength and friction angle (tensile strength is the same). This means that, in this model, plastic potential shear strength and friction angle need to be provided.
In terms of hardening law, it should be highlighted that, in CM I/ II_2D and CM I/II_3D, due to the different interaction that occurs between tangential and normal stresses, different expressions are used for the scenarios of tension and compression.
In all three CM, the evolution of the yield surface depends on the evolution of the hardening parameters, which depend on the evolution of the plastic work, W. The variation of the plastic work is considered by means of a dimensionless parameter (Eq. (10)), which translates the amount of fracture energy, G f , spent in a certain plastic work. Since CM I/II_2D and CM I/II_3D account for two fracture modes, there will be two dimensionless parameters in those CM, one for fracture mode I and other for fracture mode II.
Each dimensionless parameter is then input of a scaling function (Eq. (11)) in addition to a shape parameter (a) which can be different for each hardening variable. Eqs. (12)-(16) present the variation of each hardening variable where the indexes 0 and r, refer to the initial and residual value of the corresponding variable, respectively.
tan / ¼ tan / 0 À tan / 0 À tan / r ½ S n II ; a tan / À Á ð14Þ From the user standpoint, all the three constitutive models (CM) presented in the previous sections are included in a single global constitutive model and, depending on the type of analysis being performed, the user is allowed to set up one of the three.
To simultaneously exemplify and present the parameters adopted in the simulations further explained, Table 2 presents the required parameters in each CM.

Model validation: outline of test setups
The implemented constitutive models were validated using experimental results of direct pullout tests collected from the existing literature. In order to achieve a reliable validation, two examples were selected.
The first one, identified in this work as CaReCo, is fully described in [23] while the second, designated as GlRoTe, is documented in [24,25]. Figs. 4 and 5 show the geometry of the specimens and the test configurations of CaReCo and GlRoTe tests, respectively.
While they both consist of direct pullout tests, there are interesting differences between them, which justifies the simulation of both examples: 1. CaReCo tests uses carbon FRP (CFRP) while GlRoTe uses glass FRP (GFRP).

CaReCo specimens have rectangular FRP bars while in the case
of GlRoTe round FRP bars are used.
3. CaReCo and GlRoTe adopt test configurations which induce compression and tension, respectively, in the concrete specimens used.  4. The results of CaReCo include the post-peak slip versus pullout force curve (full range response) while the results of GlRoTe test are only up to peak pullout force. 5. In GlRoTe test strain gauges were used on the external surface of the GFRP along the bond length and their readings provided while in CaReCo test such data are not available since no strain gauges were used.

Details about the experiments and simulations
Each CaReCo specimen consisted of a plain concrete cube of 200 mm edge. On the side of the specimen a groove was made and a CFRP laminate was there inserted and fixed with an epoxy adhesive. The groove was 15 mm deep and 5 mm wide. The CFRP laminate with a rectangular cross-section with 1.4 mm thickness and 10 mm width was placed at the centre of the groove. GlRoTe specimens are prismatic plain concrete blocks (160 Â 200 Â 400 mm 3 ) to which a GFRP round bar with 8 mm of diameter was glued with an epoxy adhesive in a square groove with 14 mm cut on the concrete block.
To avoid premature failure of the specimen due to concrete cone formation near the top of the block, the anchorage length was initiated at 100 mm and 50 mm from the top, in CaReCo and GlRoTe specimens, respectively. The bond between the FRP and the concrete (L b ) was extended 60 and 300 mm downwards, in CaReCo and GlRoTe specimens, respectively.
On top of CaReCo specimen a steel plate with 20 mm of thickness was applied. To ensure negligible vertical displacement during the test this plate was fixed to the base by means of four M10 threaded steel rods. A torque of 30 N m was applied to each rod, inducing an initial state of compression in the concrete of about 2.0 MPa.
GlRoTe specimen was fixed to the base through two M20 threaded steel rods casted in the middle of the concrete block. Thus, in this test, both concrete and GFRP were in tension. Both types of test were monitored with a displacement transducer (LVDT) and a load cell. The LVDT recorded the relative displacement at the loaded-end between the FRP and the concrete (slip), while the applied force F was recorded through the load cell. Additionally, in GlRoTe test, five strain gauges were glued along the GFRP bar to measure its axial strains.
Based on the material characterisation conducted by the authors, a modulus of elasticity of 28.4/18.6 GPa, 165/51 GPa and 7.15/10.7 GPa was obtained in CaReCo/GlRoTe tests for concrete, FRP and adhesive, respectively. Since in both types of experimental tests the specimens failed by debonding at FRP/adhesive interface, all the non-linearity of the system was located at that interface. Hence, all materials were assumed linear elastic using the properties referred above (a modulus of elasticity of 200 GPa was used for the steel elements). Additionally, the interface elements were only used at the interface between FRP and adhesive, thus assuming that all the other regions of contact between different materials were fully bonded.
In order to assess the performance of the implemented interface constitutive model, two different FEM models were built for each type of test. Particularly, they differ essentially in the interface elements adopted, which were line 2D (L2D) and surface (S) interface elements.
Each FEM model was then run using either CM II or CM I/II, which resulted in four different FEM analyses for both CaReCo and GlRoTe tests. In the following paragraphs each single FEM model is described in detail. The parameters adopted in each constitutive model are presented in Table 2.

FEM model with L2D elements
In L2D FEM model, the direct pullout tests were modelled as a plane stress problem using the meshes represented in Figs. 4c and 5c for CaReCo and GlRoTe specimens, respectively. For both specimens the type of elements used was the same, namely: 4node Serendipity plane stress elements with 2 Â 2 Gauss-Legendre integration scheme for both concrete block and steel plate; 2-node frame 2D elements for both FRP bars and steel rods; 4-node interface L2D elements (see Fig. 2a) with 2 Â 1 Gauss-Lobatto integration scheme.
Both types of specimens were fixed to the corresponding testing machine by means of steel threaded rods. The only difference between them is that, while in the case of GlRoTe the rods were directly in contact with the concrete block, in CaReCo they were connected to a steel plate which in turn was in contact with the concrete block. Thus, the FEM support conditions in both types of tests consisted in fixing the bottom node of the steel rods. Additionally, unilateral contact supports were applied at the concrete block's base. Those restrain the downward movement in z direction (see Figs. 4 and 5), but allow upward free movement. In CaR-eCo test, the effect of the pre-stress in the steel rods was simulated by applying a uniform temperature variation to the rod elements equivalent to the torque applied.
In both tests the load was applied by means of a vertical prescribed displacement (direction z -see Figs. 4 and 5) in the top node of the FRP element.

FEM model with S elements
The S FEM model outlined in this section was built in order to test the S elements, thus deals with a 3D analysis with solid elements (see Fig. 4d). Due to computational costs, in each case only half of the specimen was modelled since both specimens have a symmetry on the xz plan.
In both CaReCo and GlRoTe concrete block specimens, steel plate (in the case of CaReCo), adhesive and FRP were modelled using 8-node solid elements with 2 Â 2 Â 2 Gauss-Legendre integration scheme. For the steel rods 2-node frame 3D elements were adopted. The interface elements were modelled with 8-node interface S elements (see Figs. 4 and 2 Â 2 Gauss-Lobatto integration scheme.
The test boundary conditions were simulated in a way similar to that explained in the previous section. Additionally, in these 3D simulations the displacements following y axis were also restrained along the symmetry plan.
The load was also applied by means of a vertical prescribed displacement, in this case, in all the top nodes of the FRP.

Parameters of each interface constitutive model
While in the analyses with CM II both L2D and S FEM models used the same input parameters for the interface constitutive model, in CM I/II simulations the parameters used by each FEM model (2D and 3D) were slightly different as shown in Table 2.
These differences in the parameters used in each simulation with CM I/II are related to the influence that the behaviour in the normal direction has in the global response. In fact, the behaviour in the normal direction of FEM model with S elements is affected by the stiffness of the surrounding materials (adhesive and concrete) which can be seen as a ''confinement" effect in the normal direction. Such influence does not exist when CM II is used since the behaviour in the normal direction is considered elastic in both L2D and S FEM models.

Model validation: numerical results
As previously referred, only CaReCo test results include the post-peak response while only GlRoTe test results provide FRP strains. Hence, for the sake of brevity, in the following sections the obtained results are presented and discussed only for CaReCo test. The only exceptions to this are related with the global response in terms of pullout force versus slip and the obtained FRP strains. The former is discussed for both (CaReCo and GlRoTe) in order to show the success of the FEM simulations. The later is presented and analysed in Section 4.4 for GlRoTe test only. Nevertheless it is worth to highlight that the trends and conclusions drawn in the following sections were very similar in the FEM simulation of both types of tests, thus are valid for both. Fig. 6a and b present the results of all the eight FEM analyses conducted, in terms of the relationship between pullout force and slip at the loaded-end. Each graph includes the experimental results envelope and the results for the FEM analyses with L2D and S interface elements, as well as using both CM II and CM I/II.

Experimental versus numerical results
For both types of test, in the case of L2D FEM models, the pullout force was taken from the top node of the FRP element (the node with imposed displacement) while the slip was taken from the top integration point of the top line interface element (see Figs. 4c and 5c). In the case of S FEM models, the pullout force was computed as the sum of all the forces obtained in all top nodes of the FRP element (loaded nodes), while the slip was obtained from one of the top nodes of the top surface elements (see Figs. 4d and 5d).
All the FEM analyses of CaReCo test successfully captured the three major stages of the experimental tests: the initial stage, governed by the chemical bond between concrete, adhesive and CFRP. Typically this stage is characterised by an almost linear behaviour; the second stage, corresponding to the system's stiffness degradation that occurs as a consequence of the progressive loss of chemical bond; the third stage (post-peak branch), governed by the friction that exist between the CFRP laminate and the surrounding adhesive.
The most remarkable aspect is related with the FEM model using interface S elements and considering both fracture modes (CM I/II). This is the FEM model which better captured the abrupt force decrease at the beginning of the post-peak branch. This, once again, highlights the three dimensional nature of the NSM FRP technique and the need for conducting 3D FEM analyses.
In GlRoTe, since the post-peak response was not registered, only the first two stages mentioned above were obtained. The FEM results were found to be very accurate in the first stage (up to a load level of 15-20 kN) as well as in terms of maximum pullout force prediction. Contrarily, the results in the middle region of the pullout force versus slip curve were not as accurate. The authors believe that this inaccuracy should be associated with acquisition difficulties during the experimental tests.
In addition, for GlRoTe test the beginning of the post-peak FEM curves is also included. This suggests a sudden pullout force decrease which can also justify the difficulty in capturing the post-peak response experimentally.

CM II versus CM I/II results
Figs. 7 and 8 present the graphs with the evolution of interface's slips and stresses in the simulations with both L2D and S FEM models using CM II and CM I/II, respectively. In all curves the horizontal axis corresponds to the 60 mm bonded length. For the sake of readability, the graphs only include two curves in the pre-peak phase for load levels of 5 and 15 kN, the curve for the peak load (F fmax ) and two curves in the two post-peak phase for load levels of 15 and 5 kN.
In the FEM models with L2D interfaces, slips and stresses were monitored at the integration points of the L2D interface elements which coordinates coincide with those of the interfaces' nodes.
In order to get, for each parameter, a curve comparable to that obtained in the models with L2D interfaces, in FEM models with surface elements, slips and stresses were read at the middle integration points of the two middle columns of surface interface elements. The referred reading points are inscribed inside circles in Fig. 4d.
Considering that there are differences in terms of numerical integration between 2D and 3D FEM models, the first conclusion that can be taken is that the results when using CM II are practically the same for both L2D and S FEM models, while they present some differences when using CM I/II. In fact, while in the FEM models using CM II the curves for L2D and S seem to be just slightly shifted (as a consequence of the referred differences in the numerical integration), in the FEM models using CM I/II they are actually different in terms of shape. This corroborates the previously referred influence of the normal direction behaviour.
Another conclusion that can be drawn is related with the curves of the parameters in the normal direction. Those are only presented for the FEM models using CM I/II and show that, when using L2D interface elements there is normal slip while when using S interface elements the normal slip is almost zero. As a consequence, the opposite occurs in terms of normal stress, i.e. when using L2D the normal stress is almost zero while when using S compressions are obtained in the normal direction. These findings corroborate the ''confinement" effect that only exist in 3D simulations due to the influence of the surrounding materials stiffness, as previously mentioned (see Section 3.4).

L2D versus S FEM models results
As previously mentioned, the bond phenomenon in the context of NSM FRP systems is intrinsically a three dimensional problem, even though it has been shown that such problems can be successfully simulated using 2D FEM analyses. However, there are some important aspects, like the ''confinement" effect shown before, that can only be simulated using 3D analyses. In addition to that, the type of information that can be obtained from 3D analyses is richer than that obtained in 2D analyses. As an example, Fig. 9 presents the contour plots along the S interfaces for both 3D FEM models for the peak pullout force. This figure includes all the three components of slip and stress in the three local directions of the interface elements.
As a reference it should be said that the graphs of S interfaces presented earlier in Figs. 7 and 8 correspond to the slips/stresses along the middle vertical line in each plot of Fig. 9, which coincides with the location of the CFRP and L2D interfaces in the 2D models. Now, a global picture of what happens along the entire perimeter of the interface between CFRP and adhesive can also be seen.
Firstly, this figure shows that the effect of the eccentric location of the CFRP laminate is well captured by means of the Fig. 7. Results of CaReCo FEM simulations using CM II and L2D or S interface elements: slip (a) and stress (b) along the interface in the loading direction.
interface-based modelling. For example, the slips in l 1 direction (first plot in both Figs. 9b and c) are slightly larger in the left side than in the right side, which correspond to inner and outer sides of the bond, respectively. This effect should be associated with the downwards movement of the concrete block as the CFRP is being pulled upwards, which should be smaller closer to the concrete block outer face.
Secondly it shows the different behaviour obtained when using CM I and I/II. This is more evident in the stresses along l 1 direction (fourth plot in both Figs. 9b and c). This plots show that due to the elastic behaviour in the remaining direction when using CM II, the l 1 stresses are similar in all the three sides of the groove that were simulated. Contrarily, when using CM I/II the behaviour in all the three sides is quite different. In fact, comparing the region closer to the loaded end of each groove side, values of approximately 14, 16 and 18 MPa of tangential stress can be found, at maximum pullout force, for the outer, middle and inner sides of the groove, respectively. This numerical observation can be explained by the curvature that occurs at the CFRP during the pullout. This further highlights again the different loading stage of each region of the interface. In addition, with these plots, it is easy to identify the regions where the interface remains in the elastic range and those where it already entered the softening stage.

Experimental versus numerical FRP axial strain
In the FEM models with L2D interfaces, the GFRP strains were read at the integration points of the GFRP elements which do not coincide with their elements' nodes. This is related with the adopted integration schemes. In the FEM models with surface elements, the strains in the GFRP were obtained from the integration points closer to the centreline of the concrete block front face, in order to match the position of the strain gauges in the experimental tests (inscribed in circles in Fig. 5d). Fig. 10 presents the evolution of the axial strain in the GFRP obtained in the experimental tests and the corresponding FEM results for different load levels. The FEM results include the 2D simulations using L2D interface elements and the 3D simulations with surface interface elements. Even though FEM simulations using both CM II and CM I/II were carried out, only one result is presented in order to do not overcharge the graphs. However, it should be stressed that the results were very similar in both cases.
The strain gauges provide discrete readings in the regions of the GFRP bar where they were installed. Hence, all the curves presented in Fig. 10 include a symbol to sign the regions where the strain gauges were in the experimental tests. As it can be seen, up to a load level of 20 kN the results are quite satisfactory. After this load level there are important differences between the experimental and numerical curves. However, this is true either in terms of strain (shown in Fig. 10) or pullout force (shown in Fig. 6b). In fact, as already mentioned in Section 4.1, the experimental response after the load level of 20 kN is quite unusual. Since there is a direct relation between GFRP strain and pullout force, if the later is not well captured, the former will not be well captured as well. At the peak pullout force, again the numerical model captured very well the results obtained in the experimental test.

Conclusions
In this work, the major details about the implementation of an interface constitutive model (CM) in the FEMIX software were presented. It shall be emphasised that this CM is an adaptation of three already existing CM for quasi-brittle materials. One of the CM only allows for accounting fracture mode II while the other two CM deal with considering both fractures modes I and II in 2D or 3D FEM simulations.
Hence, the main contributions of this paper were, in the first place, bring those three CM to the field of NSM FRP systems' interfaces simulation. Secondly, implementing the three CM as a single CM in order to made available in FEMIX a single and complete interface model. The third contribution corresponds to the presentation of FEM simulations with the developed model, thus highlighting its validation.
The later contribution is specially important, since this work adds to the literature examples of 2D and 3D FEM simulations of pullout bond tests, either using only mode II of fracture or combining both modes I and II together. Additionally, it was shown that the implemented model can be used with line or surface interface  elements in the framework of the well-known discrete crack analysis.
Regarding the results of the performed simulations, a good agreement between the experimental results and the numerical ones was found in all simulations performed in terms of pullout force versus slip. In addition, it was shown that further and more detailed information can be obtained when using surface interface elements. Namely, the effect of the eccentric location of the FRP bar is well captured by means of the 3D FEM simulations performed. The use of surface interfaces with CM I/II also allowed to verify that the bond behaviour varies, not only in the FRP tangential direction (load direction) but also in its normal direction (perpendicular to the loading). In fact, the value of maximum tangential stress varies from the outer to the inner regions of the interface FRP/adhesive. Finally, a comparison was made in terms of FRP axial strains where good agreement between experimental and numerical values was also obtained.