# Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination

H. Girard, A. Doitrand, B. Koohbor, R.G. Rinaldi, N. Godin, D. Long, J. Bikard

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

H. Girard, A. Doitrand, B. Koohbor, R.G. Rinaldi, N. Godin, et al., Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination. Journal of the Mechanics and Physics of Solids, 2023, pp.105498. 10.1016/j.jmps.2023.105498. hal-04309201

HAL Id: hal-04309201

https://hal.science/hal-04309201v1

Submitted on 27 Nov 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination

H. Girard<sup>a,b,*</sup>, A. Doitrand<sup>a</sup>, B. Koohbor<sup>c</sup>, R.G. Rinaldi<sup>a</sup>, N. Godin<sup>a</sup>, D. Long<sup>a</sup>, J. Bikard<sup>b</sup>

<sup>a</sup>Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France

<sup>b</sup>Axel'One (Solvay), 87 Avenue des Frères Perret, CS 70061, 69192 Saint Fons, France

<sup>c</sup>Department of Mechanical Engineering, Rowan University, 201 Mullica Hill Rd., Glassboro, NJ 08028, United States of America

## Abstract

Fiber-matrix interface debonding in two-fiber specimens under remote tensile loading is studied both experimentally and numerically by means of a coupled stress and energy criterion. Depending on its relative position, the neighboring fiber induces a perturbation of both stress and energy fields at the reference fiber interface which results in asymmetrical debonding initiation and propagation. The determination of the debonding initiation and propagation shape is addressed based on either i) stress isocontours, ii) energy isocontours or iii) the Coupled Criterion (CC). It was found that the debonding initiation configuration can be determined based on stress (respectively energy) isocontours for small (respectively large) enough interface brittleness number. For intermediate brittleness number, the debonding initiation configuration cannot be obtained using neither the stress nor the energy isocontours, but requires a coupling of both aspects. Despite different initiation debonding configurations, the corresponding initiation remote stresses do not differ much, which results in similar debonding configurations after unstable crack propagation following initiation.

**Keywords:** Fiber-matrix debonding, Finite Fracture Mechanics, Coupled Criterion, Linear Elastic Fracture Mechanics, Neighboring fiber

## 1. Introduction

Interface fracture is a challenging problem to address from both experimental and numerical points of view. First, the very definition of an interface may not be straightforward. In numerical models, it is generally defined either as a zero-thickness surface that links two different materials or possibly considering a finite thickness when its dimension is comparable or not negligible with respect to the other material dimensions. Depending on the materials involved and the manufacturing process, the transition from one material to another might not be just a surface separating two homogeneous and uniform materials, but a transition zone of a finite thickness may also be present [1, 2, 3]. Then, interface fracture is another challenge as it is generally described in models as the separation of the two lips of the surface between both materials, including or not a process

<footnote>*Corresponding author</footnote>

<footnote>Email address: hugo.girard@insa-lyon.fr (H. Girard)</footnote>
zone. In experiments, either adhesive, cohesive, or a combined adhesive-cohesive cracking may be encountered depending on the interface bonds and homogeneity (possible presence of local defects for instance). As a consequence, the viewpoint of interface fracture modeling can be regarded as a homogenized vision that gathers all these mechanisms. The direct consequence is that interface fracture properties (e.g., strength, fracture toughness) may not be straightforward. Nevertheless, numerical interface fracture models such as, e.g., Cohesive Zone Models (CZM) [4, 5] or the Coupled Criterion (CC) [6] are useful tools to describe material behavior observed experimentally. For instance, the CC enables assessing debonding initiation at the fiber-matrix interface. This approach was introduced by Leguillon [6] to extend Linear Elastic Fracture Mechanics (LEFM) [7] to crack initiation.

Mantič [8] first applied the CC to assess debonding at the interface of a circular inclusion embedded within an elastic matrix subjected to a far field tensile loading. The Hutchinson and Suo [9] relation was used to describe the interface mode mixity leading to the initiation of debonding triggered either by both stress and energy conditions or by the energy condition exclusively. García *et al.* [10] extended the embedded inclusion problem to symmetric debonding at the two opposite poles and showed that non-symmetric initiation was more favorable than symmetric crack initiation due to a larger energy released per crack surface unit. Several extensions of this approach were then proposed for bi-axial loading, for instance, by Mantič and García [11] and Muñoz-Reja *et al.* [12] as well as comparison with other numerical approaches such as Cohesive Zone Model (CZM) [13, 14, 15]. Besides, the CC can be an efficient method to perform an inverse identification of the interface properties, as demonstrated by Doitrand *et al.* [16]. A range of interface shear and opening critical energy release rates and strengths were also determined by Girard *et al.* [17] for a sample with a single glass fiber embedded into an epoxy matrix.

Previous studies have mainly focused on a single fiber in an infinite matrix. However, a more realistic composite microstructure involves numerous fibers with varying relative positions (which can be for instance quantified by inter-center distances and angles with respect to the loading direction). Several authors have investigated the influence of neighboring fibers on the debonding process. Muñoz-Reja *et al.* [12, 18] assessed debonding initiation in a two-fiber configuration under tensile loading, varying the inter-center distance between 2.25 and 2.5 times the fiber radius and compared two different inter-center angles. Asymmetric debonding locations with respect to the loading direction were obtained. It was shown that the smaller the inter-center distance, the smaller the remote tensile load required for the debonding to initiate. Sandino *et al.* [19] studied the influence of angle and distance between fibers center on debonding propagation. Inter-center angle influenced the debonding arrest angle and remote loading required for propagation. Moreover, the influence of the nearby fiber vanishes after a certain distance. Within the same framework, Velasco *et al.* [20] and Zhuang *et al.* [21] proposed a similar analysis in the presence of a pre-existing debonding and found that the debonding of another fiber may either be delayed or favored by the pre-existing debonding depending on its location. It was also found that the debonding kink into the matrix was affected by the neighboring fiber and was either damped or reinforced depending on whether there was a pre-existing debonding. Kushch *et al.* [22] used CZM for the prediction of debonding initiation with a two-fiber configuration and focused on interface stresses and debonding locations that vary with the inter-center angle, either at the top or bottom poles. However, most of the aforementioned studies do not focus on a configuration that is asymmetric with respect to the loading direction, so that the debonding location and propagation remain symmetric. Besides, there is a lack of experimental confrontation to validate the numerical
approach employed.

The asymmetric debonding shape leads to some difficulties in the implementation of the CC which requires the a priori definition of the crack shape. Such shape of the debonding can be prescribed symmetric in the case of a single inclusion [8, 23]. However, multi-fiber configurations induce a perturbation of the fields at the interface and the debonding no longer remains symmetric with respect to the loading direction [18, 22]. To overcome this shape identification problem, some authors estimated the crack shape based on stress isocontours [24, 25, 26]. This solution led to an efficient way of determining the crack shape. However, it does not take into account any energy aspects which, in some cases, can drive the initiation of debonding [8, 17, 27, 28]. In this regard, García *et al.* [29] and Doitrand *et al.* [27] proposed an iterative approach by comparing the CC solution based on a parameterized crack geometry in the vicinity of the stress concentration which, therefore, considers the energy aspects. The two above-mentioned approaches for crack shape determination do not consider all possible crack geometries. As a consequence, it is yet not established which initiation crack shape among all possible crack shapes is the most favorable, i.e. results in the minimum imposed loading at initiation.

The present study aims at determining the CC solution in the case of an asymmetric debonding by simultaneously considering all possible crack shape configurations. The configuration under investigation and related experiments are presented in Section 2. The remote loading required for a debonding to initiate is therefore evaluated with nearby two-fiber specimens having different inter-center angles and distances, as detailed in Section 3. In the following, the CC solution is compared to the stress- and energy-based debonding shape solutions to quantify the accuracy of the stress isocontours assumption in terms of debonding angles and remote loading. The CC solutions are then confronted to experimental observations [30] in Section 4. In details, the location and size of the debonding are compared, as well as the remote loading required for the debonding to initiate.

## 2. Experiments

### 2.1. *Sample geometries and testing protocol*

Two-fiber specimens under uniaxial tensile loading are considered to study the influence of a nearby fiber on fiber-matrix debonding initiation and propagation. Experimental results are taken from [30, 31]. Figure 1 describes specimen geometry and dimensions, sized accordingly to ASTM D638 with a sample thickness of 6 mm. Thirty samples are subjected to uniaxial tensile loading along the (*Oy*) direction, leading to a far-field stress (σ∞). The inter-center angle (α) from the (*Ox*) direction and distance (d) are varied in the range 0, 30, 45, 60 and 90 degrees and 2.5r, 4r and 5r, respectively, with r the fiber radius equal to 1 mm. It is worth noting that both fiber absolute positions do change between the configurations as they are rotated with respect to the inter-center midpoint. Samples are manufactured by embedding two glass macro-fibers in a room temperature curable thermoset epoxy resin whose isotropic elastic properties are detailed in Table 1.

