Implementation of interface damage model with friction to concreteFRP shear connector

Fibre-reinforced polymer (FRP) composites have become increasingly prominent across a range of industries, with significant impact in the fields of construction and bridge engineering. Their expanding utilisation in structural applications especially for the rehabilitation, strengthening, and retrofitting of existing bridge infrastructure represents a key area of advancement [1-3]. The integration of FRP materials enables the development of lightweight, highstrength, and corrosion-resistant structural systems, which are especially advantageous in the harsh environmental conditions of bridge construction and long-term maintenance $[4,5]$$[4,5]$[4,5][4,5]$[4,5]$. Ongoing research into the mechanical performance of reinforcing fibres, either synthetic or natural, continues to enhance the reliability and structural efficiency of FRP composites, ensuring their capability to sustain safety and performance under variable loading conditions [6,7]. Furthermore, the emergence of nanocomposites, characterized by superior mechanical and functional properties, opens new avenues for advanced structural applications and innovation in composite material technologies [8-10].
The integration of concrete, renowned for its superior compressive strength, with materials exhibiting high tensile properties is a well-established principle in structural engineering. Typical composite configurations include concrete-steel, concretetimber, and concrete-FRP systems, among others. Significant research efforts have been dedicated to optimising the geometry of continuous shear connectors to support efficacy and to enhance efficiency of the shear transfer mechanism at an interface with concrete. These studies aim to improve overall structural performance by refining connector configurations and characterising their mechanical behaviour under various loading conditions [11-13].
Mechanical analysis of the multi-component elements of civil engineering constructions is often accompanied by damage of materials and subsequent formation of cracks which appear either inside the components of the elements, or between them, at the interface [14-16]. It is very useful to develop computational methods for identification and analysis of development stages and interactions of cracks. These problems can be modelled using various physical approaches. One such an approach is the damage theory, which introduces an internal variable for describing a damage state [17,18].
To represent an evolution of states in a damageable structure, it is practical, especially for numerical solutions, to adopt an energy-based formulation. This formulation includes the energy balance of the damage process, which accounts for energy stored in the bodies, energy associated with newly formed defects, and dissipated energy [19]. In many cases, the processes under investigation can be treated as quasi-static, simplifying mathematical analysis and calculations. However, in general, the damage and crack growth process tend to occur rapidly, meaning inertial effects, as well as secondary damage and crack propagation caused by a kinetic energy transfer, may also be considered [20].
The interface cracks involve the contact between the components of the structural element. In the present case, the interface is modelled as an extremely thin adhesive layer. Degradation in this type of problem is referred to as cohesive or adhesive contact. Several computational models have been developed using cohesive zone models, e. g. as in references [21,22]. Some of them even account for friction which appears in interfaces which remain in contact after debonding, see [25]. From the computational point of view, a variational approach using energy principles, requires implementation of numerical
algorithms which optimise corresponding energy functionals. In the case of damage and crack problems with friction, one of the possible choices relies on application of general nonlinear programming algorithms, as in [23-26].
The presentation made here is a part of research studying slab bridges of short to medium spans with several different types of continuous shear connectors. In particular, the present paper introduces the previously developed computational cohesive models and procedures for an analysis of a shape parameter modification and a friction influence in jigsaw puzzle type shear connectors.

The currently designed and constructed deck bridges with encased filler-beams predominately contain rolled or welded steel I-sections. In the Centre for Research and Innovation in Construction of the Faculty of Civil Engineering at the Technical University of Košice, the first series of beam specimens with modified shapes of steel section was developed and tested. The individual variants differed in the method of ensuring composite action between the steel and concrete. Smooth and comb-like web edges of the T-section (Fig. 1) [27] were also compared. The efforts to better employ the steel section have led the authors to the idea of designing a GFRP section that could act mainly in the tension region of a future composite bridge structure. Different types of FRP section have been considered with the goal to design and experimentally verify deck bridges with modified FRP sections and achieve major economies.
A FRP beam makes use of a variety of methods of composite action using specially adjusted strip connectors. Various shapes of strip connectors can be seen in Fig. 2 [27].

In FRP-concrete composite girders, rib shear connectors can be used to transfer shear forces between the concrete slab and the FRP section. Rib shear connectors are produced in cutting regular recesses into the webs of FRP members. After encasing with concrete, the vertically embedded FRP dowels and the interstitial concrete dowels ensure a structural, interlocked connection (Fig. 3, a, b) [28].
The primary failure modes associated with adhesive bonding are [29]:
a) Adhesive failure,
b) Cohesive failure of adhesive, i.e. rupture of the glue,
c) Cohesive failure of adherent, i.e. rupture of the composite.

