Inter-fibre failure under biaxial loads in glass–epoxy composite materials: Effect of the presence of a nearby fibre

Theoretical and Applied Fracture Mechanics 127 (2023) 103956 Available online 22 June 2023 0167-8442/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Inter-fibre failure under biaxial loads in glass–epoxy composite materials: Effect of the presence of a nearby fibre Carlos Sandino * , Elena Correa , Federico París Group of Elasticity and Strength of Materials, School of Engineering, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Seville, Spain ARTICLE INFO Keywords: Composites Micromechanics Interfacial Fracture Mechanics Matrix/inter-fibre failure Biaxial loads Two-fibre model ABSTRACT Fibre-reinforced composite materials are especially prone to transverse failure. It appears at the lamina level following the mechanism known as matrix/inter-fibre failure. This mechanism of damage is associated with the appearance of fibre–matrix debonds (interface cracks), as shown in previous numerical micromechanical studies. After the nucleation, growth and kinking into the matrix, these interface cracks give rise to the final macrofailure. When compared to uniaxial loading, the growth stages of this mechanism of damage (analysed in light of Interfacial Fracture Mechanics) show some alterations under different combinations of biaxial loads. This work gives a step forward and focuses on the micromechanical BEM study of the evolution of an interface crack in the presence of a neighbouring fibre. Thus, after considering a transverse tensile load (nominally responsible for the failure) a secondary transverse load is also applied (tensile or compressive, perpendicular to the primary load). When considering the two-fibre BEM model, the results obtained lead to identifying the neighbouring fibre locations that act as accelerative agents on failure progression and establishing the effect of the biaxial load on them. Specifically, when a secondary tensile load is applied, the presence of the nearby fibre (for most of its positions) confirms the slight inhibition of the mechanism of failure for biaxial tensile loads already referred to in previous single-fibre studies by the authors. As the secondary tensile load increases, it tends to mitigate the effect of the presence of the neighbouring fibre that was previously observed for uniaxial tensile load. The opposite effects are found when a secondary compressive load is considered, which intensifies the alterations of the presence of the neighbouring fibre on the interface crack growth. Experimental evidence on some aspects is provided confirming the associated conclusions derived from the numerical models. 1. Introduction The study of the damage mechanism taking place in long fibrereinforced composites at the micromechanical level is fundamental for the improvement of industrial applications involving these materials. In this regard, transverse failure, which starts with small interfacial debonds that lead to a transverse crack that produces the failure at the ply level, plays a key role. This failure is generated by the damage mechanism known as inter-fibre or matrix failure, which is due to loads perpendicular to the direction of the fibres in unidirectional laminates or impact loads in multidirectional laminates. The importance of this mechanism lies in the fact that it usually defines the damage initiation in any laminate. At the micromechanical level, the inter-fibre failure has been numerically analysed under uniaxial tension [1,2] and under uniaxial compression [3–5], by means of BEM single-fibre models. In these studies, based on the hypothesis that transverse failure starts with the appearance of debonds at the interface between the fibres and the surrounding matrix, the different stages of the damage mechanism that lead to the macro-failure were identified. In the case of the initiation of the failure due to a transverse tensile load, three stages were considered, as shown in Fig. 1: Stage 1. Crack nucleation: the initiation of the damage mechanism is assumed to be controlled by the radial stress at the interface. For a single-fibre case, the maximum value of the radial stress responsible for the failure is detected at 0◦from the direction of the uniaxial tensile load, resulting in the appearance of the first debond. Stage 2. Interface crack growth: starting from the initial debond, the * Corresponding author. E-mail address: [email protected] (C. Sandino). Contents lists available at ScienceDirect Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec https://doi.org/10.1016/j.tafmec.2023.103956 Received 26 March 2023; Received in revised form 21 April 2023; Accepted 5 June 2023 Theoretical and Applied Fracture Mechanics 127 (2023) 103956 2 crack grows unstably and symmetrically along the interface until it reaches a certain length (60◦-70◦) from the direction of the external load. The end of this stage coincides with the appearance of a finite contact zone at the crack tips. Stage 3. Kinking: the crack growth becomes stable, which favours a change in the direction of propagation, the crack kinking into the matrix in the direction perpendicular to the tensile load applied. This crack continues growing through the matrix, where the coalescence of cracks leads to failure at the lamina level. After the identification of the differentiated stages of the damage mechanism and the use of single-fibre models, some authors have included a greater number of fibres in their numerical tools (BEM or FEM). Such is the case of the studies presented in [6,7] for uniaxial tensile and compressive loads respectively. These works are focused on the BEM analysis of the interface crack growth in the presence of a nearby undamaged fibre. Thus, a two-fibre model is employed for analysing the influence of the presence of the nearby fibre and its position on the propagation of the crack. Multifibre FEM models are presented in [8], using cohesive elements, and in [9], applying Linear Elastic Interface Models (LEBIM). The study of this damage mechanism has also been analysed from the point of view of Finite Fracture Mechanics. Such is the case of [10] for the tensile uniaxial case. In [11], the authors proposed a model based on Finite Fracture Mechanics for the study of the interface debonding under a combination of plane tension and out-of-plane shear. Some other numerical studies have focused on the initial stage of the damage, such is the case of [12,13], whereas other authors have analysed the state previous to the appearance of the macro-crack [14,15]. In [16,17], different models of Representative Volume Elements have been implemented, in order to study the influence of the fibre volume fraction, the distance between interface cracks, the distance to the free surface or the distance to the 0◦/90◦interface in cross-ply laminates, on the evolution of the Energy Release Rate. All these approaches are especially appropriate for the study of fibre-reinforced composite laminates containing ultra-thin plies, in order to have a better understanding of the phenomenon known as the scale effect ([18–20]). Regarding the experimental approaches of this problem, some authors have studied the damage mechanism at the micromechanical level. These studies included the outstanding works in the observation of interface cracks presented in [21] for static loads, and in [22] for dynamic loads. The presence of thermal residual stresses, which represent an important aspect to be considered in the curing process of fibrous composites, is evaluated in [23]. Other studies include the effects of the fibre type [24], high and low curing temperatures [25], the thickness of 90◦plies in cross-ply laminates [26], or the use of X-ray microtomography in single-fibre specimens to analyse the initial debond and the propagation of the interface crack [27]. Recent approaches include the works presented in [28], which presents optical microscopy measurements of the different stages of the damage mechanism under tension in cross-ply laminates, and in [29], which follows the same approach for the compressive damage mechanism. With respect to the aforementioned stages of the damage mechanism, some authors have considered the possible alterations that biaxial transverse loading could cause. It has been historically assumed that, when considering biaxial transverse loads, inter-fibre failure occurs in a plane governed by the components of the stress vector associated with the referred plane [30]. This fact implies that, for instance, a secondary load transverse to a main tensile load would not be relevant in the generation of the interface crack. In order to confirm this hypothesis in the case of a single-fibre model under biaxial loads, the studies carried out in [31,32] analysed the possible effects of a secondary transverse load on the initiation and propagation of an interface crack governed by tensile or compressive loads. In order to approach the problem from an experimental point of view, the design of cruciform composite specimens for biaxial testing was carried out in [33]. More specifically, four different designs of cruciform specimens consisting of 90◦unidirectional composite laminates were proposed and studied for different load levels, using FEM models. This analysis was validated by the results of the study presented in [34], in which, the design of a biaxial device adapted for use in a universal uniaxial testing machine was carried out to perform experimental biaxial tests. Subsequently, based on the study presented in [33], the design and manufacture of cruciform specimens and the development of an experimental testing campaign were carried out in [35]. In the present study, the authors follow the approach of the two-fibre models presented in [6,7] for uniaxial tension and compression and the analysis of biaxial loading for a single-fibre case presented in [31,32]. To this aim, a two-fibre BEM model is used to observe the influence of a transverse (tensile or compressive) secondary load on the interface crack growth when a tensile main load is considered. Thus, the alterations produced by the presence of a neighbouring undamaged fibre are analysed by modifying the value of the secondary load and comparing the results with the uniaxial tensile case [6]. The analysis of all the cases under consideration is focused on the evolution of the Energy Release Rate during the interface crack growth and its use to predict some aspects such as the length and morphology of the crack. This work is structured as follows: Section 2 is devoted to the description of the BEM numerical model; in Section 3 the tension–tension biaxial cases (T-nT) are studied; and Section 4 corresponds to the tension–compression biaxial cases (T-nC); in Section 5, the most important results are summarised and discussed, and they are compared with experimental evidence. Finally, Section 6 is devoted to the main conclusions of this study. 2. Numerical model In this work, the numerical study was carried out using a tool based on BEM [36] and developed in [37]. This BEM code allows the numerical analysis of plane elastic problems for contact and interface cracks to be performed. A scheme representing the model employed in this work is shown in Fig. 2. It represents the case of a crack that grows (under the plane strain hypothesis) along the interface between the matrix and the so-called primary fibre, in the presence of a secondary undamaged fibre. The appearance of the first debond is assumed to be associated with the direction of the primary tensile load for all cases under consideration, even though the secondary load and the position of the secondary fibre could alter its position. This assumption is reinforced by the results obtained in [35] where in all T-nT and T-nC cases the macro-failure always occurred perpendicular to the main external tension. In order to compare the results of this study with the biaxial single-fibre case [31] and the uniaxial two-fibre case [6], a symmetrical model is established. Thus, the presence of the secondary fibre produces the same effects at both crack tips. Notice that now the entire model consists of three fibres. Solids ‘M’, ‘PF’ and ‘SF’ represent the matrix and the primary/secondary fibres respectively. The fibres’ size is defined by the radius r 1 Fig. 1. Micromechanical phases of the inter-fibre failure identified with a single-fibre model under uniaxial tension. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 3 (r 1 =7.5⋅10 -6 m) and the matrix is large enough in order to avoid the border effect (b=100⋅r 1 and a=b). The position of the secondary fibre is defined by parameters r 2 and θ 2 and is modified to analyse different relative positions between the fibres. For this study, the initial distance employed between the fibres is r 2 =r 2 0 =2.416⋅r 1 , which corresponds to the distance calculated for a fibre volume fraction V f =62% and hexagonal packing. With regard to the loads applied, the different cases follow the notation T-nT/T-nC, where the first T represents the transverse tension nominally responsible for the failure ( σ 0 , aligned with the x-axis) and the term nT/nC corresponds to the secondary tensile or compressive load (aligned with the y-axis). This secondary load is a fraction of that applied on the x-axis, this fraction being defined by coefficient n. The n values considered in this study are n =0, 0.25, 0.5, 0.75, 1. The materials chosen for the analysis correspond to a glass fibreepoxy matrix system. The elastic properties are included in Table 1. To describe the problem from the Fracture Mechanics point of view, the Energy Release Rate, G, is used. It is calculated following Eq. (1) (based on [38]), for a circular crack that propagates from a certain debonding angle, θ d , Fig. 2, to θ d +Δθ d (Δθ d ≪θ d ): G(θd,Δθd) = 1 2Δθd∫θd+Δθd θd [ σ rr(θd+θ)Δur(θd−Δθd+θ) + σ rθ(θd+θ)Δuθ(θd−Δθd+θ) ]dθ (1) where θ is the circumferential coordinate with reference to the x-axis. σ rr and σ rθ represent the radial and shear stresses along the interface, and Δu r and Δu θ represent the relative displacements of the crack faces. Modes of fracture I and II are obviously considered in Eq. (1), associated with σ rr and σ rθ respectively. For this study, the value of Δθ d employed is 0.5◦. Dimensionless results for G will be presented. These values are obtained, following [39,40], by dividing the values of G by G 0 : G0=1+κm 8 μ m σ 2 0r1 π (2) where κ m =3-4 ν m and μ m is the shear modulus of the matrix. 3. T-nT Cases 3.1. Energy Release Rate In this Section, the Energy Release Rate, G, is calculated from Eq. (1) for a growing interface crack. This allows the problem to be characterised from the Interface Fracture Mechanics point of view. Considering the model presented in Fig. 2, the damage was assumed to appear as an initial debond of θ d =10◦[6,31]. Then, G was obtained for a debonding interval 10◦≤θ d ≤150◦, using the single-fibre model (reference case) and the two-fibre model (with different relative positions between the fibres), and applying the biaxial tensile cases T-nT. Regarding the twofibre case, for any biaxial loading configuration and any relative angle between the fibres θ 2 , the appearance of the initial debond follows the same tendency of the single-fibre case, even though the alterations of the radial stress along the undamaged fibre–matrix interface could modify the position of the initial debond. However, the experimental evidence presented in [35] (for biaxial T-nT and T-nC cases) reveals that, regardless of the position of the initial debond and the propagation of the interface crack, the latter stages of the mechanism of damage and the final macro-crack are contained in the failure plane perpendicular to the main tensile load. This evidence implies that the main external tension is responsible for the failure and suggests that the location of the initial debond would in any case lead to the interfacial growth assumed here. Therefore, for all the cases considered in this study, the position of the initial debond is the same, in order to obtain a better understanding of the effects of a growing secondary transverse load. The evolution of G calculated with the single-fibre model is represented in Fig. 3. These results obviously coincide with those obtained in [31] and were employed as a reference for the two-fibre cases. As can be observed, for T-0.25 T and T-0.5 T (Fig. 3a), the global G level slightly decreases with respect to the uniaxial case (T-0), whereas for T-0.75 T and T-T (Fig. 3b), it tends to recover the G level of T-0. Regarding the qualitative point of view, propagation under Mode I becomes more relevant as the secondary load increases: the G I level rises above the uniaxial case and the mixed mode propagation ends at larger debonds, which translates the G maxima for T-0.75 T and T-T to higher values of θ d . Thus, Mode I is still active for larger debonds because the secondary transverse load is approximately aligned with the zones where the contact zone was supposed to grow for the uniaxial case [1,6], which produces an inhibition of the contact and pure Mode II propagation as the secondary load increases. However, the G value associated with the first debond considered (Gθd=10◦), which is mostly associated with Mode I propagation, slightly decreases as the secondary load increases, as shown in Fig. 4. This fact could be related to the level of load needed to initiate the propagation, as explained in previous studies [6,31]. Thus, when compared to the uniaxial case, the growth of the incipient crack would require a higher level of external load to be applied. Therefore, increasing values of the Fig. 2. Model including a primary fibre with an interface crack and an undamaged secondary fibre, under two different biaxial configurations. Table 1 Elastic properties of the bi-material system. Material Young modulus, E (Pa) Poisson coefficient, ν Matrix (epoxy) E m =2.79⋅10 9 ν m =0.33 Fibre (glass) E f =7.08⋅10 10 ν f =0.22 C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 4 secondary tensile load involve a slight inhibition of the mechanism of damage, representing a protective effect versus failure with respect to the uniaxial case. Another interesting analysis regarding the effect of the biaxial tensile load versus failure can be carried out if the maximum value of the Energy Release Rate (G max ) is now selected as a representative parameter with respect to the inhibition or promotion of failure, as commented in [41]. In view of the results represented in Fig. 3, and taking the uniaxial case as a reference, the biaxial loading cases tend to inhibit the propagation of the interface crack, as the value of G max for these cases is smaller than the reference case. Thus, the cases T-0 and T-T would be the most prone to failure, whereas the cases T-0.5 T and T-0.75 T would be the less likely to cause crack propagation. Once the single-fibre case was analysed, the same approach was implemented for the two-fibre model. The evolution of G during the propagation (10◦≤θ d ≤150◦) is calculated for the biaxial tensile cases T- nT. The propagation occurs in the presence of the undamaged secondary fibre for a fixed value of the distance between the centres of the fibres (r 2 =r 2 0 ) and a range of relative angles (25◦≤θ 2 ≤155◦). In Fig. 5, these results are plotted for a selection of θ 2 values (θ 2 =25◦, 40◦, 60◦, 90◦, 120◦, 155◦). The G evolution of the interface crack for the angles θ 2 =25◦and θ 2 =155◦(fibres approximately aligned with the primary load) are represented in Fig. 5a and 5f. For these values of θ 2 , the same tendency is observed for increasing values of the secondary tensile load, in particular, the global G level decreases until θ d ≈80◦-100◦and increases for larger debonding angles. Nevertheless, for θ 2 =40◦(Fig. 5b), the biaxial load does not considerably alter the evolution of G, except for the increase produced for large debonding angles as the secondary tension increases. Regarding the intermediate positions, they follow a different tendency. As shown in Fig. 5c, 5d and 5e for θ 2 =60◦, θ 2 =90◦and θ 2 =120◦, the global G level remains almost unaltered for n =0.25, 0.5 with respect to the uniaxial case and increases considerably for n =0.75, 1. It is important to remark the displacement of the position of the maximum value of G for θ 2 =90◦(θ d =100◦), which is related to the fact that the relative position of the fibres is aligned with the secondary load. This tends to amplify the radial displacement of the crack and, consequently, the disappearance of Mode I is reached at greater debonding angles than those of the single-fibre case, as observed in [6,7] (this fact will be further discussed). Another interesting analysis is the comparison between the results of the two-fibre model for each biaxial case T-nT and the corresponding results of the single-fibre model. Considering the angles θ 2 =25◦, 40◦, 60◦, 90◦, 120◦, 155◦, this analysis is shown for the cases n =0, 0.25, 0.5, and n =0.5, 0.75, 1 in Figs. 6 and 7, respectively. In general, in view of the results for a given value of θ 2 , as the secondary load increases, G evolutions for all positions of the secondary fibre tend to approach the G level of its corresponding single-fibre case, which becomes more relevant when n≥0.5 (Fig. 7). This effect means that the biaxial tensile load inhibits the alterations brought about by the secondary fibre, as the differences between all θ 2 cases among themselves and with respect to the single-fibre case become less pronounced for increasing values of the secondary load. In other words, as n increases, the evolution of G for all the θ 2 cases under consideration tends to converge toward the corresponding single-fibre case. When Figs. 6 and 7 are analysed in more detail, it is worth highlighting that for the range 40◦≤θ 2 ≤90◦, (e.g. θ 2 =60◦, Fig. 6c and 7c) and increasing values of the secondary load, the evolution of G tends to match the single-fibre case for large debonding angles (θ d >90◦). On the contrary, for the range 90◦≤θ 2 ≤140◦, the evolution of G matches the single-fibre case for small debonding angles (e.g. θ 2 =120◦, Fig. 6e and 7e). It is also interesting to evaluate the effect of the secondary fibre for the same biaxial configuration. Considering the coefficient n=0.5, Fig. 8 shows the evolution of G, G I and G II , calculated with the two-fibre model for a selection of θ 2 values, also including the single-fibre curves associated with the cases T-0 and T-0.5 T. When the two fibres are approximately aligned with the primary load (25◦≤θ 2 ≤40◦and 140◦≤θ 2 ≤155◦), the evolution of G is totally or partially above the single-fibre case. For the rest of the positions, G reaches lower values than those of the single-fibre case. This issue was also observed and explained in [6] for uniaxial tension, but the differences between the G evolutions with reference to the single-fibre case are less remarkable for the T-0.5 T case. Regarding Mode I propagation, for 25◦≤θ 2 ≤50◦G I evolutions present significant differences with respect to the single-fibre case, both in the level of the curves and the position of their maxima, which tends to occur for θ d ≈θ 2 . This fact is produced by an increase in the radial displacement as the fibres tend to be aligned [6,7]. For 50◦≤θ 2 ≤120◦, G I remains below the single-fibre case evolution, whereas for 130◦≤θ 2 ≤155◦G I stays above. The presence of the secondary fibre also affects the disappearance of G I . For 40◦≤θ 2 ≤60◦, it occurs at larger debonding angles than in the biaxial single-fibre case, with all G I Fig. 3. G/G 0 , G I /G 0 and G II /G 0 versus θ d (single-fibre case), for T-nT: a) n =0, 0.25, 0.5; b) n =0.5, 0.75, 1. Fig. 4. Detail of G/G 0 versus θ d (single-fibre case) for T-nT (n =0, 0.25, 0.5, 0.75, 1). C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 5 evolutions vanishing for higher θ d values than in the uniaxial single-fibre case. With reference to the G II distributions, the maximum is reached at θ d =70◦for the single-fibre case. Nevertheless, for 40◦≤θ 2 ≤70◦, the maximum takes place at a larger θ d value than the single-fibre case and for 80◦≤θ 2 ≤140◦the maximum occurs at a smaller θ d value. The results represented in Fig. 9 for the evolution of G for the case T- T and a selection of θ 2 values confirm the previously commented tendencies. As the value of the secondary load increases (in view of Figs. 8 and 9), and with respect to the uniaxial case [6], the differences between the evolution of G for the two-fibre cases and their corresponding singlefibre case become smaller. For T-T, these differences are more relevant for 60◦≤θ 2 ≤90◦, due to the fact that the disappearance of Mode I is produced for larger debonding angles (radial displacements increase as the fibres tend to be aligned with the secondary load), as well as the G II maximum. The analysis of the G value associated with the first debond considered (Gθd=10◦)was also carried out for the two-fibre and biaxial loading cases. Fig. 10 represents the variation of Gθd=10◦for the two-fibre cases with respect to the single-fibre case. It is observed that, as could be expected from Fig. 4, when the secondary load was taken into consideration, Gθd=10◦is lower than the single-fibre uniaxial case value for most of the positions of the secondary fibre in the T-0.5 T case. For T-T, this fact occurs for any position of the secondary fibre, confirming the inhibition of the mechanism of failure provoked by the biaxial loading, previously commented in this section. The location of the most protective effect produced by the secondary fibre is also modified as n increases: it was found at θ 2 =60◦for T-T, whereas for n =0, 0.5, it was located at θ 2 =40◦. If the results of the two-fibre and single-fibre cases for each T-nT case are compared, taking the corresponding single-fibre case as a reference, Gθd=10◦is above the single-fibre case for a wider range of θ 2 values as the secondary load increases. However, the variations of Gθd=10◦with respect to the corresponding single-fibre case become smaller, confirming the hypothesis that the secondary tensile load reduces the influence of the secondary fibre on the interface crack growth. In order to quantify these effects on failure, in [31] the authors explained the relation between the parameters from Interfacial Fracture Mechanics stress and energetic approaches for a small debond that grows in pure Mode I. More specifically, it was confirmed that the ratio between the uniaxial case and the biaxial cases for the maximum value of the radial stress along the interface ( σ max rr (T−0) σ max rr (T−nT))and for the square root of the G value associated with the first debond considered ( Gθd=10◦(T−0) Gθd=10◦(T−nT) √)were equivalent. Thus,  Gθd=10◦(T−0) Gθd=10◦(T−nT) √can be used as an Fig. 5. G/G 0 versus θ d for T-nT (n =0, 0.25, 0.5, 0.75, 1): a) θ 2 =25◦, b) θ 2 =40◦, c) θ 2 =60◦, d) θ 2 =90◦, e) θ 2 =120◦, f) θ 2 =155◦. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 6 indicator of the increase in the load that has to be applied to initiate the propagation of the interface crack. In the case of the single-fibre model, the results found (see Table 2) confirm a very close agreement with those obtained in [31] for θ d =5◦, for all the T-nT cases considered. This ratio can be modified for the two-fibre cases, which results in an interesting way to quantify the inhibition or promotion of the failure produced by the secondary fibre. Thus, the ratio  Gθd=10◦(T−0,SF) Gθd=10◦(T−nT) √is proposed as an indicator of the effects produced by the presence of the secondary fibre. It is necessary to point out that, the reference value Gθd=10◦corresponds to the uniaxial single-fibre case. Table 2 represents the values of this ratio for all θ 2 values and the T-nT cases. The results highlight that, as the secondary load increases, the differences in the increase of needed load to initiate the propagation become smaller and the accelerative effect (initially found for some θ 2 values) of the presence of the secondary fibre against failure ( Gθd=10◦(T−0,SF) Gθd=10◦(T−nT) √<1)tends to disappear. 3.2. Propagation of the interface crack The end of the unstable crack growth can be predicted as explained in [31], using the estimation of the critical value of the Energy Release Rate, G c , based on [42], and comparing it with the corresponding G evolution: Gc( ψ k) = G1c[1+tan2[(1−λ) ψ k]](3) where G 1c is the critical value of G c for Mode I and λ is the fracture mode sensitivity parameter (in this study, λ=0.25, a typical value for the bimaterial system under consideration). ψ k is the local phase angle and represents the evolution of the fracture mode mixity. When the finite contact zone at the crack tip is developed and Mode I disappears, ψ k= 90◦is reached. ψ k can be calculated following [43]: ψ k(G) = 0.5arccos[F( ε )−1GI−GII GI+GII ](4) where F( ε ) = 1+( π 2 3−2) ε 2+O( ε 4)and ε is the oscillatory index defined as: ε =1 2 π ln 1−β 1+β(5) β being Dundurs parameter [44]. For the bi-material system under consideration ε = − 0.074 and F( ε )−1=1.0071. Fig. 6. Comparison of G/G 0 versus θ d for T-nT (n =0, 0.25, 0.5), with reference to the single-fibre model: a) θ 2 =25◦, b) θ 2 =40◦, c) θ 2 =60◦, d) θ 2 =90◦, e) θ 2 = 120◦, f) θ 2 =155◦. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 7 In the absence of direct experimental data, G 1c has been adjusted as detailed in [31], making the G c value associated with the first debond considered coincide with the initial G value and establishing the uniaxial single-fibre G 1c value as the reference for all the cases under consideration. Thus, all G evolutions start at the same value, and therefore, a different level of load is needed to initiate the propagation for each case. The comparison between G and G c is shown in Fig. 11 for the singlefibre case and the biaxial tensile loads T-nT, with n =0, 0.25, 0.5. The n ≥0.5 cases are not included in Fig. 11 for the sake of clarity and simplicity, these cases following the same tendency shown in the figure. The crack length obtained for the end of the unstable growth is determined from the comparison between G and G c at θ d ≈55◦for T-0, θ d ≈60◦ for T-0.25 T and θ d ≈70◦for T-0.5 T, coinciding with the appearance of a finite size contact zone at the crack tip. Therefore, the results obviously confirm the tendency observed in [31]: for increasing values of the secondary tensile load the interface crack reaches greater lengths of unstable growth. The same approach was applied to the two-fibre cases. As a summary of the results of the prediction of the interface crack growth, a scheme representing the crack length at the end of the unstable growth is plotted in Fig. 12. This figure includes the single-fibre case and a selection of two-fibre cases, for all the biaxial T-nT loads (n =0, 0.25, 0.5, 0.75, 1). For n =0, 0.25 and the range 25◦≤θ 2 ≤50◦, the crack reaches larger lengths of unstable growth than the single-fibre case. For n=0.5, this range of larger lengths is 25◦≤θ 2 ≤60◦. As the level of the secondary load increases, this θ 2 range becomes wider, it being 40◦≤θ 2 ≤100◦for n=0.75 and 50◦≤θ 2 ≤120◦for n=1. On the contrary, the positions where the lengths of unstable growth are smaller than in the single-fibre case are 60◦≤θ 2 ≤90◦for n=0, 60◦≤θ 2 ≤80◦for n=0.25, 80◦≤θ 2 ≤90◦for n=0.5, θ 2 =130◦and 150◦≤θ 2 ≤155◦for n=0.75 and 25◦≤θ 2 ≤30◦and 140◦≤θ 2 ≤155◦(positions aligned with the primary load) for n=1. Moreover, the end of the unstable growth is associated with the appearance of a finite size contact zone at the crack tip (or even before the appearance of the contact zone), except for the cases 30◦≤θ 2 ≤50◦(n =0, 0.25), since the crack stops growing after a considerable length of contact zone has been achieved. It can also be observed in Fig. 12 that, for increasing values of the secondary load, the lengths of the debonding angles for the end of the unstable growth are larger, but there are fewer differences between all θ 2 cases and the single-fibre results. This fact suggests that, as mentioned in Section 3.1 in relation to the G evolutions of the two-fibre and singlefibre cases, when the secondary load increases, the effects produced by the biaxial load tend to mitigate the effects of the presence of the secondary fibre. Fig. 7. Comparison of G/G 0 versus θ d for T-nT (n =0.5, 0.75, 1), with reference to the single-fibre model: a) θ 2 =25◦, b) θ 2 =40◦, c) θ 2 =60◦, d) θ 2 =90◦, e) θ 2 = 120◦, f) θ 2 =155◦. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 8 Fig. 8. G/G 0 , G I /G 0 and G II /G 0 versus θ d for T-0.5 T and θ 2 =25◦, 40◦, 60◦, 90◦, 130◦, 155◦. Fig. 9. G/G 0 , G I /G 0 and G II /G 0 versus θ d for T-T and θ 2 =25◦, 40◦, 60◦, 90◦, 130◦, 155◦. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 9 4. T-nC Cases 4.1. Energy Release Rate The same approach explained in the previous section was employed for the cases T-nC, i.e. tensile load T responsible for the failure and a secondary transverse compressive load, whose value is determined by factor n (n =0, 0.25, 0.5, 0.75, 1). Thus, the Energy Release Rate, G, is calculated from Eq. (1) for a growing interface crack, considering an initial debond θ d =10◦(which appears aligned with the direction of the tensile load) and a debonding interval 10◦≤θ d ≤150◦. The evolution of G for the single-fibre case is represented in Fig. 13. The results obviously agree with those presented in [31]. Contrary to the results represented in Fig. 3 for the T-nT cases, the effects of the compressive secondary load follow the same tendency for all the values of factor n: the values of G increase as the secondary load increases, especially within the range 40◦≤θ d ≤100◦, with no remarkable alterations from a qualitative point of view. This increase in the Energy Release Rate affects the values of Gθd=10◦and G max , selected as representative parameters of the effect of the secondary load on the inhibition or promotion of failure (as mentioned in Section 3.1). Thus, the presence of the secondary compressive load always represents an accelerative effect on the propagation of the interface crack, as expected in view of the results presented in [31]. This increase in the evolution of G is mainly due to the increase in G II , as the compressive secondary load does not alter considerably the level of G I . Moreover, the range of the propagation under mixed mode is slightly shortened by the secondary load: the disappearance of Mode I propagation occurs at θ d =60◦for T-0 and at θ d =50◦for T-C. When including a secondary fibre in the model, the tendency on the G evolution is the same. The level of G increases as the secondary load increases, as can be observed in Fig. 14 for all the relative positions between the fibres represented (θ 2 =25◦, 40◦, 60◦, 90◦, 120◦, 155◦). This is especially relevant for θ 2 =25◦, 155◦(fibres aligned with the tensile load), θ 2 =90◦(fibres aligned with the compressive load) and θ 2 =60◦ (intermediate position of the secondary fibre, in front of the crack growth range). The secondary load affects to a lesser extent other positions like θ 2 =40◦or θ 2 =120◦(intermediate position, further from the Fig. 10. Variation of G θd=10◦for the two-fibre cases with respect to the singlefibre case (T-nT, n =0, 0.5, 1). Table 2  Gθd=10◦(T−0,SF) Gθd=10◦(T−nT) √for the T-nT cases and the single-fibre (SF) and two-fibre (TF) models. Case SF TF θ 2 =25◦TF θ 2 =30◦TF θ 2 =40◦ T-0 1.000 0.975 1.307 1.770 T-0.25 T 1.024 1.009 1.254 1.609 T-0.5 T 1.050 1.042 1.201 1.472 T-0.75 T 1.077 1.079 1.155 1.356 T-T 1.105 1.121 1.113 1.254 Case TF θ 2 =50◦TF θ 2 =60◦TF θ 2 =70◦TF θ 2 =80◦ T-0 1.483 1.226 1.125 1.106 T-0.25 T 1.456 1.259 1.168 1.138 T-0.5 T 1.426 1.295 1.211 1.173 T-0.75 T 1.394 1.327 1.254 1.206 T-T 1.358 1.358 1.301 1.237 Case TF θ 2 =90◦TF θ 2 =100◦TF θ 2 =110◦TF θ 2 =120◦ T-0 1.121 1.127 1.109 1.062 T-0.25 T 1.133 1.129 1.106 1.065 T-0.5 T 1.150 1.131 1.102 1.067 T-0.75 T 1.164 1.131 1.098 1.068 T-T 1.177 1.132 1.094 1.070 Case TF θ 2 =130◦TF θ 2 =140◦TF θ 2 =150◦TF θ 2 =155◦ T-0 1.003 0.943 0.892 0.875 T-0.25 T 1.015 0.968 0.927 0.911 T-0.5 T 1.027 0.992 0.962 0.952 T-0.75 T 1.039 1.015 1.000 0.997 T-T 1.052 1.042 1.039 1.046 Fig. 11. G/G 0 and G c /G 0 versus θ d for the single-fibre case and T-nT (n =0, 0.25, 0.5). Fig. 12. Scheme of the crack length at the end of the unstable growth (T-nT). C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 16 and θ 2 =115◦, 157◦for Fig. 25b). The configuration presented in Fig. 25a could be compared with the prediction of the two-fibre model for the end of the unstable growth for the case T-0.75C. This prediction gives the values θend d(θ2=155◦)≈50◦and θend d(θ2=140◦)≈45◦, which corresponds to a total crack length of 95◦. This prediction is relatively close to the microscopical observation (total crack length of 80◦). The difference in the crack length and its lack of symmetry could be associated with the influence of the rest of the fibres and the slight difference between both loading cases. For the uniaxial case represented in Fig. 25b, the prediction of the two-fibre model gives θend d(θ2=155◦)≈50◦and θend d(θ2=110◦,120◦)≈40◦, and a total length of 90◦, in full agreement with the total length of the microscopical observation (89◦). Considering all the information provided by Figs. 22-24, the influence of the secondary fibre is more important for the T-nC cases and increasing values of n. This fact would imply that, experimentally, if each biaxial configuration is applied progressively, and observing for a similar applied load, interface cracks should appear before and they could be more frequently observed for T-nC and the range 150◦≤θ 2 ≤155◦, even though the crack lengths would be considerably shorter than in the uniaxial case. Interface cracks should be rarely observed also for θ 2 =40◦and T-C. In general, and for the same applied load, for T-nT and increasing values of n, interface cracks observation should be less frequent than for T-nC, even though the extent of cracks Fig. 23. G max versus θ 2 and all the biaxial configurations. Fig. 24. θ end d versus θ 2 and all the biaxial configurations. Fig. 25. Microscopical observations of interface cracks at cruciform specimens subjected to biaxial loading a) T-0.7C and b) T-0. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 17 should be greater. These predictions cannot be easily checked in terms of the comparison between the T-nT and the T-nC cases, as the experimental results reveal that the external load at the moment of rupture is lower for the T-nC cases, and the failure is very explosive [35]. However, this fact agrees with the increase in the accelerative effects found for the T-nC cases and the generalised protective effect of the T-nT cases. 6. Conclusions In this work, the propagation of an interface crack subjected to biaxial loading has been studied by means of a two-fibre BEM model. The objective was to analyse the effect of the presence of a nearby fibre on the interface crack growth under different loading configurations. The biaxial configuration consists in a main tensile load, nominally responsible for the failure and a secondary tensile or compressive load, perpendicular to the main load. The results have been analysed under the light of Interfacial Fracture Mechanics. Thus, this work has been focused on the study of the Energy Release Rate for different configurations of the two-fibre model, in order to compare the results with previous uniaxial and biaxial single-fibre studies. This comparison is based on the influence of three main parameters: the type of secondary load (tensile or compressive), the level of load and the relative position between the fibres. Some of the most important conclusions are summarised next. With regard to the effects of the biaxial loading and the secondary fibre on the propagation of the interface crack, some effects were observed as the secondary load increased: •With reference to the uniaxial case, when considering a secondary tensile load and its value increased, the end of the unstable growth was produced at greater debonding angles for a fixed value of θ 2 . Nevertheless, a secondary compressive load implied, for most of the cases, shorter crack lengths at the end of the unstable growth. These effects were already detected for a single-fibre model in [31], but the intensity of the effects depends on the position of the secondary fibre. •For the T-nT biaxial configuration, as the secondary load increased the differences between the θ 2 cases and their corresponding singlefibre case became smaller. On the contrary, for the T-nC cases, these differences were more accentuated as the secondary load increased. •For increasing values of the transverse tensile load, the θ 2 ranges for which the length of unstable growth was greater than its corresponding single-fibre case were wider. However, as the secondary compressive load increased, this fact only occurred for a few relative positions between the fibres. Additionally, micromechanical experimental evidence has been presented supporting the aforementioned conclusions. According to the study of the Energy Release Rate for a given biaxial configuration and a relative position between the fibres, the maximum and initial values (G max and Gθd=10◦) were used as representative parameters of the inhibitions or promotion of failure, with respect to the uniaxial single-fibre case. In this regard, when a secondary tensile load was included (T-nT), for most of the relative positions between the fibres the presence of the nearby fibre intensified the inhibition and reduced the promotion of the mechanism of failure observed in the uniaxial tensile case [6]. For T-nC, the biaxial loads implied the opposite effect, intensifying the most protective or accelerative effects found in [6]. On a final note, considering the combination of both effects, biaxial loading and the surrounding fibre, on the aspects analysed in this study (with reference to the uniaxial single-fibre case), as the level of the secondary tensile load increases, the biaxial loading tends to inhibit the presence of the undamaged secondary fibre. On the contrary, when a secondary compressive load is considered, as its value increases, all the alterations induced by the secondary fibre are intensified. CRediT authorship contribution statement Carlos Sandino: Methodology, Validation, Investigation, Writing – review & editing. Elena Correa: Conceptualization, Methodology, Validation, Writing – review & editing, Supervision, Funding acquisition. Federico París: Conceptualization, Writing – review & editing, Supervision, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgement The authors would like to thank Prof. E. Graciani whose BEM code was used, and Dr P.L. Zumaquero for the help in the search of experimental evidence. This work was supported by: •Consejería de Economía y Conocimiento Junta de Andalucía. Project P18-FR-3360 corresponding to the I +D +I Programa Operativo FEDER 2014-2020 and Junta de Andalucía and European Social Fund. •Spanish Ministry of Universities: ’Margarita Salas’ grants for the training of young doctors. Funded by the Spanish Recovery, Transformation and Resilience Plan and Next Generation EU. •Spanish Ministry of Science and Innovation. Project PID2021- 126279OB-I00 corresponding to the State Plan 2021-2023: Knowledge Generation Projects. References [1] F. Parıs, E. Correa, J. Ca˜ nas, Micromechanical view of failure of the matrix in fibrous composite materials, Compos. Sci. Technol. 63 (7) (2003) 1041–1052. [2] F. París, E. Correa, V. Mantiˇ c, Kinking of transversal interface cracks between fibre and matrix, J. App. Mech. 74 (4) (2007) 703–716. [3] E. Correa, V. Mantiˇ c, F. París, Numerical characterisation of the fibre-matrix interface crack growth in composites under transverse compression, Eng. Fract. Mech 75 (14) (2008) 4085–4103. [4] E. Correa, V. Mantiˇ c, F. París, A micromechanical view of inter-fibre failure of composite materials under compression transverse to the fibres, Compos. Sci. Technol. 68 (9) (2008) 2010–2021. [5] E. Correa, F. París, V. Mantiˇ c, Effect of thermal residual stresses on the matrix failure under transverse compression at micromechanical level – a numerical and experimental study, Compos. Part A- Appl. S 43 (1) (2012) 87–94. [6] C. Sandino, E. Correa, F. París, Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load, Eng. Fract. Mech 168 (2016) 58–75. [7] C. Sandino, E. Correa, F. París, A study of the influence of a nearby fibre on the interphase crack growth under transverse compression in composite materials, Eng. Fract. Mech. 193 (2018) 1–16. [8] V.I. Kushch, S.V. Shmegera, P. Brøndsted, L. Mishnaevsky, Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading, Int. J. Engng. Sci. 49 (1) (2011) 17–29. [9] L. T´ avara, V. Mantiˇ c, E. Graciani, F. París, Modelling interfacial debonds in unidirectional fibre-reinforced composites under biaxial transverse loads, Compos. Struct. 136 (2016) 305–312. [10] V. Mantiˇ c, Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. application of a coupled stress and energy Criterion, Int. J. Solids Struct. 46 (6) (2009) 1287–1304. [11] P.A. Carraro, M. Quaresimin, Modelling fibre–matrix debonding under biaxial loading, Compos. A Appl. Sci. Manuf. 61 (2014) 33–42. [12] M. Hojo, M. Mizuno, T. Hobbiebrunken, T. Adachi, M. Tanaka, S.K. Ha, Effect of fiber array irregularities on microscopic interfacial normal stress states of transversely loaded UD-CFRP from viewpoint of failure initiation, Compos. Sci. Technol. 69 (11–12) (2009) 1726–1734. C. Sandino et al. Theoretical and Applied Fracture Mechanics 127 (2023) 103956 18 [13] F. Sharifpour, J. Montesano, R. Talreja, Assessing the effects of ply constraints on local stress states in cross-ply laminates containing manufacturing induced defects, Compos. B Eng. 199 (2020). [14] L. Zhuang, R. Talreja, J. Varna, Transverse crack formation in unidirectional composites by linking of fibre/matrix debond cracks, Compos. A 107 (2018) 294–303. [15] L. Zhuang, A. Pupurs, J. Varna, R. Talreja, Z. Ayadi, Effects of inter-fiber spacing on fiber-matrix debond crack growth in unidirectional composites under transverse loading, Compos. A 109 (2018) 463–471. [16] L. Di Stasio, J. Varna, Z. Ayadi, Energy release rate of the fiber/matrix interface crack in UD composites under transverse loading: effect of the fiber volume fraction and of the distance to the free surface and to non-adjacent debonds, Theor. Appl. Fract. Mech. 103 (2019) 102251. [17] L. Di Stasio, J. Varna, Z. Ayadi, Effect of the proximity to the 0◦/90◦interface on Energy Release Rate of fiber/matrix interface crack growth in the 90◦-ply of a cross-ply laminate under tensile loading, J. Compos. Mater. 54 (21) (2020) 3021–3034. [18] H. Saito, H. Takeuchi, I. Kimpara, A study of crack suppression mechanism of thinply carbon-fiber-reinforced polymer laminate with mesoscopic numerical simulation, J. Compos. Mater. 48 (17) (2014) 2085–2096. [19] Arteiro A, Catalanotti G, Melro AR, Linde P, Camanho PP. Micro-mechanical analysis of the effect of ply thickness on the transverse compressive strength of polymer composites. Compos Part A-Appl Sci Manuf 2015; 79: 127–37. [20] París F, Velasco ML, Correa E. The scale effect in composites: An explanation physically based on the different mechanisms of damage involved in failure. Compos Struct 2021;257:113089. [21] H. Zhang, M.L. Ericson, J. Varna, L.A. Berglund, Transverse single-fibre test for interfacial debonding in composites: 1. experimental observation, Compos. Part A- Appl. Sci. Manuf. 28 (4) (1997) 309–315. [22] E.K. Gamstedt, B. Sj¨ ogren, Micromechanics in tension-compression fatigue of composite plies containing transverse plies, Compos. Sci. Technol. 59 (1999) 167–178. [23] M. Gentz, D. Armentrout, P. Rupnowski, L. Kumosa, E. Shin, J.K. Sutter, M. Kumosa, In-plane shear testing of medium and high modulus woven graphite fiber reinforced/polyimide composites, Compos. Sci. Technol. 64 (2) (2004) 203–220. [24] N. Baral, H. Guezenoc, P. Davies, C. Baley, High modulus carbon fibre composites: correlation between transverse tensile and Mode I interlaminar fracture properties, Mater. Lett. 62 (6-7) (2008) 1096–1099. [25] T. Hobbiebrunken, M. Hojo, T. Adachi, C. De Jong, B. Fiedler, Evaluation of interfacial strength in CF/epoxies using FEM and in-situ experiments, Compos. A 37 (12) (2006) 2248–2256. [26] H. Saito, H. Takeuchi, I. Kimpara, Experimental evaluation of the damage growth restraining in 90◦layers of thin-ply CFRP cross-ply laminates, Adv. Compos. Mater 21 (1) (2012) 57–66. [27] K. Martinyuk, B.F. Sørensen, P. Modregger, E.M. Lauridsen, 3D in situ observations pf glass fibre/matrix interfacial debonding, Compos. Part A- Appl. Sci. Manuf. 55 (2013) 63–73. [28] E. Correa, M.I. Valverde, M.L. Velasco, F. París, Microscopical observations of inter-fibre failure under tension, Compos. Sci. Technol. 155 (2018) 213–220. [29] P.L. Zumaquero, E. Correa, J. Justo, F. París, Microscopical observations of interface cracks from inter-fibre failure under compression in composites laminates, Compos Part A-Appl Sci Manuf 110 (2018) 76–83. [30] Z. Hashin, A. Rotem, A fatigue failure criterion for fibre reinforced materials, J. Comp. Mat. 7 (1973) 448–464. [31] E. Correa, F. París, V. Mantiˇ c, Effect of the presence of a secondary transverse load on the inter-fibre failure under tension, Eng. Fract. Mech. 103 (2013) 174–189. [32] E. Correa, F. París, V. Mantiˇ c, Effect of a secondary transverse load on the interfibre failure under compression, Compos. B 65 (2014) 57–68. [33] E. Correa, A. Barroso, M.D. P´ erez, F. París, Design for a cruciform coupon used for tensile biaxial transverse test son composite materials, Compos. Sci. Technol. 145 (2017) 138–148. [34] A. Barroso, E. Correa, J. Freire, F. París, A device for biaxial testing in uniaxial machines. design, manufacturing and experimental results using cruciform specimens of composite materials, Exp. Mech. 58 (1) (2018) 49–53. [35] P. Lucía Zumaquero, E. Correa, J. Justo, F. París, Transverse biaxial tests on long fibre reinforced composites, Comput. Struct. 297 (2022). [36] París F, Ca˜ nas J. Boundary element method. Fundamentals and applications. Oxford: OUP; 1997. [37] E. Graciani, V. Mantiˇ c, F. París, A. Bl´ azquez, Weak formulation of axi-symmetric frictionless contact problems with boundary elements, Application to interface cracks. Comput. Struct. 83 (10-11) (2005) 836–855. [38] G.R. Irwin, Analysis of stresses and strain near the end of a crack transversing a plate, J. App. Mech. 24 (1957) 361–364. [39] M. Toya, A crack along the interface of a circular inclusion embedded in an infinite solid, J. Mech. Phys. Solids 22 (5) (1974) 325–348. [40] Murakami Y. Stress intensity factor handbook. Oxford: Pergamon Press; 1988. [41] M.L. Velasco, E. Correa, F. París, Interaction between fibres in the transverse damage in composites, Eng. Fract. Mech. 239 (2020). [42] J.W. Hutchinson, Z. Suo, Mixed mode cracking in layered materials, Adv. App. Mech. 29 (1992) 63–191. [43] V. Mantiˇ c, F. Par´ Is, Relation between SIF and ERR based measures of fracture mode mixity in interface cracks, Int. J. Fracture 130 (2) (2004) 557–569. [44] J. Dundurs, Discussion of a paper by DB Bogy, J Appl Mech 36 (1968) 650–652. C. Sandino et al.