In the following, the fibers are referred to as respectively reference and neighboring fibers and all quantities will be referring to the reference fiber in the sequel (see Figure 1). The angular position at the interface of the reference fiber is denoted θf. In order to study the field perturbation near the fiber and the debonding process, a speckle pattern was applied on the samples surfaces to enable performing 2D Digital Image Correlation (DIC). A camera equipped with a high magnification
<table><caption>Table 1: Linear elastic properties of the epoxy matrix and glass fiber. Fiber properties are provided by the manufacturer and matrix properties are based on 3 in-house measurements from [30].</caption><thead><tr><th>Property</th><th>Epoxy</th><th>Glass Fiber</th></tr></thead><tbody><tr><td>Young Modulus [GPa]</td><td>2.36 ± 0.10</td><td>63</td></tr><tr><td>Poisson's ratio [-]</td><td>0.40 ± 0.04</td><td>0.30</td></tr></tbody></table>

lens was used to capture debonding initiation and propagation in situ. The camera recording was synchronized with the load cell for a similar extraction time.

<figure><img src="image_3.png" /><figcaption>Figure 1: Specimen geometry and dimensions including two fibers (radius $r$), where $d$ is the inter-center distance and $\alpha$ is the angle between the ($Ox$) axis and the direction passing through both fiber centers. The reference fiber interface is divided into a top half circle (T) and a bottom half circle (B) where $\theta_f$ denotes the angular position.</figcaption></figure>

## 2.2. Experimental measurements

The numerical simulation of fiber-matrix debonding initiation and propagation presented in the following enables determining the following physical quantities:

* Debonding initiation location (top half circle, bottom half circle or both fiber half circles);

* Range of debonding initiation angles along the interface ($\theta_d$);

* Remote applied stress ($\sigma^\infty$) associated with debonding initiation.
On the one hand, the extraction of the debonding locations and angles is based on DIC residuals [32]. The debonding leads to the appearance of a new black area at the fiber-matrix interface. Since the emerging black pixel area was not present in the reference image, high correlation residuals can be detected and thus related to the debonding location and size. On the other hand, the determination of the debonding initiation loading is based on previous work by Livingston et al. [30] and Girard et al. [31] where virtual gauges are placed across the fiber-matrix interface for local displacement and strain measurements. The DIC solution is used to extract the gauge length variation as a function of the applied loading. Debonding initiation thus results in a sudden change in slope which allows the determination of the far-field load (and stress by simply dividing the load by the pristine transverse cross section area). All the experimental results obtained for each configuration are detailed in Section 4.

### 3. Debonding shape prediction

#### 3.1. Coupled Criterion implementation

The initiation of debonding is predicted using the CC. A 2D plane strain Finite Element (FE) model is implemented under the assumption of linear elasticity and small deformations. As a consequence, the stress is proportional to the applied load whereas the elastic strain energy is proportional to the square applied load. In the sequel, the calculations were performed with a $\sigma^{\infty} = 6.85$ MPa remote stress. The FE model consists of two 2 mm diameter fibers embedded in the matrix. The material properties presented in Table 1 are assigned to each component. A 0.02 mm mesh size is chosen at the fiber-matrix interface, which ensures that the influence on the stress and elastic strain energy is smaller than 1% for a finer mesh. The geometry of the model is limited to 50 mm and 13 mm along the longitudinal (Oy) and transverse (Ox) directions with limited nodal displacements so that the configuration is close to the experimental setup.

Debonding initiation is assessed using the CC, which combines two conditions to provide the initiation loading and debonding angle range. On one hand, the stress condition must be satisfied at any location over the entire debonding path $\Gamma$ before initiation. As a consequence, the stress requirement is not a point stress but a nonlocal condition since it must be fulfilled all over the debonding path before initiation. It can thus be employed even in the case of a nonuniform stress field. On the other hand, the Incremental Energy Release Rate (IERR) must be larger than the interface average critical Energy Release Rate (ERR). These two conditions are summarized in Equation (1).