The present paper describes a computational model for interface damage which introduces the interface as a zone of adhesion. The interface damage finally causes debonding between material components. Debonded surfaces however partially remain in contact which is supposed to account for friction. Inside its procedures, the developed in-house MATLAB computational code [ $21,22,25$$21,22,25$21,22,2521,22,25$21,22,25$ ] discusses the influence of shape changes of the connector and effects of friction.

The computational model for interface damage developed in [25] was applied in a computational analysis of the connector, whose particular form is described below in Fig.4. The model is expressed as a problem of a quasi-static energy evolution in a multi-material domain $\mathrm{\Omega }$$\mathrm{\Omega }$Omega\Omega$\mathrm{\Omega }$. A brief description of the energy ingredients to the balance includes stored energy, energy of external forces and dissipated energy. The stored energy contains strain energy of the material domains, and it is described in terms of the displacement field $\mathit{u}$$\mathit{u}$u\boldsymbol{u}$\mathit{u}$. The interface is described a very thin elastic layer ${\mathrm{\Gamma }}_{\mathrm{c}}$${\mathrm{\Gamma }}_{\mathrm{c}}$Gamma_(c)\Gamma_{\mathrm{c}}${\mathrm{\Gamma }}_{\mathrm{c}}$, so its deformation is written as a function of a gap of the interface displacement $[[\mathit{u}]]$$[[\mathit{u}]]$[[u]]\llbracket \boldsymbol{u} \rrbracket$[[\mathit{u}]]$. Simultaneously, an interface damage variable $\zeta$$\zeta$zeta\zeta$\zeta$ characterising degradation of the interface is introduced together with a degradation function. The stored energy $E$$E$EE$E$ is than considered as follows

for an admissible displacement field $\mathit{u}$$\mathit{u}$u\boldsymbol{u}$\mathit{u}$ satisfying the displacement boundary conditions (constraints and prescribed displacement $\mathit{u}(t)$$\mathit{u}(t)$u(t)\boldsymbol{u}(t)$\mathit{u}(t)$ in Fig. 4) and admissible damage parameter $\zeta$$\zeta$zeta\zeta$\zeta$ ranged within the interval $[0,1]$$[0,1]$[0,1][0,1]$[0,1]$. Violate the state variables the constraints, the value of $E$$E$EE$E$ be unbounded. The introduced material parameters include stiffness matrix $\mathbf{C}$$\mathbf{C}$C\mathbf{C}$\mathbf{C}$ of the materials, stiffnesses $k_n,k_s,k_{\tau }$$k_n,k_s,k_{\tau }$k_(n),k_(s),k_(tau)k_{n}, k_{s}, k_{\tau}$k_n,k_s,k_{\tau }$ (normal, tangential, and transversal, respectively) of the interface layer and the compressive stiffness $k_g$$k_g$k_(g)k_{g}$k_g$ to penalise the contact between the material domains. The choice of the degradation function $\phi$$\phi$varphi\varphi$\phi$ affects the form of the interface stress-strain relation, generally it decreases from $\phi (0)=1$$\phi (0)=1$varphi(0)=1\varphi(0)=1$\phi (0)=1$ to $\phi (1)=0$$\phi (1)=0$varphi(1)=0\varphi(1)=0$\phi (1)=0$, but in the analysis a special form is considered: $\phi (\zeta )=\frac{\beta (1-\zeta )}{\beta +\zeta }$$\phi (\zeta )=\frac{\beta (1-\zeta )}{\beta +\zeta }$varphi(zeta)=(beta(1-zeta))/(beta+zeta)\varphi(\zeta)=\frac{\beta(1-\zeta)}{\beta+\zeta}$\phi (\zeta )=\frac{\beta (1-\zeta )}{\beta +\zeta }$, which leads to a bilinear stress-strain relation, with $\beta$$\beta$beta\beta$\beta$ determining the slope of its decreasing part. The term containing $k_{\mathrm{g}}$$k_{\mathrm{g}}$k_(g)k_{\mathrm{g}}$k_{\mathrm{g}}$ adds compressive stiffness (active only in compression where $[[\mathit{u}{]]}_n<0$$[[\mathit{u}{]]}_n<0$[[u]]_(n) < 0\llbracket \boldsymbol{u} \rrbracket_{n}<0$[[\mathit{u}{]]}_n<0$ ), which in a penalised form expresses normal contact condition.
The external forces, if present (c.f. the vertical pressure $p$$p$pp$p$ in Fig. 4), include the energy into balance. If force boundary conditions $\mathit{p}=(\mathbf{C}\mathbf{e}(\mathit{u}))\cdot \mathit{n}=\mathit{h}$$\mathit{p}=(\mathbf{C}\mathbf{e}(\mathit{u}))\cdot \mathit{n}=\mathit{h}$p=(Ce(u))*n=h\boldsymbol{p}=(\mathbf{C e}(\boldsymbol{u})) \cdot \boldsymbol{n}=\boldsymbol{h}$\mathit{p}=(\mathbf{C}\mathbf{e}(\mathit{u}))\cdot \mathit{n}=\mathit{h}$ are given on ${\mathrm{\Gamma }}_{\mathrm{N}}$${\mathrm{\Gamma }}_{\mathrm{N}}$Gamma_(N)\Gamma_{\mathrm{N}}${\mathrm{\Gamma }}_{\mathrm{N}}$, the pertinent functional for such a contour load is given by

The processes which dissipate the energy include the damage propagation and crack nucleation, especially due to its unidirectionality (damage only increases) expressed as $\dot{\zeta }\ge 0$$\dot{\zeta }\ge 0$zeta^(˙) >= 0\dot{\zeta} \geq 0$\dot{\zeta }\ge 0$. Then in problems of contact, there may appear friction, here considered in a standard Coulomb model, and, mainly due to computational reasons but also may be important physically, some simple rheology, e.g. Kelvin-Voigt. All these ingredients may be summarised in a dissipation (pseudo)potential as

where for any vector $\mathit{u}$$\mathit{u}$u\boldsymbol{u}$\mathit{u}$ the total tangential part is ${\mathit{u}}_t=(u_s,u_z)$${\mathit{u}}_t=\left(u_s,u_z\right)$u_(t)=(u_(s),u_(z))\boldsymbol{u}_{t}=\left(u_{s}, u_{z}\right)${\mathit{u}}_t=(u_s,u_z)$. The interface fracture energy $G_c$$G_c$G_(c)G_{c}$G_c$ determines the crack formation together with the condition of unidirectionality, otherwise the value of $R$$R$RR$R$ would be infinite. Additionally, the form of the fracture energy supposes its mode dependence, introducing the values $G_{\mathrm{c}}^{\mathrm{I}},G_{\mathrm{c}}^{\mathrm{II}},G_{\mathrm{c}}^{\mathrm{III}}$$G_{\mathrm{c}}^{\mathrm{I}},G_{\mathrm{c}}^{\mathrm{II}},G_{\mathrm{c}}^{\mathrm{III}}$G_(c)^(I),G_(c)^(II),G_(c)^(III)G_{\mathrm{c}}^{\mathrm{I}}, G_{\mathrm{c}}^{\mathrm{II}}, G_{\mathrm{c}}^{\mathrm{III}}$G_{\mathrm{c}}^{\mathrm{I}},G_{\mathrm{c}}^{\mathrm{II}},G_{\mathrm{c}}^{\mathrm{III}}$ inside the formula

The second term introduces dissipation due to friction. The matrix $M$$M$MM$M$ defines an orthotropic friction, where the coefficient of friction may depend on mutual orientation of material axes of symmetry due to their inner structure. The dependence on damage may reflects how the state of the adhesiveness changes the friction coefficient. The orientation of the tangential stress ${\mathit{p}}_t$${\mathit{p}}_t$p_(t)\boldsymbol{p}_{t}${\mathit{p}}_t$ is determined by the vector $\mathrm{M}(\zeta )[[\dot{\mathit{u}}{]]}_t$$\mathrm{M}(\zeta )[[\dot{\mathit{u}}{]]}_t$M(zeta)[[u^(˙)]]_(t)\mathrm{M}(\zeta) \llbracket \dot{\boldsymbol{u}} \rrbracket_{t}$\mathrm{M}(\zeta )[[\dot{\mathit{u}}{]]}_t$ and while sliding, the magnitude be proportional to the normal compressive stress $\left|p_n\right|=\left|k_{\mathrm{g}}\left([[\mathit{u}{]]}_n^{-}\right)\right|$$\left|p_n\right|=\left|k_{\mathrm{g}}\left([[\mathit{u}{]]}_n^{-}\right)\right|$|p_(n)|=|k_(g)(([[)u]]_(n)^(-))|\left|p_{n}\right|=\left|k_{\mathrm{g}}\left(\llbracket \boldsymbol{u} \rrbracket_{n}^{-}\right)\right|$\left|p_n\right|=\left|k_{\mathrm{g}}\left([[\mathit{u}{]]}_n^{-}\right)\right|$.
Finally, the last term appears due to the Kelvin-Voigt rheology of the material, providing stress as $\mathit{\sigma }=\mathbf{C}\mathit{e}+\mathbf{D}\dot{\mathit{e}}$$\mathit{\sigma }=\mathbf{C}\mathit{e}+\mathbf{D}\dot{\mathit{e}}$sigma=Ce+De^(˙)\boldsymbol{\sigma}=\mathbf{C} \boldsymbol{e}+\mathbf{D} \dot{\boldsymbol{e}}$\mathit{\sigma }=\mathbf{C}\mathit{e}+\mathbf{D}\dot{\mathit{e}}$, to be considered so that $\mathbf{D}={\tau }_{\mathrm{r}}\mathbf{C}$$\mathbf{D}={\tau }_{\mathrm{r}}\mathbf{C}$D=tau_(r)C\mathbf{D}=\tau_{\mathrm{r}} \mathbf{C}$\mathbf{D}={\tau }_{\mathrm{r}}\mathbf{C}$ for a relaxation time parameter ${\tau }_{\mathrm{r}}$${\tau }_{\mathrm{r}}$tau_(r)\tau_{\mathrm{r}}${\tau }_{\mathrm{r}}$.
The quasi-static evolution is then controlled by the following system of nonlinear variational inclusions.