$$ \left\{ \begin{array}{l} \sigma_{\text{eq}} = \sqrt{\sigma_{\text{nn}}^2 + \frac{1}{\mu^2} \tau_{\text{nt}}^2} \geq \sigma_c \quad \forall \vec{x} \in \Gamma, \\ G_{\text{inc}}(\theta_d, \sigma^{\infty}) \geq \bar{G}_c(\theta_d). \end{array} \right. \qquad (1) $$

The parameter $\mu$ denotes the $\tau_c$ to $\sigma_c$ ratio, which are respectively the shear and tensile strengths. The IERR is denoted by $G_{\text{inc}}$ and $\bar{G}_c$ is the average critical ERR. The normal and shear stresses of the interface are $\sigma_{\text{nn}}$ and $\tau_{\text{nt}}$, respectively. The calculation of all these quantities is detailed in the sequel. Contrary to the CC application to single-fiber configuration where the initiation debonding location is a priori known and symmetric with respect to (Oy) direction [8, 10, 17], the initiation debonding configuration may now depend on the nearby fiber (intercenter distance and angle, see Figure 1). As a consequence, solving the CC requires to determine not only the remote imposed stress causing debonding initiation but also the initiation debonding configuration, which is addressed in next sections.
## 3.2. Methodology

Previous studies have shown that for multi-fiber samples, the interface stress fields are influenced by their respective neighboring fibers [19, 20, 22]. They affect both the stress level at the interface and the elastic strain energy released during debonding initiation. Consequently, the debonding trajectory can no longer be considered symmetric with respect to the loading configuration and must now be determined.

Based on the CC, debonding can initiate provided the stress and energy criteria given in Equation 1 are fulfilled. As a consequence, it is likely that the potential debonding configurations do not lie far from the ones maximizing either the equivalent stress or the elastic strain energy variation.

Figure 2 shows the variation of the equivalent stress ($\sigma_{eq}$) and elastic strain energy variation ($-\Delta W$) computed as a function of the angular position along the fiber-matrix interface ($\theta_f$). The equivalent stress variation as a function of $\theta_f$ is obtained based on a calculation without debonding. The evaluation of $-\Delta W$ is performed by releasing one interface node at a given angular position $\theta_f$. A constant mesh size along the interface is used to ensure that the angular location maximizing $-\Delta W$ also maximizes the IERR. For a single-fiber configuration under remote tensile loading, the equivalent stress and $-\Delta W$ profiles are symmetric and both their maxima are attained for the same angle (90 deg. angular position, see Figure 2a). Asymmetry in the equivalent stress and elastic strain energy variation nevertheless occurs for a two-fiber configuration, as shown in Figure 2b in a case where $\alpha = 45$ deg. and $d = 2.5r$. More precisely, asymmetric fields prevail and both quantities reach a maximum at a similar 100 deg. angular position, different from the single fiber configuration. In the following $\theta_i$ is used to define the angular position $\theta_f$ maximizing the stress or the elastic strain energy variation, where debonding could most likely initiate from. Table 2 lists the angles maximizing the strain energy released ($\theta_i^{-\Delta W}$) or the equivalent stress criterion ($\theta_i^{\sigma_{eq}/\sigma_c}$) for different inter-center angles at a fixed inter-center distance of 2.5r, distance for which the fields are most affected. The ratio between the two angles are presented in Table 2. All ratios remain close to 1 whatever the configuration meaning that both angles are close to each other, even for highly perturbed configurations.

<figure><img src="image_4.png" /><figcaption>Figure 2: Equivalent stress and elastic strain energy variation as a function of the angular position along the fiber-matrix interface for (a) single-fiber configuration (b) two-fiber configuration for a constant $\mu = 1.5$.</figcaption></figure>
<table><caption>Table 2: $θ_i$ variation for various inter-center angles where $d = 2.5r$ and $μ = 1.5$, evaluated at the reference fiber top half circle. The upperscripts refer to the energy or stress conditions used for the angular position determination. Ratio between the two angles are added to highlight the small discrepancy between the two approaches.</caption><thead><tr><th>α [deg.]</th><th>0</th><th>30</th><th>45</th><th>60</th><th>90</th></tr></thead><tbody><tr><th>θ<sub>i</sub><sup>σ<sub>eq</sub>/σ<sub>c</sub></sup> [deg.]</th><td>87</td><td>93</td><td>98</td><td>63</td><td>90</td></tr><tr><th>θ<sub>i</sub><sup>-ΔW</sup> [deg.]</th><td>89</td><td>95</td><td>102</td><td>60</td><td>90</td></tr><tr><th>θ<sub>i</sub><sup>-ΔW</sup>/θ<sub>i</sub><sup>σ<sub>eq</sub>/σ<sub>c</sub></sup> [-]</th><td>1.02</td><td>1.00</td><td>1.03</td><td>0.96</td><td>1.02</td></tr></tbody></table>

Figure 3 illustrates the field perturbation by pointing out the equivalent stress maximum location $|θ_i|$ for several fiber inter-center distances and angles. This angle is calculated at both top and bottom half circles of the reference fiber. The larger the difference from 90 degrees, the higher the perturbation of the fields with respect to the single fiber configuration. Noticeably, the stress fields at the fiber bottom half circle are only slightly influenced since the angles remain within the 90 deg. vicinity (+ or - 3%) whatever the inter-center distance. On the contrary, at the top half circle of the fiber, the smaller the inter-center distance, the larger the perturbation. Additionally, $θ_i$ differs the most from 90 degrees for an inter-center angle close to 60 degrees. This configuration may thus result in strongly asymmetric crack initiation. In contrast, $θ_i$ reverts to 90 deg. when $α = 90$ deg. as the direction through the two fiber centers is parallel to the load direction.

<figure><img src="image_5.png" /><figcaption>Figure 3: Angle $θ_i$ for which the maximum equivalent stress is attained at the top and bottom fiber half circles (T and B respectively) as a function of the fiber inter-center distance and angle for $μ = ∞$.</figcaption></figure>

Starting from an initiation debonding angular position $θ_i$, the debonding shape can then be described by increments of angles on each side of this initial location, denoted by $Δθ^+$ and $Δθ^-$, as detailed in Figure 4. The sum of these two debonding increments corresponds to the total debonding angle $θ_d$. In the sequel, we evaluate different approaches to determine the possible initiation debonding configurations, based on stress, energy or both stress and energy.

### 3.3. Debonding configuration based on stress criterion

A convenient crack shape determination approach consists in defining the possible crack configurations based on the stress isocontours [24, 25, 26]. The stress isocontour-based crack definition consists in considering areas which strictly fulfill the stress criterion for a given imposed loading and specified interface strengths. As a consequence, it reverts to determining the debonding configuration described by $θ_i$, $Δθ^+$ and $Δθ^-$ for which the stress criterion is fulfilled before initiation
<figure><img src="image_6.png" /><figcaption>Figure 4: Schematic of the debonding shape including the initiation debonding angular position $\theta_i$ and debonding increments $\Delta\theta^+$ and $\Delta\theta^-$ on each side of the initiation location at the reference fiber top interface. The sum of both debonding increments correspond to the total debonding angle $\theta_d$.</figcaption></figure>

(see Figure 4). The two debonding tips respectively correspond to angles $\theta_i + \Delta\theta^+$ and $\theta_i - \Delta\theta^-$. The advantage of this approach is that, for given shear and tensile strengths, the possible crack surfaces can be determined based on a single calculation without crack and that the stress criterion is strictly fulfilled for each possible crack. Figure 5 shows an example of the variation of the equivalent stress along the fiber-matrix interface at the reference fiber top half circle for a given sample geometry and varying strength ratios. Two configurations can be encountered depending on the tensile and shear strength magnitudes (that is a variation in the $\mu$ ratio):

* Tension-induced crack initiation leading to one equivalent stress maximum per half circle;

* Shear-induced crack initiation leading to two equivalent stress maxima per half circle.

More precisely, for $\mu$ larger than 1.5, the equivalent stress decreases monotonically from the maximum peak, inducing a single initiation site. On the contrary, for $\mu$ smaller than 1.5, the equivalent stress exhibits two local yet uneven maxima thus corresponding to two possible (yet non equiprobable) initiation sites. It is worth mentioning that experimental observations often show a single initiation site at one or both half circles of the fiber, which could indicate tension induced debonding initiation for the studied fiber-matrix interface. Nevertheless, we remind that without a high speed camera that would enable capturing either tension- or shear-induced crack initiation, it is not possible to draw a clear conclusion on this statement. Additionally, the many stochastic effects resulting from imperfections or heterogeneity could also influence debonding location. As a consequence, configurations with smaller $\mu$, and thus shear-induced crack initiation, are also studied for comparison purposes.

### 3.3.1. Tension-induced crack initiation

In case of tension-induced crack initiation ($\mu > 1.5$), the equivalent stress reaches a single maximum, which leads to a straightforward definition of the possible debonding configurations based on equivalent stress isocontours. Figure 6a shows equivalent stress to normal strength over the debonding angle corresponding to $\theta_i + \Delta\theta^+$ and $\theta_i - \Delta\theta^-$ according to the stress condition introduced by Leguillon [6]. This curve is obtained as the minimum stress attained over the
<figure><img src="image_7.png" /><figcaption>Figure 5: Equivalent stress to tensile strength ratio variation as a function of the angular position along the reference fiber top interface for various $\mu$ values.</figcaption></figure>

whole crack path prior to debonding. The equivalent stress to tensile strength ratio decreases with increasing both angle increments from $\theta_i$. The possible initiation debonding configurations based on the stress isocontours are depicted by the red line in Figure 6. The admissible initiation debonding angles correspond to configurations where $\sigma_{eq}$ is larger than $\sigma_c$. Figure 6b shows the debonding shape based on the stress criterion isocontours for several values of $\mu$ larger than or equal to 1.5. $\Delta\theta_{TC}^+$ indicates the transition angle between the areas undergoing tension or compression. Note that the range of $\mu$ presented in Figure 6b results in $\theta_i$ differences smaller than 3 degrees. The debonding shape is thus slightly influenced by the interface shear to tensile strength ratio. Ultimately, the increase of $\mu$ results in a debonding shape close to that obtained with a normal stress criterion (corresponding to $\mu = \infty$).

The average difference between the shapes for a normal stress criterion and the equivalent stress criterion for various $\mu$ values is evaluated. By comparing for fixed $\Delta\theta^+$, the difference in $\Delta\theta^-$ between debonding shapes obtained using the equivalent stress ($\Delta\theta_{\sigma_{eq}}^-$) or the normal stress criterion ($\Delta\theta_{\sigma_{nn}}^-$) is evaluated using Equation 2.

$$ \overline{\text{diff}} = \left( \frac{1}{\Delta\theta_{\text{TC}}^{+}} \int_{0}^{\Delta\theta_{\text{TC}}^{+}} \frac{|\Delta\theta_{\sigma_{\text{eq}}}^{-} - \Delta\theta_{\sigma_{\text{nn}}}^{-}|}{\Delta\theta_{\sigma_{\text{nn}}}^{-}} \right) \times 100 \quad (2) $$

For all the configurations covered, the average difference turns out to be no larger than 6%. Only debonding angles located in the tensile zone are considered because a crack is considered not to initiate under compressive stress. All in all, whatever the shear to tensile strength ratio, the normal stress-based debonding configuration provides a good estimate of the stress-based debonding shape obtained based on the equivalent stress criterion. Importantly, the advantage of using the debonding shape based on the normal stress rather than on the equivalent stress is that it allows defining possible initiation configurations that do not depend on fracture properties. As a consequence, it is numerically more efficient for inverse identification approach for instance. The same study is performed for several angles and inter-center distances to ensure that this observation is valid for any fiber inter-center distances and angles and the maximum differences never exceeded
<figure><img src="image_8.png" /><figcaption>Figure 6: (a) Equivalent stress to tensile strength variation as a function of angle increments (b) Stress isocontours for several values of $\mu$ where tension zone ($\sigma_{nn} > 0$) and compression zone ($\sigma_{nn} < 0$) are highlighted. The dashed lightgray line corresponds to symmetric configuration.</figcaption></figure>

10%. It is worth mentioning that the difference vanishes for a symmetric configuration, i.e. inter-center angle equal to 90 degrees.

### 3.3.2. Shear-induced crack initiation

For $\mu$ smaller than 1.5, the equivalent stress shows two local maxima (Figure 5). This configuration could then lead to various initiation scenarios:

* Single debonding initiation at one stress maximum;

* Two debonding initiations at both maxima.

García *et al.* [10] showed that two simultaneous debonding initiations at both half circles of a single fiber were not favorable compared to a single initiation site. Similar approach is employed here to compare the IERR of a single or double initiation located on the same half circle of the fiber. The same isovalues are determined based on each stress maximum and the nodes are released on the corresponding angles. Table 3 summarizes the IERR obtained for a randomly selected isovalue for $\mu = 0.5$ where $\alpha = 45$ deg. and $d = 2.5r$.

<table><caption>Table 3: IERR for a given initiation configuration considering single initiation site separately or both sites.</caption><thead><tr><th></th><th>Both maxima</th><th>$\theta_f \simeq 40$ deg. maximum</th><th>$\theta_f \simeq 132$ deg. maximum</th></tr></thead><tbody><tr><th>$G_{\text{inc}}$ [N/mm]</th><td>$3.97 \cdot 10^{-3}$</td><td>$4.13 \cdot 10^{-3}$</td><td>$3.7 \cdot 10^{-3}$</td></tr></tbody></table>

The debonding configuration maximizing the IERR can be considered as the favorable initiation site. Thus, a single initiation site is energetically more favorable. Interestingly, it corresponds to the location where the equivalent stress is maximum (see Figure 5). Finally, in case of shear-induced crack initiation, possible debonding configurations can be determined following stress
isocontours by considering only one debonding location that corresponds to the global equivalent stress maximum.

### 3.4. Debonding configuration based on energy criterion

The second approach consists in determining the debonding shape based on the configurations maximizing the IERR. To this end, nodes at the fiber-matrix interface are successively released on both sides of $\theta_i$ in increments of $\Delta\theta$ to compute the IERR for a given debonding angle ($\theta_d$) using Equation 3, where W denotes the elastic strain energy of the model and r the fiber radius:

$$G_{\text{inc}}(\theta_d, \sigma^\infty) = \frac{W(0, \sigma^\infty) - W(\Delta\theta^+ + \Delta\theta^-, \sigma^\infty)}{r(\Delta\theta^- + \Delta\theta^+)} \quad (3)$$

The local critical ERR is evaluated at both crack tips by means of the Hutchinson and Suo relation [9]. This relation links the interface opening ($G_{\text{IC}}$) and shear ($G_{\text{IIC}}$) critical energy release rates to the local critical ERR ($G_c$) using the mode mixity defined as $\psi(\theta) = \arctan(|\tau_{nt}|/\sigma_{nn})$ (Equation 4). The mode mixity is assessed at a distance of 0.015 mm from the two crack tips for each configuration according to [8] and the maximum of the two $G_c$ is retained.

$$G_c(\psi(\theta_d)) = G_{\text{IC}} [1 + \tan^2 [(1 - \lambda) \psi(\theta_d)]] \quad (4)$$

The parameter $\lambda$ is linked to the ratio between $G_{\text{IC}}$ and $G_{\text{IIC}}$. The average interface critical ERR, $\overline{G}_c$, corresponds to the average of the local interface critical ERR over the complete debonding shape, corresponding to $\theta_d$ (Equation 5).

$$\overline{G}_c = \frac{1}{\theta_d} \int_{\theta_i - \Delta\theta^-}^{\theta_i + \Delta\theta^+} G_c(\psi(\theta_d)) d\theta_d \quad (5)$$

Figure 7a shows the IERR to critical ERR ratio as a function of debonding increments on each side of the debonding initiation location. For a small debonding angle, the energy released is too small to satisfy the criterion. By increasing the debonding angle, the criterion is fulfilled and the surface exhibits a maximum, corresponding to an optimal angle maximizing the criterion. Therefore, this maximum corresponds to the favorable condition for the debonding to initiate, meaning that the loading to apply in order to fulfill the energy criterion is minimal. The debonding shape can be evaluated by maximizing the surface gradient until the maximum is reached, denoted by the blue marker in Figure 7.

The $G_{\text{IIC}}$ to $G_{\text{IC}}$ ratio, related to the parameter $\lambda$, influences the shape of the debonding as well as the configuration maximizing $G_{\text{inc}}/\overline{G}_c$. Figure 7b shows the shapes obtained from three different values of $\lambda$ for the studied configuration. The debonding initiation angle remains the same for a fixed value of $\lambda$. Setting $\lambda = 1$ results in $G_{\text{IC}} = G_{\text{IIC}}$ and a smaller value of $\lambda$ corresponds to $G_{\text{IC}} < G_{\text{IIC}}$. Decreasing $\lambda$ leads to an increase in the critical ERR and $G_{\text{IIC}}$ for a fixed $G_{\text{IC}}$ and mode mixity. Therefore, no shape change is observed for a fixed value of $\lambda$ as the ratio between $G_{\text{IC}}$ and $G_{\text{IIC}}$ remains the same.

### 3.5. Coupled stress and energy criterion

Coupling both stress and energy conditions allows the determination of the minimum loading to apply in order to initiate a debonding as well as the corresponding debonding angle. When both conditions are involved, the initiation debonding is sought as the angle for which both criteria are
<figure><img src="image_9.png" /><figcaption>Figure 7: (a) IERR to critical ERR variation. (b) Debonding shape for several values of λ and the dashed lightgray line corresponds to symmetric opening around θᵢ.</figcaption></figure>

fulfilled for a minimum applied loading. As a consequence, it is likely that the initiation debonding configuration lies between the debonding shapes determined based on either the stress isocontours or maximizing the IERR. Figure 8 shows the applied remote stress required to fulfill both conditions ($\sigma_{req}^{\infty}$) for a range of angle increments. The loading required to fulfill the equivalent stress condition increases monotonically with increasing debonding angle whereas the one required to fulfill the energy condition reaches a minimum. Combining both criteria thus results in determining the debonding configurations and corresponding applied stress for which the CC is fulfilled. Among these debonding configurations, the initiation debonding configuration is determined as the one minimizing the applied stress.

The debonding configuration for which both criteria are fulfilled for a minimum loading (i.e. the initiation debonding configuration) is highlighted on Figure 8a and 8b with a green marker. The same approach is employed to determine the minimum required remote stress that follows either the stress-based or energy-based debonding configurations, also highlighted in Figure 8a and 8b. These minima differ slightly from the CC solution, both in terms of the required remote stress and debonding angle.

For given $\lambda$ and $\mu$, the debonding configuration minimizing the loading required to fulfill the CC actually depends on the interface brittleness number $\gamma$ which was introduced by Mantič [8] as follows (Equation 6):

$$ \gamma = \frac{1}{\sigma_c} \sqrt{\frac{G_{IC}E^*}{r}} \qquad (6) $$

Each marker in Figure 8b corresponds to the CC solution with a different $\gamma$, where $E^*$ denotes the harmonic mean of the effective elastic moduli of the two constituents. For small enough $\gamma$, the debonding configuration obtained using the CC is close to the debonding shape predicted based on stress isocontours, which is also close to the one determined by maximizing $G_{inc}/\overline{G_c}$. How-
<figure><img src="image_10.png" /></figure>

$\alpha = 45 \text{ deg.}, d = 2.5r$

<figure><img src="image_11.png" /></figure>

(a)

<figure><img src="image_12.png" /></figure>

$\gamma = 0.66$

<figure><img src="image_13.png" /></figure>

$\alpha = 45 \text{ deg. } d = 2.5r$

<figure><img src="image_14.png" /></figure>

(b)

<figure><img src="image_15.png" /></figure>

$\gamma = 1.16$

<figure><img src="image_16.png" /></figure>

(c)

<figure><img src="image_17.png" /></figure>

(e)

<figure><img src="image_18.png" /></figure>

(d)

Figure 8: (a) Required loading to satisfy stress or energy conditions where solid lines indicate the configurations corresponding to either stress-based or energy-based possible debonding shapes and circles indicates the corresponding initiation configurations minimizing the imposed loading. (b) Isocontours of the loading required to fulfill the CC emphasizing the minimum indicated in Fig. 8a and minima obtained for various other properties range, the smaller the gray intensity, the higher $\gamma$ is. Visualization of $\theta_i \simeq 100$ deg. and associated debonding angle at initiation, before any potential unstable growth, taken from Fig. 8b for several values of $\gamma$, (c) $\gamma = 0.66$ with $\theta_d = 33$ deg. (yellow circle), (d) $\gamma = 1.16$ with $\theta_d = 85$ deg. (green marker) and (e) $\gamma = 1.45$ with $\theta_d = 102$ deg. (orange circle).

ever, for large enough $\gamma$ the debonding configuration retrieves the debonding shape determined by maximizing $G_{\text{inc}}/\bar{G}_c$. Whatever the interface brittleness number, the debonding initiation configuration actually lies between the debonding shapes determined based on either stress isocontours or maximizing $G_{\text{inc}}/\bar{G}_c$.

Figure 9 shows the stress-based, energy-based and CC-based debonding configurations obtained for inter-center angles of 90, 60 and 0 degrees with $d = 2.5r$. For the 90 degrees inter-center angle, a symmetric configuration is observed with respect to the loading direction. The stress-based, energy-based and CC-based debonding shapes are similar since there is no asymmetric perturbation of the fields. Similarly to the previous example, for inter-center fiber angles different from 90 degrees
(Fig. 9b and Fig. 9c), the debonding configuration predicted using the CC always lies between the stress-based and energy-based debonding configurations. Moreover, it is also consistent with the stress-based debonding shape for small $\gamma$ and with the energy-based debonding shape for large $\gamma$. For intermediate $\gamma$, there is a transition zone for which the debonding configuration does not correspond to neither shapes. Figure 9b shows the predicted debonding configurations for a 60 deg. inter-center angle for which the transition from stress-based to energy-based occurs abruptly, i.e. all debonding configurations predicted by the CC are close either to the stress-based or the energy-based debonding configuration whatever $\gamma$. On the contrary, Figure 9c shows a smoother transition, where the debonding configuration predicted by the CC does not correspond to stress-based or energy-based debonding configuration for a significant range of $\gamma$ (from 0.96 to 1.36).

<figure><img src="image_19.png" /><figcaption>Figure 9: Comparison between debonding shapes for configuration where $d = 2.5r$ and (a) $\alpha = 90$ deg. (b) $\alpha = 60$ deg. and (c) $\alpha = 0$ degrees. The markers denotes the CC solutions for various $\gamma$ values and dashed lightgray lines correspond to symmetric openings on both sides from $\theta_i$ whose values are given in Table 2.</figcaption></figure>

It is noteworthy that even in the case of a strongly influenced initiation location (when $\theta_i$ is much larger or smaller than 90 deg.), the debonding initiation angle tends to recover a symmetric configuration with respect to the loading direction. Figure 9b highlights this last observation with a $\alpha = 60$ deg. then the increment of angle $\Delta\theta^-$ is small compared to the increment of angle $\Delta\theta^+$ tending towards an initiation angle centered at about 90 degrees. The same observation was obtained by Mantič et al. who numerically [33] found that, wherever debonding initiation might occur, propagation will tend towards the loading direction as the amount of available elastic strain energy is larger.

Stress isocontours allow efficient numerical prediction of the debonding shape because they can be derived from a purely linear elastic calculation without node release. The difference induced by the stress-based debonding shape from the favorable configuration can be evaluated. The stress-based solution is determined by locating the minimum on surface $\sigma_{\text{req}}^{\infty}/\sigma^{\infty}$ (Fig. 8a) that follows the stress-based debonding shape and directly comparing it to surface global minimum, corresponding to the favorable CC configuration. Figure 10a shows the variation of the relative difference in terms of remote loading between the two solutions as a function of $\gamma$ for inter-center angles of 60 and 0 degrees. The relative difference shows increasing trends with $\gamma$. The maximum relative
difference remains smaller than 5.1 % and the maximum is reached for the largest values of $\gamma$ ($\gamma = 1.54$) when the configuration tends towards energy-based debonding shape. It is noticeable that the difference induced by the stress remains small and arguably not worth the computation time cost it requires. Therefore, basing the crack shape on the stress isocontours could be a good compromise to determine the debonding initiation remote loading. Figure 10b shows the location of the debonding based on the CC and the stress at initiation for $\alpha = 60$ deg. and $\gamma = 1.54$. The debonding angles at initiation predicted by the two shapes are identical, but their locations differ slightly with respect to $\theta_i$. The location of the debonding based on the CC solution is equivalent to the energy solution in this case, tending to maximize the energy released by being symmetric with respect to the loading direction.

<figure><img src="image_20.png" /><figcaption>Figure 10: (a) Relative difference between CC and stress-based debonding shape solution as a function of $\gamma$. (b) Schematic comparison of the debonding location and angle at initiation predicted based on both stress ($\theta_d \simeq 100$ deg.) and CC ($\theta_d \simeq 97$ deg.) shapes with $\theta_i \simeq 60$ degrees.</figcaption></figure>

### 3.6. Unstable crack propagation

After debonding initiation, it is possible to assess its propagation based on Griffith's criterion $G(\theta_d) = -dW/d\theta_d \ge G_c(\psi(\theta_d))$. The ERR $G$ can actually be obtained from $G_{inc}$ based on the following relation (Equation 7), where $S$ denotes the debonding surface i.e. $r\theta_d$ or $r(\Delta\theta^+ + \Delta\theta^-)$:

$$G(\theta_d) = G_{inc}(\theta_d) + S \times \frac{dG_{inc}}{dS} \qquad (7)$$

When $G_{inc}/\overline{G_c}$ maximum is reached, $d(G_{inc}/\overline{G_c})/dS = 0$ so that $G_{inc}/\overline{G_c} = G/G_c$ (see Appendix A). Therefore, two configurations can be encountered after the initiation of debonding. On the one hand, if $d(G_{inc}/\overline{G_c})/dS \ge 0$, then $G/G_c$ is larger than 1 which implies further unstable debonding phase without increasing the applied load until $G/G_c = 1$. On the other hand, $d(G_{inc}/\overline{G_c})/dS = 0$ for the initiation debonding configuration so that $G_{inc}/\overline{G_c} = G/G_c = 1$, resulting in no further propagation of debonding after initiation, as $d/dS(G/G_c) < 0$ (see Appendix B).

The debonding has two different crack ends, which may separately or simultaneously propagate. One approach to identify the unstable debonding shape is to consider the direction that maximizes
the ERR to critical ERR ratio. This enables the determination of the arrest angle corresponding to an angle for which $G/G_c$ becomes smaller than 1.

Figure 11a shows the unstable propagation phase of the crack in the configuration driven by both stress and energy for $\gamma = 1.1$. The debonding propagates in an unstable manner from the debonding angle solution, at initiation, obtained using the CC to the angle where $G/G_c$ becomes smaller than 1, after unstable propagation phase. Figure 11b shows the unstable propagation shape of the debonding for several values of $\gamma$. Higher values of $\gamma$ result in an energy-driven configuration where no further propagation of the debonding is possible after initiation ($\gamma = 1.5$). On the opposite, a large unstable debonding phase is encountered when $\gamma$ decreases, leading to a larger arrest debonding angle ($\gamma = 1$). Although the initiation debonding angle differs from the location maximizing $G/G_c$, the unstable debonding retrieves an identical shape regardless of the interface properties.

<figure><img src="image_21.png" /><figcaption>Figure 11: (a) ERR to critical ERR variation. (b) Isocontours of the critical ERR ($G_c$) and unstable debonding shapes obtained for several values of $\gamma$.</figcaption></figure>

Configuration with 0 deg. inter-center angle highlights a transition zone where neither stress-based nor energy-based debonding shapes describe the CC solution (see Figure 9c). By computing the favorable solution that minimizes the loading necessary for initiation based on either stress, energy or CC shape, several starting angles can be used to assess the unstable debonding phase. Figure 12 shows the three corresponding debonding shapes. As the required remote loadings for the three solutions have relative differences smaller than 1%, the ERR to critical ERR ratio surface presented is based on the CC solution for illustration purposes. However, the ERR calculation differs because the three solutions have different initiation remote loadings, so that the shape could be slightly different. All cases provide a similar debonding arrest angle with the same location regarding $\theta_i$ even with a slightly different initiation remote loading. Moreover, all shapes retrieve the same trajectory that maximizes the ERR to critical ERR ratio. Regardless of the type of shape used (stress, energy, CC), a configuration driven by both stress and energy will propagate to
approximately the same debonding arrest angle, as the remote loading required for both solutions are close. Some limitations exist for energy-driven configuration, where a stress-based debonding shape solution could lead to an unstable propagation phase while the energy-based solution leads to a stable propagation. However, the final solution can be obtained using the local ERR released for each configuration, that is computationally expensive to determine.

$$ \alpha = 0 \text{ deg. } d = 2.5r \gamma = 1.16 $$

<figure><img src="image_22.png" /><figcaption>Figure 12: Unstable debonding shape obtained from the stress-based, energy-based and CC-based initiation debonding angle. Table shows the associated debonding initiation angle and arrest angle obtained after unstable propagation phase.</figcaption></figure>

## 4. Nearby fiber influence on debonding initiation

Perturbation induced by the neighboring fiber influences the debonding location along the matrix-fiber interface. However, fields intensity are also influenced by the neighboring fiber, both in terms of stress and energy, depending on the geometrical parameters *d* and *α*. Stress fields at the fiber-matrix interface are therefore either intensified or damped and the same observation can be made concerning the elastic strain energy released by the debonding initiation. For similar interface properties, debonding initiation location (top or bottom) and debonding angle range vary as a function of both inter-center angle and distance.

Figure 13a shows the variation of the maximum normal to far field stress ratio evaluated at both top and bottom half circles of the reference fiber as a function of the neighboring fiber position (varying *d* and *α*). The bottom half circle is less affected than the top half circle. The latter is more affected by the neighboring fiber, highlighting a larger variation in the stress field. Also, with the distance between the top half circle and the neighboring fiber reduces, larger stress perturbations occur. The stress intensification reaches a maximum for 90 deg. configuration. For a configuration close to 45 deg., the top half circle stress decreases, resulting in a higher stress at the bottom half circle. Figure 13b shows the variation of the maximum IERR to critical ERR ratio which is influenced similarly to the stress fields. The latter is obtained by releasing debonding based on
the tensile stress isocontours. Higher IERR to critical ERR ratio is observed for angles close to 90 degrees, where the stress intensity is maximum. However, for inter-center angle range from 20 to 60 degrees, IERR to critical ERR ratio appears to be higher at the bottom half circle. Based on the latest observations, it is more likely that debonding would occur at the top half circle for inter-center distances close to 90 degrees at a lower imposed loading. The smaller the distance between the centers, the smaller the load required for debonding to initiate. Conversely, for an inter-center angle close to 45 degrees, debonding could most likely occur at the bottom half circle.

<figure><img src="image_23.png" /><figcaption>Figure 13: (a) Maximum normal to far-field stress ratio variation and (b) maximum IERR to average critical ERR ratio variation as a function of the inter-center angle and distance, for an arbitrary remote loading.</figcaption></figure>

The perturbation of the stress and energy fields at the fiber-matrix interface therefore has a significant influence on:

* The debonding initiation location;

* The required loading for debonding to initiate.

In the following, the CC solution is compared to experimental observations for different geometrical configurations.

### 4.1. *Comparison with experimental observations*

The interface properties of a single fiber sample made of the same constituents were previously determined by inverse identification [17]. The normal and shear strengths were selected within a range of acceptable strengths. Table 4 summarizes the identified interface fracture properties. An application of the CC for all configurations is carried out in the following. The debonding initiation stresses at the top and bottom fiber half circles obtained using the CC and debonding shapes based on stress isocontours are compared with experimental measurements. Furthermore, considering the relevance of the CC, the comparisons will allow to challenge or validate the accuracy of the interface properties determined previously [17].

The debonding initiation sites (reference fiber top or bottom half circle) are compared with the experimental observations in Table 5. For various geometrical configurations, debondings that
<table><caption>Table 4: Fiber-matrix interface properties from [17].</caption><thead><tr><th>G<sub>IC</sub> [N/mm]</th><th>G<sub>IIC</sub> [N/mm]</th><th>λ</th><th>σ<sub>c</sub> [MPa]</th><th>τ<sub>c</sub> [MPa]</th></tr></thead><tbody><tr><td>0.0037</td><td>0.09</td><td>0.13</td><td>3.0</td><td>5.0</td></tr></tbody></table>

<table><caption>Table 5: Debonding location (T for top and B for bottom) for several inter-center distances (columns) and angles (lines). Within a single cell, the upper right term corresponds to the experimental result and the lower left value to the computed CC one.</caption><thead><tr><th>d</th><th>α</th><th>0 deg.</th><th>30 deg.</th><th>45 deg.</th><th>60 deg.</th><th>90 deg.</th></tr></thead><tbody><tr><td rowspan="2">2.5r</td><td></td><td>T</td><td>T</td><td>T</td><td>T</td><td>B</td></tr><tr><td></td><td>T or B</td><td>B</td><td>B</td><td>B</td><td>T</td></tr><tr><td rowspan="2">4r</td><td></td><td>B</td><td>T</td><td>B</td><td>T</td><td>B</td></tr><tr><td></td><td>T or B</td><td>T</td><td>T or B</td><td>B</td><td>T or B</td></tr><tr><td rowspan="2">5r</td><td></td><td>T</td><td>B</td><td>T or B</td><td>B</td><td>B</td></tr><tr><td></td><td>T or B</td><td>T</td><td>T or B</td><td>B</td><td>T or B</td></tr></tbody></table>

are found to be favorable at the top half circle are denoted as “T” and at the bottom half circle as “B”. Concerning the simulations, when the difference between the top and bottom half circles is less than 1% relative error, the location of the debonding is designated as “T or B” since the position of the debonding may take place at either half circle. The debonding locations show no particular trend as a function of the angle and distance inter-center. The experimental results show poor agreement with the CC solution, especially for an inter-center angle of 90 degrees. The difference observed may be associated with the idealized loading employed numerically whereas the experimental loading can differ from a simplified far-fields loading. Additionally, 2D DIC allows one single free edge to be observed at a time, while the initiation of debonding may occur on the other side due to geometric imperfections for instance.

Figure 14 shows the load required to initiate debonding as a function of the inter-center angles. These values are evaluated experimentally and using the CC for three different inter-center distances. Uncertainties of the experimental results are obtained by tracking a sudden change of slope in the extension of the virtual gauges located at each fiber pole [31]. Two samples per configuration are studied, resulting in 8 different values. The results of each pole are averaged and the two extremes are then used for the estimation of the uncertainties.

For all angles and inter-center distances, the CC solutions are indicated by round and square markers for the top and bottom position, respectively. The CC solution that requires a minimum load is the favorable configuration. The latter depends on the inter-center angle. For angles close to 0 degrees, both the top or bottom locations are favorable (difference smaller than 1%). For angles between 30 and 60 degrees, the favorable initiation site is located at the bottom half circle. Finally, for angles close to 90 degrees, the top half circle seems to be the most favorable initiation site. The higher the inter-center distance, the smaller the difference between the top and bottom solutions, due to lower field perturbation, see Figures 13 and 14.

The overall comparison leads to a good agreement between the CC solution and the experimental observation. This observation attests to the accuracy of the interface properties determined previously. However, a larger difference is observed for an inter-center angle of 90 degrees for all three inter-center distances, possibly due to the aforementioned difficulties in the identification of
the debonding locations.

<figure><img src="image_24.png" /><figcaption>Figure 14: Loading required for debonding initiation as a function of the inter-center angle for (a) $d = 2.5r$ (b) $d = 4r$ (c) $d = 5r$. The scatter bars associated with the experimental data points represent variations between several independent tests.</figcaption></figure>

The experimental determination of the debonding angle at initiation after the potential unstable propagation is relatively difficult. The main difficulty lies in the identification of the locations of the two crack tips, possibly masked by the DIC speckles or considered as a process zone. Moreover, experimental capture of the unstable propagation phase remains impossible without the use of a rapid camera. In the sequel, experimental debonding angles are therefore compared to arrest angles obtained numerically. The CC solutions provide an almost identical debonding angle after initiation, close to 77 degrees, whatever the inter-center angle and distance. Figure 15a shows the variation of the debonding angles as a function of the inter-center angles. The experimental values correspond to the debonding angles averaged and their standard deviation observed at initiation for all inter-center distances with fixed inter-center angle. Only small variations are
seen and the experimental angles are close to 62 degrees. The latter is close to the averaged CC solution, with 24% relative difference, as the CC solutions slightly overestimate the experimental observations. Figure 15b shows the location and size of the debonding angles predicted by the CC using the stress isocontours and the corresponding experimental results, for the smallest distance ($d = 2.5r$), extracted from a half circle of one of the two samples. The CC solution accurately describes the location of the debonding by providing an upper limit to the experimental angle as observed on the overall trend (Figure 15a). The CC solution and the experimental angle locations and sizes are in good agreement, endorsing the stress isocontour as an effective approximation to describe the debonding shape. Moreover, even in the presence of the neighboring fiber, the experimental debonding is approximately centered at 90 degrees. The same result is obtained with the CC solution which validates the previous numerical observation (see Section 3.5) concerning the debonding initiation angle obtained with a strongly influenced $\theta_i$.

## 5. Conclusion

Debonding initiation and propagation were assessed using both the CC and Griffith criterion on a two-fiber sample subjected to remote uniaxial tensile loading. The debonding shape was predicted on the basis of stress- and energy-based isocontours compared to the CC solution. Finally, the results of the CC solution were compared with the corresponding experimental results.

Using normal stress-based or equivalent stress-based debonding configurations lead to an almost identical initiation configuration. Normal stress-based debonding configurations can be efficiently extracted without introducing any fracture properties and therefore appear as a relevant compromise to describe the stress-based debonding configurations.

The debonding angle size and location predicted by the CC lies between possible debonding configurations based on energy or stress isocontours. Depending on the interface properties, e.g. the interface brittleness number $\gamma$, three different configurations can be encountered:

* Small $\gamma$: the stress-based debonding prediction provides a reasonable approximation for the possible initiation debonding shapes;

* Intermediate $\gamma$: the possible initiation debonding shapes lie in-between the stress-based and energy-based debonding predictions;

* Large $\gamma$: the energy-based debonding prediction provides a good approximation, as the solution tends to be energy driven.

Despite the difference observed between the three possibilities, similar initiation remote loadings are obtained using either the stress-based, the energy-based or the CC-based debonding configurations. The stress-based debonding configurations therefore appear to be a good compromise as it is numerically much more efficient since it can be extracted based on a single calculation without crack. However, it is not possible to accurately describe the unstable debonding phase based on stress isocontours because the stress-based and energy-based debonding configurations may differ widely for some configurations. Therefore, the final location of the angle described by the two shapes could be slightly different. However, basing the debonding configuration on stress isocontours and assessing the unstable debonding phase on energy isocontours leads to the same debonding arrest angle, but is computationally more expensive.
<figure><img src="image_25.png" /><figcaption>Figure 15: (a) Experimental and numerical (CC) debonding arrest angles as a function of the inter-center angles. (b) Representation of the debonding angle and location predicted by the CC and obtained experimentally with $d = 2.5r$ from a half circle of one of the two samples. The corresponding DIC residual is shown for $\alpha = 60$ deg. to illustrate the experimental angle extraction process, more details can be found in [31].</figcaption></figure>

Blind application of the CC and comparison with experimental results showed good overall agreement both in terms of required remote load and debonding angle after propagation using previously identified interfacial shear and tensile strengths and critical opening ERR. The debonding angle after initiation is slightly overestimated by the CC, which could be explained by an un-
derestimation of the critical shear ERR that strongly influences the debonding arrest angle. The CC results obtained for a 90 degree inter-center configuration differ from the experimental results, which could be explained by the idealized loading applied numerically or by an experimental inter-center angle slightly different from 90 degrees. It is also important to mention that the properties used were determined using a 2D simulation without considering the singularity acting at the free edge. Therefore, the 3D simulation could allow a better modeling of the problem, especially of the 90 deg. geometrical configuration where the CC fails to accurately reproduce the experimental results.

## Appendix A. Proof of equality $G/G_c = G_{inc}/\overline{G_c}$ when $d/dS(G_{inc}/\overline{G_c}) = 0$

Recalling that $G_{inc}$, $\overline{G_c}$, $G$ and $G_c$ are depending on the debonding surface $S$. When, a debonding initiates for a debonding surface maximizing $G_{inc}/\overline{G_c}$:

$$ \frac{d}{dS} \frac{G_{inc}}{\overline{G_c}} = 0 \qquad (A.1) $$

Recalling A.1:

$$ \frac{d}{dS} \frac{G_{inc}}{\overline{G_c}} = \frac{\frac{dG_{inc}}{dS} \overline{G_c} - \frac{d\overline{G_c}}{dS} G_{inc}}{\overline{G_c}^2} \qquad (A.2) $$

It derives from the right part:

$$ \frac{dG_{inc}}{dS} \overline{G_c} = \frac{d\overline{G_c}}{dS} G_{inc} \qquad (A.3) $$

The IERR can be expressed as:

$$ G_{inc} = G - S \frac{dG_{inc}}{dS} \qquad (A.4) $$

Integrating A.4 into A.3:

$$ \frac{dG_{inc}}{dS} \overline{G_c} = \frac{d\overline{G_c}}{dS} \left( G - S \frac{dG_{inc}}{dS} \right) \qquad (A.5) $$

$$ \frac{dG_{inc}}{dS} \overline{G_c} = \frac{d\overline{G_c}}{dS} G - S \frac{d\overline{G_c}}{dS} \frac{dG_{inc}}{dS} \qquad (A.6) $$

$$ \frac{dG_{inc}}{dS} \overline{G_c} + S \frac{d\overline{G_c}}{dS} \frac{dG_{inc}}{dS} = \frac{d\overline{G_c}}{dS} G \qquad (A.7) $$

$$ \frac{dG_{inc}}{dS} \left( \overline{G_c} + S \frac{d\overline{G_c}}{dS} \right) = \frac{d\overline{G_c}}{dS} G \qquad (A.8) $$

The local critical ERR can be expressed as:

$$ G_c = \overline{G_c} + S \frac{d\overline{G_c}}{dS} \qquad (A.9) $$

Replacing into A.8:

$$ \frac{dG_{inc}}{dS} G_c = \frac{d\overline{G_c}}{dS} G \qquad (A.10) $$
$$ \frac{dG_{\text{inc}}}{dS} = \frac{d\overline{G_c} G}{dS \overline{G_c}} \qquad (A.11) $$

Remembering that:

$$ \frac{d G_{\text{inc}}}{d S \overline{G_c}} = 0 \qquad (A.12) $$

So it can be stated that:

$$ \frac{dG_{\text{inc}}}{dS} = \frac{d\overline{G_c} G_{\text{inc}}}{dS \overline{G_c}} \qquad (A.13) $$

Finally, the following equality is obtained by replacing A.13 into A.11:

$$ \overline{G_c} = \frac{G}{G_{\text{inc}}} \qquad (A.14) $$

## Appendix B. Proof of inequality $d/dS(G/G_c) < 0$ when $S = S_{\min}$

Recalling that $S_{\min}$ denotes the debonding surface for which $d/dS(G_{\text{inc}}/\overline{G_c}) = 0$.

$$ \frac{d G}{d S \overline{G_c}} = \frac{\frac{dG}{dS} G_c - \frac{dG_c}{dS} G}{G_c^2} \qquad (B.1) $$

As $G = G_c$ for $S = S_{\min}$:

$$ \frac{d G}{d S \overline{G_c}} = \frac{\frac{dG}{dS} G_c - \frac{dG_c}{dS} G}{G_c} \frac{1}{G_c} = \left( \frac{dG}{dS} - \frac{dG_c}{dS} \right) \frac{1}{G_c} \qquad (B.2) $$

As demonstrated by Mantič [8], $\frac{dG}{dS} < \frac{dG_c}{dS}$ for $S = S_{\min}$, so the following equality is obtained:

$$ \frac{d G}{d S \overline{G_c}} = \left( \frac{dG}{dS} - \frac{dG_c}{dS} \right) \frac{1}{G_c} < 0 \qquad (B.3) $$

## References

[1] E. Mäder, K. Mai, E. Pisanova, Interphase characterization in polymer composites - monitoring of interphasial behavior in dependence on the mode of loading, Composite Interfaces 7 (3) (2000) 133–147, publisher: Taylor & Francis. doi:10.1163/156855400300185270. URL https://doi.org/10.1163/156855400300185270

[2] J. Karger-Kocsis, H. Mahmood, A. Pegoretti, Recent advances in fiber/matrix interphase engineering for polymer composites, Progress in Materials Science 73 (2015) 1–43. doi:10.1016/j.pmatsci.2015.02.003. URL https://www.sciencedirect.com/science/article/pii/S0079642515000328

[3] Z. Fritah, C. Drouet, F. Salles, O. Marsan, M. Aufray, Influence of Water on an Epoxy/Amine-Metal Interphase: A Combined DFT and Mixing Calorimetry Approach, ACS Applied Materials & Interfaces 15 (8) (2023) 11342–11352, publisher: American Chemical Society. doi:10.1021/acsami.2c22431. URL https://doi.org/10.1021/acsami.2c22431
[4] G. Barenblatt, The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks, Journal of Applied Mathematics and Mechanics 23 (3) (1959) 622-636. doi:10.1016/0021-8928(59)90157-1.
URL https://www.sciencedirect.com/science/article/pii/0021892859901571

[5] D. Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids 8 (2) (1960) 100–104. doi:10.1016/0022-5096(60)90013-2.
URL https://www.sciencedirect.com/science/article/pii/0022509660900132

[6] D. Leguillon, Strength or toughness? A criterion for crack onset at a notch, European Journal of Mechanics - A/Solids 21 (1) (2002) 61–72. doi:https://doi.org/10.1016/S0997-7538(01)01184-6.
URL https://www.sciencedirect.com/science/article/pii/S0997753801011846

[7] A. A. Griffith, VI. The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 221 (582-593) (1921) 163–198, publisher: The Royal Society. doi:10.1098/rsta.1921.0006.
URL https://doi.org/10.1098/rsta.1921.0006

[8] V. Mantič, Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion, International Journal of Solids and Structures 46 (6) (2009) 1287–1304. doi:10.1016/j.ijsolstr.2008.10.036.
URL https://www.sciencedirect.com/science/article/pii/S0020768308004630

[9] J. Hutchinson, Z. Suo, Mixed Mode Cracking in Layered Materials, in: J. W. Hutchinson, T. Y. Wu (Eds.), Advances in Applied Mechanics, Vol. 29, Elsevier, 1991, pp. 63–191. doi:10.1016/S0065-2156(08)70164-9.
URL https://www.sciencedirect.com/science/article/pii/S0065215608701649

[10] I. García, V. Mantič, E. Graciani, Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?, European Journal of Mechanics - A/Solids 53 (2015) 75–88. doi:10.1016/j.euromechsol.2015.02.007.
URL https://www.sciencedirect.com/science/article/pii/S0997753815000200

[11] V. Mantič, I. García, Crack onset and growth at the fibre-matrix interface under a remote biaxial transverse load. Application of a coupled stress and energy criterion, International Journal of Solids and Structures 49 (17) (2012) 2273–2290. doi:10.1016/j.ijsolstr.2012.04.023.
URL https://www.sciencedirect.com/science/article/pii/S0020768312001709

[12] M. Muñoz-Reja, L. Távara, V. Mantič, P. Cornetti, Crack onset and propagation at fibre-matrix elastic interfaces under biaxial loading using finite fracture mechanics, Composites Part A: Applied Science and Manufacturing 82 (2016) 267–278. doi:10.1016/j.compositesa.2015.09.023.
URL https://www.sciencedirect.com/science/article/pii/S1359835X15003449

[13] T. Gentieu, J. Jumel, A. Catapano, J. Broughton, Size effect in particle debonding: Comparisons between finite fracture mechanics and cohesive zone model, Journal of Composite Materials 53 (14) (2019) 1941–1954, publisher: SAGE Publications Ltd STM.
doi:10.1177/0021998318816471.

URL https://doi.org/10.1177/0021998318816471

[14] I. García, M. Paggi, V. Mantič, Fiber-size effects on the onset of fiber-matrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models, Engineering Fracture Mechanics 115 (2014) 96–110. doi:10.1016/j.engfracmech.2013.10.014.
URL https://www.sciencedirect.com/science/article/pii/S0013794413003408

[15] L. Távara, I. G. García, R. Vodička, C. G. Panagiotopoulos, V. Mantič, Revisiting the Problem of Debond Initiation at Fibre-Matrix Interface under Transversal Biaxial Loads - A Comparison of Several Non-Classical Fracture Mechanics Approaches, Key Engineering Materials 713 (2016) 232–235, publisher: Trans Tech Publications, Ltd. doi:10.4028/www.scientific.net/kem.713.232.
URL https://doi.org/10.4028/www.scientific.net/kem.713.232

[16] A. Doitrand, R. Henry, H. Saad, S. Deville, S. Meille, Determination of interface fracture properties by micro- and macro-scale experiments in nacre-like alumina, Journal of the Mechanics and Physics of Solids 145 (2020) 104143. doi:10.1016/j.jmps.2020.104143.
URL https://www.sciencedirect.com/science/article/pii/S002250962030377X

[17] H. Girard, A. Doitrand, B. Koohbor, R. Rinaldi, N. Godin, D. Long, J. Bikard, L. Trouillet-Fonti, Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties, Engineering Fracture Mechanics (2023) 109254.
doi:10.1016/j.engfracmech.2023.109254.
URL https://www.sciencedirect.com/science/article/pii/S0013794423002126

[18] M. Muñoz-Reja, L. Távara, V. Mantič, P. Cornetti, Influence of a neighbour fibre on the onset and growth of a fibre-matrix debond under biaxial loading. A study by Finite Fracture Mechanics at linear elastic interfaces, 21st European Conference on Fracture, ECF21, 20-24 June 2016, Catania, Italy 2 (2016) 2022–2029. doi:10.1016/j.prostr.2016.06.254.
URL https://www.sciencedirect.com/science/article/pii/S2452321616302657

[19] 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, Modeling of fracture and damage in composite materials 168 (2016) 58–75. doi:10.1016/j.engfracmech.2016.01.022.
URL https://www.sciencedirect.com/science/article/pii/S0013794416000503

[20] M. Velasco, E. Correa, F. París, Interaction between fibres in the transverse damage in composites, Engineering Fracture Mechanics 239 (2020) 107273. doi:10.1016/j.engfracmech.2020.107273.
URL https://www.sciencedirect.com/science/article/pii/S0013794420308560

[21] L. Zhuang, R. Talreja, J. Varna, Transverse crack formation in unidirectional composites by linking of fibre/matrix debond cracks, Composites Part A: Applied Science and Manufacturing 107 (2018) 294–303. doi:10.1016/j.compositesa.2018.01.013.
URL https://www.sciencedirect.com/science/article/pii/S1359835X18300137
[22] V. Kushch, S. Shmegera, P. Brøndsted, L. Mishnaevsky, Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading, Recent Advances in Micromechanics of Materials 49 (1) (2011) 17–29. doi:10.1016/j.ijengsci.2010.06.020.
URL https://www.sciencedirect.com/science/article/pii/S0020722510001278

[23] F. París, E. Correa, V. Mantič, Kinking of Transversal Interface Cracks Between Fiber and Matrix, Journal of Applied Mechanics 74 (4) (2006) 703–716. doi:10.1115/1.2711220.
URL https://doi.org/10.1115/1.2711220

[24] A. Doitrand, D. Leguillon, 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending, International Journal of Solids and Structures 143 (2018) 175–182. doi:10.1016/j.ijsolstr.2018.03.005.
URL https://www.sciencedirect.com/science/article/pii/S0020768318301136

[25] A. Doitrand, D. Leguillon, Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints, International Journal of Adhesion and Adhesives 85 (2018) 69–76. doi:10.1016/j.ijadhadh.2018.05.022.
URL https://www.sciencedirect.com/science/article/pii/S0143749618301416

[26] N. Carrere, A. Doitrand, E. Martin, D. Leguillon, Theoretical study based on 2D assumptions of the influence of small pores on crack initiation in adhesively bonded joints, International Journal of Adhesion and Adhesives 111 (2021) 102979. doi:10.1016/j.ijadhadh.2021.102979.
URL https://www.sciencedirect.com/science/article/pii/S0143749621001779

[27] A. Doitrand, C. Fagiano, N. Carrère, V. Chiaruttini, M. Hirsekorn, Damage onset modeling in woven composites based on a coupled stress and energy criterion, Engineering Fracture Mechanics 169 (2017) 189–200. doi:10.1016/j.engfracmech.2016.11.021.
URL https://www.sciencedirect.com/science/article/pii/S0013794416306348

[28] E. Martin, B. Poitou, D. Leguillon, J. M. Gatt, Competition between deflection and penetration at an interface in the vicinity of a main crack, International Journal of Fracture 151 (2) (2008) 247–268. doi:10.1007/s10704-008-9228-0.
URL https://doi.org/10.1007/s10704-008-9228-0

[29] I. García, B. Carter, A. Ingraffea, V. Mantič, A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics, Composites Part B: Engineering 95 (2016) 475–487. doi:10.1016/j.compositesb.2016.03.023.
URL https://www.sciencedirect.com/science/article/pii/S1359836816300476

[30] R. Livingston, B. Koohbor, Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements, Composites Part C: Open Access 7 (2022) 100229. doi:10.1016/j.jcomc.2022.100229.
URL https://www.sciencedirect.com/science/article/pii/S2666682022000020

[31] H. Girard, B. Koohbor, A. Doitrand, R. Livingston, Experimental characterization of in-plane debonding at fiber-matrix interface using single glass macro fiber samples, Composites Part A: Applied Science and Manufacturing 171 (2023) 107573. doi:10.1016/j.compositesa.2023.107573.
URL https://www.sciencedirect.com/science/article/pii/S1359835X23001495
[32] Z. Tomičević, S. Roux, F. Hild, Evaluation of fatigue crack network growth in cast iron for different biaxial loading paths via full-field measurements, International Journal of Fatigue 92 (2016) 281–303. doi:10.1016/j.ijfatigue.2016.07.013.
URL https://www.sciencedirect.com/science/article/pii/S0142112316302171

[33] V. Mantič, L. Távara, A. Blázquez, E. Graciani, F. París, A linear elastic-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads, International Journal of Fracture 195 (1) (2015) 15–38. doi:10.1007/s10704-015-0043-0.
URL https://doi.org/10.1007/s10704-015-0043-0