The form of initial conditions assumes that damage starts from an intact state. The first relation determines the stress equilibrium and assumes also solution with friction. The other relation provides the flow rule for damage evolution. In particular case of energy form Eq. (1) and dissipation potential (3) it can provide the damage condition pointwisely by an inequality

which initiates degradation, when, at a pristine state ( $\zeta =0$$\zeta =0$zeta=0\zeta=0$\zeta =0$ ), the interface strain expressed in terms of the displacement gap reaches the equation, modified to provide a stress limit for interface damaging

The computational code derived for the described model is written within an in-house MATLAB code, as described in [25]. See the details therein, though it should be stated that the code advantageously utilises discretisation by boundary element method as the unknowns in the nonlinear evolution process are gathered along the boundaries, and in particular along the interface. The evolution is expressed by a time stepping algorithms which is described by a staggered approach. Supposing at least some amount of rheology, allows to consider friction from the previous time step which simplifies calculation so that the relations in Eq. (4) can be seen as minimisation conditions for convex functionals. It allows to use sequential quadratic algorithms for such minimisation problems.

The model and its MATLAB implementation as described in [25] were used in calculations. Fig. 4 shows a scheme which served for a simple parametric study performed on the connector. The parameters which were being varied in the calculation include the shape of the connector and friction coefficient along the interface after debonding the concrete and FRP parts. The shape of the connector includes three values: width $d$$d$dd$d$ and two radii of curvature. To simplify the geometry $d$$d$dd$d$ was kept fixed $d=50\mathrm{mm}$$d=50\mathrm{mm}$d=50mmd=50 \mathrm{~mm}$d=50\mathrm{mm}$, and the two radii were kept the same. Thus, there remained one geometric parameter $r$$r$rr$r$. For analysing the friction influence we supposed that the material orientation is along the $x_1$$x_1$x_(1)x_{1}$x_1$ axis, in the direction of the load. So that only a scalar parameter was sufficient, with matrix $M$$M$MM$M$ from Eq. (3) defined as $\mathbf{M}(\zeta )=\mu {\zeta }^4\mathrm{I}$$\mathbf{M}(\zeta )=\mu {\zeta }^4\mathrm{I}$M(zeta)=muzeta^(4)I\mathbf{M}(\zeta)=\mu \zeta^{4} \mathrm{I}$\mathbf{M}(\zeta )=\mu {\zeta }^4\mathrm{I}$ , with the unit matrix I. This choice stresses that the Coulomb friction does not play a role while the interface is fully adhesive. The scheme also shows the prescribed displacement $\mathit{u}$$\mathit{u}$u\boldsymbol{u}$\mathit{u}$ linearly increasing during loading as $\mathit{u}(t)=0.001t(1,0,0)[\mathrm{mm}]$$\mathit{u}(t)=0.001t(1,0,0)[\mathrm{mm}]$u(t)=0.001 t(1,0,0)[mm]\boldsymbol{u}(t)=0.001 t (1,0,0)[\mathrm{mm}]$\mathit{u}(t)=0.001t(1,0,0)[\mathrm{mm}]$ with pseudo time step $t=1$$t=1$t=1t=1$t=1$. To simulate the weight of the above material, the constant vertical pressure for Eq. (2) $\mathit{h}=\mathit{p}=(0,5,0)\mathrm{kPa}$$\mathit{h}=\mathit{p}=(0,5,0)\mathrm{kPa}$h=p=(0,5,0)kPa\boldsymbol{h}=\boldsymbol{p}=(0,5,0) \mathrm{kPa}$\mathit{h}=\mathit{p}=(0,5,0)\mathrm{kPa}$ was applied. Both loads are schematically shown in Fig. 4 (right), too.
The material parameters (as a part of the stiffness matrix $\mathbf{C}$$\mathbf{C}$C\mathbf{C}$\mathbf{C}$ ) are considered as follows: the upper concrete part has the following data (considered as an isotropic material) $E=32\mathrm{GPa},v=0.2$$E=32\mathrm{GPa},v=0.2$E=32GPa,v=0.2E=32 \mathrm{GPa}, v=0.2$E=32\mathrm{GPa},v=0.2$, the lower FRP part (glass fibres) is supposed as a transversally isotropic material with characteristics: $E_1=39\mathrm{GPa},E_2=E_3=8.6\mathrm{GPa},G_{12}=G_{13}=3.8\mathrm{GPa},v_{12}=v_{13}=0.28,v_{23}=0.35$$E_1=39\mathrm{GPa},E_2=E_3=8.6\mathrm{GPa},G_{12}=G_{13}=3.8\mathrm{GPa},v_{12}=v_{13}=0.28,v_{23}=0.35$E_(1)=39GPa,E_(2)=E_(3)=8.6GPa,G_(12)=G_(13)=3.8GPa,v_(12)=v_(13)=0.28,v_(23)=0.35E_{1}=39 \mathrm{GPa}, E_{2}=E_{3}=8.6 \mathrm{GPa}, G_{12}=G_{13}=3.8 \mathrm{GPa}, v_{12}=v_{13}= 0.28, v_{23}=0.35$E_1=39\mathrm{GPa},E_2=E_3=8.6\mathrm{GPa},G_{12}=G_{13}=3.8\mathrm{GPa},v_{12}=v_{13}=0.28,v_{23}=0.35$. The interface data are adjusted according to empirical data obtained by [30] which assessed different adhesion conditions between GFRP and concrete. Relying on that the interface fracture energy took the value $G_{\mathrm{c}}=5\mathrm{mN}{\mathrm{m}}^{-1}$$G_{\mathrm{c}}=5\mathrm{mN}{\mathrm{m}}^{-1}$G_(c)=5mNm^(-1)G_{\mathrm{c}}=5 \mathrm{mN} \mathrm{m}^{-1}$G_{\mathrm{c}}=5\mathrm{mN}{\mathrm{m}}^{-1}$, and the stiffness was set to fit the material stiffness by putting $k=1\mathrm{TPa}{\mathrm{m}}^{-1}$$k=1\mathrm{TPa}{\mathrm{m}}^{-1}$k=1TPam^(-1)k=1 \mathrm{TPa} \mathrm{m}^{-1}$k=1\mathrm{TPa}{\mathrm{m}}^{-1}$. The degradation function $\phi$$\phi$varphi\varphi$\phi$, accounting for the effect of fibre bridging, was considered with $\beta =0.1$$\beta =0.1$beta=0.1\beta=0.1$\beta =0.1$. These adjustments provided the critical normal interface stress $p_n=9.5\mathrm{kPa}$$p_n=9.5\mathrm{kPa}$p_(n)=9.5kPap_{n}=9.5 \mathrm{kPa}$p_n=9.5\mathrm{kPa}$, based on the relation (5) expressed in stress variables. The friction coefficient was varied within the range $[0,2]$$[0,2]$[0,2][0,2]$[0,2]$. Thought, the actual measurements confirm the friction somewhere in the range $[0.5,0.8]$$[0.5,0.8]$[0.5,0.8][0.5,0.8]$[0.5,0.8]$, the demonstration is also done for extremely high friction, for small frictions, and for the frictionless case.

The changes of geometry assume $r$$r$rr$r$ between 7.25 mm and 20 mm . The actual changes in the shape of the connector due to modification of $d$$d$dd$d$ are shown in Fig. 5. It might be expected that for a large $r$$r$rr$r$ the upper part tends to slide along the interface, while for a small $r$$r$rr$r$ such sliding is inhibited and at least at some points the interfaces would not be damaged.

The curves in Fig. 6 reflect that expectation and also demonstrate other key features. They show the resultant horizontal force applied at the right face of the upper material. The graphs are pertinent to two friction options: no friction and friction with $\mu =0.8$$\mu =0.8$mu=0.8\mu=0.8$\mu =0.8$. The first part of the relation pertains to full adhesion, and it is almost the same for all friction and geometry options. Small differences were influenced by variations in the overall stiffness due to connector modifications (various $r$$r$rr$r$ ). Then, up to the first local minimum (or change of slope in the frictional case), debonding at the right end of the interface occurs. That interface part slides (friction helps to transfer the horizontal force), and finally (as shown below in Figs. 8,9) a crack is opened along the left side of the connector. Subsequently, the data start to be different as the next debonding is affected not only by the friction but also by the connector's shape. Anyhow, the second weakening part is connected with debonding at the left end of the interface. Finally, at least in the frictionless case, the whole interface is fully debonded and the two material domains hold together due to contact forces. With friction, as shown below in Figs. 11, 12, a part of interface around the compression zone is undamaged, basically caused right by the friction. Particular effect of the friction for a fixed connector shape is shown in Fig. 7. The hardening effect of friction is natural. Although, the friction does not prevent from opening cracks, it significantly helps the connector to improve the load bearing capacity of the compound structural element [31].

To support the previous statements, graphs displaying deformation of the structural element and some interface variables are provided. Debonding observed in deformations is shown in Fig. 8. Some instants which were important in Fig. 6 are used now to see the sliding initialisation of debonding at $u=0.1\mu \mathrm{m}$$u=0.1\mu \mathrm{m}$u=0.1 mumu=0.1 \mu \mathrm{~m}$u=0.1\mu \mathrm{m}$, interface crack opening along the right part of the connector at $u=0.14\mu \mathrm{m}$$u=0.14\mu \mathrm{m}$u=0.14 mumu=0.14 \mu \mathrm{~m}$u=0.14\mu \mathrm{m}$, opening the left end of the interface at $u=0.3\mu \mathrm{m}$$u=0.3\mu \mathrm{m}$u=0.3 mumu=0.3 \mu \mathrm{~m}$u=0.3\mu \mathrm{m}$, and decohesion along whole but compressive part of the interface (close to the point $B$$B$BB$B$ of the interface according to Fig. 5) at $u=0.51\mu \mathrm{m}$$u=0.51\mu \mathrm{m}$u=0.51 mumu=0.51 \mu \mathrm{~m}$u=0.51\mu \mathrm{m}$. As friction affects the force relation, it also has an effect on the deformation. The final instant of the applied load pertains to $u=0.8\mu \mathrm{m}$$u=0.8\mu \mathrm{m}$u=0.8 mumu=0.8 \mu \mathrm{~m}$u=0.8\mu \mathrm{m}$. The deformation at this instant for three friction options is shown in Fig. 9. Though it is not so large, the difference may be observed at the deformation close to the topmost face of the connector (notice the interference with the text label). It is also observable at the vicinity of the right end point of the interface: the shape change of the bottom domain.

Graphs in Figs. 10-12 contain interface distributions of contact variables. Three friction options (taken of those used in Fig. 7) document how interface damage evolves documented by the same selected instants (plus an initial one) as shown above. In all cases, distribution of damage variable shows where the interface is broken and how friction affects the tractions along the interface where broken. The location of particular interface point is referred to by the letters introduced in Fig. 5.

First, compare the propagation of interface damage in all figures. Naturally, initiation is very similar, but later as friction is initialised, differences appear so that higher friction blocks or even prevents propagation of damage along the whole interface. Of course, it plays the role only in the zones of compression (check the situation close to the point $B$$B$BB$B$ ).

The effects on the contact forces can be compared in the left graphs of the figures. The range of the tangential stress is scaled by friction (in cases of nonzero friction), so that after debonding and during sliding the normal and tangential components, which obey the Coulomb rule, graphically coincide. While in the surroundings of the point $E$$E$EE$E$ the stress distributions are very similar due to opening, the situation in sliding zones is different. Check e. g. the right interface end: initially the stress distribution is friction independent, anyhow, after damage triggering, friction plays a significant role in distribution of tangential stresses. In the frictionless case the tangential stress vanishes, in frictional cases it obeys its own rule. The compressive stress in more frictional case is not so big as in the no-friction case as the horizontal components are more distributed.
The instants and strength of influence of all these phases are affected by the value of the parameter $r$$r$rr$r$ as can be guessed in. The result which concludes the parametric study is that small $r$$r$rr$r$ guarantees higher stiffness of the connector also after debonding between FRP and concrete. Higher friction than plays a positive role on structural stiffness. It was interesting to see that the design of the connector made without this computational analysis, and which supposed $r=8.5\mathrm{mm}$$r=8.5\mathrm{mm}$r=8.5mmr=8.5 \mathrm{~mm}$r=8.5\mathrm{mm}$ is confirmed computationally, where the best force transfer in the cases of moderate friction was obtained for $r$$r$rr$r$ approximately 9 mm .

Ashear jigsaw-puzzle shaped connector experimentally analysed in published materials was studied computationally using an interface damage model accounting also for frictional response. A simple parametric study covered two distinct aspects of the connector. One of them was a shape parameter - geometrical. The other was physical and included the influence of friction as a characteristic of roughness of the contacting surfaces. An in-house MATLAB code provided the implementation of the cohesive interface model. The results demonstrate how the compound materials may be improved by adjusting dimensions of the connector and by appropriate preparing of surfaces to result in a reliable structural component. Simultaneously it showed that these calculations are in accordance with experimental observations. Anyhow, for anisotropic materials with fibre structure, it is also important to know how the analysed properties vary with material orientation axes, which was not tested here. Therefore, the future plans for computational analysis of material connectors cover this aspect of anisotropy.

We would like to acknowledge the support by the Scientific Grant Agencies of the Slovak Republic under projects VEGA 1/0307/23, VEGA 1/0365/25 and by Slovak Research and Development Agency under the project APVV-23-0204.

[1] Sonnenschein, R., Gajdosova, K. and Holly, I. (2016). FRP composites and their using in the construction of bridges, Proc. Eng., 161, pp. 477-482. DOI: https://doi.org/10.1016/j.proeng.2016.08.665.
[2] Siwowski, T. Rajchel, M. and Kulpa, M. (2019). Design and field evaluation of a hybrid FRP composite-lightweight concrete road bridge, Comp. Struct., 230, 111504. DOI: https://doi.org/10.1016/j.compstruct.2019.111504.
[3] Chróscielewski, J., Miskiewicz, M., Pyrzowski, Ł., Sobczyk, B. and Wilde, K. (2017). A novel sandwich footbridgePractical application of laminated composites in bridge design and in situ measurements of static response, Comp. Part B, 126, pp. 153-161. DOI: https://doi.org/10.1016/j.compositesb.2017.06.009.
[4] Naser, M. Z., Hawileh, R. A. and Abdalla, J. (2021). Modeling strategies of Finite Element simulation of reinforced concrete beams strengthened with FRP: A review, Journal of Composites Science, 5(19).
DOI: https://doi.org/10.3390/jcs5010019.
[5] Kormaníková, E., Źmindák, M. and Sabol, P. (2021). Tensile properties of carbon fiber reinforced polymer matrix composites: Application for the strengthening of reinforced concrete structure. Composite Structures, 275, pp. 1-12. DOI: https://doi.org/10.1007/978-3-030-27053-7_64.
[6] Dixit, S., Goel, R., Dubey, A., Ahivhare, P.R. and Bhalavi, T. (2017). Natural fibre reinforced polymer composite materials. A review. Polym. Renew. Resour., 8, pp. 71-78. DOI: https://doi.org/10.1177/204124791700800203.
[7] Rajak, D. K., Pagar, D. D., Kumar, R. and Pruncu, C. (2019). Recent progress of reinforcement materials: A comprehensive overview of composite materials. J. Mater. Res. Technol. 8, 6354-6374.
DOI: https://doi.org/10.1016/j.jmrt.2019.09.068.
[8] Aziz, T., Fan, H., Haq, F., Khan, F.U., Numan, A., Iqbal, M., Raheel, M., Kiran, M. and Wazir, N. (2020). Adhesive properties of poly (methyl silsesquioxanes)/bio-based epoxy nanocomposites. Iran. Polym. J. 29, pp. 911-918. DOI: https://doi.org/10.1007/s13726-020-00849-x.
[9] Jamil, M. I., Wang, Q., Ali, A., Hussain, M., Aziz, T., Zhan, X. and Zhang, Q. (2021). Slippery photothermal trap for outstanding deicing surfaces. J. Bionic Eng., 18, pp. 548-558. DOI: https://doi.org/10.1007/s42235-021-0046-7.
[10] Ahmed, N., Fan, H., Dubois, P., Zhang, X., Fahad, S., Aziz, T. and Wan, J. (2019). Nano-engineering and micromolecular science of polysilsesquioxane materials and their emerging applications. J. Mater. Chem. A, 7(1), pp. 21577-21604. DOI: https://doi.org/10.1039/C9TA04575A.
[11] Zou, Y., Guo, J., Zhou, Z., Wang, X., Yu, Y. and Zheng, K. (2022). Evaluation of shear behavior of PCSC shear connection for the construction of composite bridges with prefabricated decks. Engineering Structures, 257, 113870. DOI: https://doi.org/10.1016/j.engstruct.2022.113870.
[12] Vanova, P. and Kvocak, V. (2021). Experimental and FEM analysis of innovative bridge design. Transportation Research Procedia, 55, pp. 1188-1193. DOI: https://doi.org/10.1016/j.trpro.2021.07.099.
[13] Zhang, J., Hu, X., Kou, L., Zhang, B., Jiang, Y. and Yu, H. (2018). Experimental study of the short-term and longterm behavior of perfobond connectors. Journal of Constructional Steel Research, 150, pp. 462-474.
DOI: https://doi.org/10.1016/j.jcsr.2018.09.004.
[14] Feklistova, E. V., Mugatarov, A. I., Wildemann, V. E. and Agishev, A. A. (2024). Fracture processes numerical modeling of elastic-brittle bodies with statistically distributed subregions strength values, Frattura ed Integrità Strutturale, 68, pp. 325-339. DOI: https://doi.org/10.3221/IGF-ESIS.68.22.
[15] Markides, Ch. F. and Kourkoulis, S. K. (2025). Revisiting classical concepts of Linear Elastic Fracture Mechanics Part III: The stress field in a double-edge notched finite strip by means of the "stress-neutralization" technique, Fracture and Structural Integrity, 71, pp. 302-316. DOI: https://doi.org/10.3221/IGF-ESIS.71.22.
[16] Kormanikova, E. and Kotrasova, K. (2025). Numerical study of delamination process of the CFRP composite, Fracture and Structural Integrity, 71, pp. 182-192. DOI: https://doi.org/10.3221/IGF-ESIS.71.13.
[17] Capozucca, R., Magagnini, E., Vecchietti, M.V. and Khatir, S. (2021). RC beams damaged by cracking and strengthened with NSM CFRP/GFRP rods, Frattura ed Integrità Strutturale, 58, pp. 386-401.
DOI: https://doi.org/10.3221/IGF-ESIS.58.28.
[18] Maugin, G. A (2015). The saga of internal variables of state in continuum thermo-mechanics (1893-2013). Mechanics Research Communications 69, pp. 79-86. DOI: https://doi.org/10.1016/j.mechrescom.2015.06.009.
[19] Bourdin, B., Francfort, G. A. and Marigo, J.-J. (2008). The variational Approach to Fracture, J. Elasticity, 91, pp. 5148. DOI: https://doi.org/10.1007/s10659-007-9107-3.
[20] Roubíček, T. (2019). Topics in Applied Analysis and Optimisation: Partial Differential Equations, Stochastic and Numerical Analysis, CIM Series in Math. Sci., pp. 363-396, Springer.
DOI: https://doi.org/10.1007/978-3-030-33116-0.
[21] Vodička, R. (2016). A quasi-static interface damage model with cohesive cracks: SQP-SGBEM implementation, Engineering Analysis with Boundary Elements, 62(3), pp. 123-140.
DOI: https://doi.org/10.1016/j.enganabound.2015.09.010.
[22] Vodička, R. and Mantič, V. (2017). An energy-based formulation of a quasi-static interface damage model with a multilinear cohesive law, Discrete and Cont. Dynam. Syst.-S, 10(6), pp. 1539-1561.
DOI: https://doi.org/10.3934/dcdss. 2017079.
[23] Nocedal, J., Wright, S.J. (2006). Numerical Optimisation, Springer Series in Operation research, Second Edition.
[24] Kormaníková, E., Kšiňan, F. and Vodička, R. (2024). Computational and experimental approaches for Mode I delamination problems, International Journal of Solids and Structures, 300, 112926, DOI: https://doi.org/10.1016/j.ijsolstr.2024.112926.
[25] Vodička, R., Kormaníková, E. and Kšiňan, F. (2018). Interfacial debonds of layered anisotropic materials using a quasistatic interface damage model with Coulomb friction, International Journal of Fracture, 211 (1-2), pp. 163-182. DOI: https://doi.org/10.1007/s10704-018-0281-z.
[26] Blázquez, A., Mantič, V., París, F. and McCartney, L.N. (2008). Stress state characterization cracks in [0/90] symmetric laminates by BEM. International Journal of Solids and Structures, 45, pp. 1632-1662.
DOI: https://doi.org/10.1016/j.ijsolstr.2007.10.013.
[27] Dubecký, D. and Kvočák, V. (2021). Research and Development of Deck Bridges, Springer Nature, ISBN 978-3-030-66924-9.
[28] Classen, M., Herbrand, M., Adam, V., Kueres, D. and Sarac, M. (2018). Puzzle-shaped rib shear connectors subjected to combined shear and tension, Journal of Constructional Steel Research, 145, pp. 232-243.
DOI: https://doi.org/10.1016/j.jcsr.2018.02.036.
[29] Zoghi, M. (2013). The international handbook for FRP composite in civil engineering. CRC Press Taylor & Francis Group, U.S, pp. 706.
[30] Kvočák, V., Dubecký, D. and Weissová, M. (2024). Experimental measurement of friction coefficient between concrete and composite material (in Slovak), In: Scientific and research activities of the Institute of Structural Engineering and Transportation Structures, Technical University of Košice, pp. 65-68.
[31] Weissova, M., Kvocak, V., Dubecky, D., Vanova P. (2022). Numerical analysis of shear resistance on composite beams. IOP Conf. Series: Materials Science and Engineering, 1252, 012035.
DOI: https://doi.org/10.1088/1757-899X/1252/1/012035.