# Comparison between 2D and 3D fiber-matrix debonding simulation for inverse identification of interface fracture properties

Hugo Girard, Aurélien Doitrand, Behrad Koohbor, Renaud G Rinaldi, Nathalie Godin, Jérôme Bikard

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

## ► To cite this version:

Hugo Girard, Aurélien Doitrand, Behrad Koohbor, Renaud G Rinaldi, Nathalie Godin, et al. Comparison between 2D and 3D fiber-matrix debonding simulation for inverse identification of interface fracture properties. Journal of Theoretical, Computational and Applied Mechanics, 2025, 10.46298/jt-cam.12997. hal-04431332v2

HAL Id: hal-04431332

https://hal.science/hal-04431332v2

Submitted on 3 Jul 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.12997

HAL hal-04431332v2

**History**

Received Feb 5, 2024

Accepted Jan 27, 2025

Published Jun 30, 2025

**Associate Editor**

Julien RÉTHORÉ

**Reviewers**

Anonymous

Anonymous

**Open Review**

HAL hal-04996881

**Licence**

CC BY 4.0

©The Authors

# Comparison between 2D and 3D fiber-matrix debonding simulation for inverse identification of interface fracture properties

Hugo GIRARD¹,², Aurélien DOITRAND¹, Behrad KOOHBOR³, Renaud RINALDI¹,

Nathalie GODIN¹, and Jérôme BIKARD²

¹ Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, Villeurbanne, France

² Axel'One (Solvay), Saint Fons, France

³ Department of Mechanical Engineering, Rowan University, United States of America

Fiber-matrix debonding in single transverse fiber specimen is studied experimentally and numerically based on the coupled criterion for which various 2D and one 3D configurations are used. Debonding initiation and propagation are mainly due to normal opening stresses in a 3D model whereas shear stresses play a minor role contrary to a 2D front model, i.e. in a plane normal to the fiber main axis. The 3D model enables describing the free surface singularity similarly to a 2D side model, i.e. along the fiber main axis. The latter cannot represent the debonding arrest and stable propagation after initiation. Overall, a 2D front model under plane strain assumption provides the best description of debonding initiation loading level compared to the 3D model, yet for a larger debonding opening. Experimental debonding openings are determined using Digital Image Correlation, providing the debonding initiation remote loading and corresponding opening. Tensile strengths and critical energy release rates respectively slightly higher and in the same order of magnitude are identified in 3D, based on the debonding opening, compared to a 2D front model.

**Keywords:** fiber-matrix debonding, finite fracture mechanics, inverse identification, 3D, coupled criterion

## 1 Introduction

Determining the properties of fiber-matrix interface is crucial for characterizing the overall damage behavior of composite materials as observed in (Dève and Schmauder 1992; Tripathi and Jones 1998; Park et al. 2006; Johnson et al. 2012) in the case of long fiber system. These interfaces experience significant loading and often undergo debonding due to their relatively lower fracture properties compared to the matrix and the fibers. Consequently, fiber-matrix debonding is frequently the initial form of damage that occurs. Once debonding initiates at the fiber-matrix interface, it propagates and eventually kinks into the matrix. This process leads to the coalescence of micro-cracks in the matrix between the fibers, creating a weak zone in the microstructure. As a result, this process tends to govern the overall failure of the composite part. Subsequent mechanisms can result in critical damage to the composite structure, such as delamination between plies, fiber breakage, or even complete failure of a composite lamina, see (Singh and Talreja 2012). In this regard, the mechanical characterization of the properties of the fiber-matrix interface is necessary to effectively control damage in composite structures.

Several experimental methods were developed to assess interface properties. These approaches are reviewed in (AhmadvashAghbash et al. 2023) and primarily focus on shear interface characterization. Among the large set of tests available, one can mention: pull-out (Zarges et al. 2018) and microbond tests (Lee et al. 2024), single fiber fragmentation test (Meer et al. 2019), and micro-indentation test as described in (Mandell et al. 1980). However, the aforementioned methods lead to a significant variability in measured properties, as illustrated in (Pitkethly et al.
1993). Furthermore, all the approaches focus on shear properties and employ analytical formulas that do not account for the real geometry and, consequently, the realistic stress fields developing in the examined samples.

Other characterization approaches have emerged to determine tensile fracture properties. Indeed, transverse loading of unidirectional plies implies a tensile loading mode rather than a shear loading mode. In this regard, it is crucial to determine both interface opening and shear fracture properties. Single embedded fiber is usually employed under transverse loading, as seen in experiments performed by Koyanagi et al. (2009), coupled with Cohesive Zone Model (CZM) to evaluate the interface fracture properties. However, such specimen geometries involve non-negligible free edge effects where debonding initiation is often observed, requiring 3D numerical modeling. Ogihara et al. 2009 tested either *straight* or *cruciform* sample geometries to investigate and confirm this phenomenon. They observed that the free edges induce smaller loading at initiation compared to the cruciform specimen. Debonding effectively initiates from the free edge and propagates towards the center of the specimen. Cruciform specimens have been employed by several authors to determine tensile interface properties (Gundel et al. 1995; Tandon 2000; Tandon et al. 2002). The Broutman test can also be employed to characterize tensile interface properties, thanks to a specific specimen geometry. Some applications are detailed in (Broutman 1966; Ageorges et al. 1999). These approaches require the use of inverse identification to determine the interface properties, as seen in the following papers (Meurs et al. 1998; Koyanagi et al. 2009; Totten et al. 2016).

Several approaches exist in order to model fiber-matrix debonding initiation and propagation. CZM are usually employed (Koyanagi et al. 2009; García et al. 2014; Gentieu et al. 2019). They are particularly advantageous as they predict the debonding process autonomously and allow for several simultaneous debonding events (Kushch et al. 2011). However, such advantages imply a larger computational cost. The Coupled Criterion (CC), introduced by Leguillon (2002), can also assess debonding initiation, propagation and predict the associated loading and size, see (Doitrand et al. 2024). This approach was first applied to fiber-matrix debonding by Mantič (2009), and the approach was then extended to various configurations (García et al. 2015; Muñoz-Reja et al. 2016). The CC approach is numerically efficient, making it suitable for inverse identification. The CC reverts to Linear Elastic Fracture Mechanics (LEFM) when employed to assess further debonding propagation, as employed in (Sandino et al. 2016; Velasco et al. 2020).

Fiber-matrix interface properties were recently identified by Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti (2023) using the CC and an inverse identification approach. Experimental results showed an abrupt debonding initiation at a finite angle, followed by stable debonding propagation. Similar debonding processes were predicted by the CC. However, the inverse identification was based on the debonding angle observed at the free edges whereas the employed 2D plane strain model provided somehow results representative of a debonding at the sample middle plane but lacked the ability to capture the stress singularity at the free edges. Such singularity can be represented through 3D modeling of fiber-matrix debonding initiation and propagation. The 3D application of the CC raises new challenges, such as crack shape determination. Leguillon 2014 attempted to determine the crack shapes based on stress isocontours and proposed a comparison between 2D and 3D applications of the CC. Properties derived from both models were in the same order of magnitude. Similarly, Doitrand and Leguillon 2018b proposed a CC application in both 2D and 3D for scarf joint failure prediction, where the 3D approach more accurately predicted the failure loading. The 3D applications of the CC were also extended to several configurations such as laminate failure in (García et al. 2016) and woven composite in (Doitrand et al. 2017) and, more recently, to adhesive pores in (Carrere et al. 2021) or anisotropic fracture in nacre-like alumina (Duminy et al. 2023). In this regard, a 3D model of a single fiber specimen is developed herein to perform the inverse identification of interface fracture properties and quantify the accuracy of 2D models. The presentation of the 3D model is provided in Section 2, along with the associated challenges related to the modeling. In Section 3, the stress and energy fields extracted from different 2D models are compared to the results from the 3D model to assess their validity. The debonding openings are determined experimentally in Section 4 using DIC and provide the initiation remote loading and corresponding opening, facilitating inverse identification of the interface fracture
properties. Finally, the inverse identification is performed in Section 5, based on the debonding opening, using the 3D model, and the identified properties are compared to the ones obtained in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023).

## 2 3D numerical model

### 2.1 Model presentation

Fiber-matrix debonding is studied by means of tensile tests on a single 2 mm diameter glass fiber embedded in an epoxy matrix sample, see Figure 1. The sample gauge section is 13 mm wide and 6 mm thick. A thorough description of the sample preparation can be found in (Livingston and Koohbor 2022) and further details are also given in Section 5. The specimen is subjected to tensile displacement (U) along the (Oy) direction resulting in a far-field stress $\sigma^\infty$ in the gauge section of the specimen. The elastic properties of the constitutive materials used are summarized in Table 1.

<table><caption>**Table 1** Elastic properties of the matrix and fiber. The fiber properties are provided by the manufacturer and the matrix properties are averaged from 3 in-house measurements at room temperature from Livingston and Koohbor 2022.</caption><thead><tr><td>Property</td><td>Epoxy</td><td>Glass Fiber</td></tr></thead><tbody><tr><td>Young's 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.3</td></tr></tbody></table>

Along with the geometry of the sample (left), Figure 1 depicts the dimensions of the corresponding 3D numerical model (right), where the proximity of the free edge to the fiber-matrix interface induces a 2% difference in the stress fields compared to a fiber embedded in infinite medium.

<figure><img src="image_3.png" /><figcaption>**Figure 1** Tensile specimen dimensions and associated numerical models where *r* denotes the fiber radius, $\sigma^\infty$ the remote stress and $\theta_f$ the angular position from the top fiber pole along the interface.</figcaption></figure>

Both 2D front and side simplifications of the problem are also displayed. To reduce computational costs, the actual sample geometry can be indeed simplified. Symmetry boundary
conditions are applied to the (x = 0), (y = 0) and (z = 0) faces, with the origin (0, 0, 0) being located at the fiber center point. Displacements are applied to the top face to replicate the remote loading. The calculations are performed assuming linear elasticity and small deformations using linear hexahedral elements. In the sequel, debonding initiation is assessed based on both stress and energy criteria. Therefore, two mesh convergence campaigns are performed to account for both stress and energy aspects. Due to the elastic properties mismatch between the fiber and the matrix, there is a singularity at the free edge interface, so that the stress locally tends to infinity. To ensure that the mesh size is small enough to accurately describe the stress gradient, the interface stress isocontours location are compared for each mesh size. Convergence is achieved when the element size does not affect the stress isocontours thereby ensuring that the smallest debonding surface (see Section 2.2) is not influenced. A subsequent campaign is conducted to assess the mesh size used to compute the variation in the elastic strain energy. The size chosen ensures a difference in released elastic strain energy for fixed areas is smaller than 1.5% compared with a converged solution. It results in ≈ 15 µm mesh size along the fiber matrix interface.

## 2.2 Debonding shape determination

The CC implementation requires the definition of the debonding shape along which the crack is likely to initiate. Considering a 2D simulation results in symmetric initiation from the fiber pole which can be thus described by a single parameter, i.e., the debonding angle, see (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023). However, a 3D simulation involves a third dimension to describe the crack shape, i.e., the crack front topology, which is initially unknown. Monitoring the debonding shape along the thickness is possible experimentally (see (Koyanagi et al. 2009), Martyniuk et al. 2013). It only yields observations of already propagated debondings after initiation. Approaches have emerged to overcome this challenge. Using stress isocontours to determine the crack shape yielded promising results, as shown for instance by Leguillon (2014); Doitrand and Leguillon (2018a); Doitrand and Leguillon (2018b); Carrere et al. (2021). Recently, Girard et al. (2024) showed that, in 2D, determining the possible initiation shapes based on the stress isocontours represents a relevant solution. Indeed, it yields the optimal debonding shape (i.e., the debonding shape minimizing the initiation loading) for small enough interface characteristic lengths. For larger characteristic lengths, it yields differences in the initiation loading no larger than 5% compared to the optimal initiation shape. Since it is much more efficient numerically than determining debonding shape based on the energy isocontours, stress isocontours based debonding shape determination is adopted in the sequel.

Figure 2(a) shows the variation in normalized normal stress as a function of angular position ($r\theta_f$) and through-thickness direction $z$. The normal stress component increases and tends towards infinity when approaching the free edge ($-z = 3$ mm).

Likewise, maximum values are attained at the fiber pole ($r\theta_f = 0$ mm.rad.). The normal stress component stabilizes at a constant value when moving away from the free edge for a given angular position. Figure 2(b) shows the normal stress isocontours and highlights the singularity in the vicinity of the free edge. It should be mentioned that whether remote stress or displacement is imposed as a boundary condition, it has no significant influence on the variation of interface stresses, with differences smaller than 3% on the equivalent stress $\sigma_{eq}$, see Equation (1). Thus, in the sequel, only displacement boundary conditions will be considered. Stress isocontours delineate different potential debonding shapes where the normal stress is strictly larger than a specified stress level. Figure 3(a) shows a 3D representation of a stress isocontour, for a normal stress level of $\sigma_{nn}/\sigma^\infty = 1.5$, on the fiber interface represented by the grey surface. An isocontour, e.g. the solid yellow line, encompasses a S surface, e.g. the shaded yellow area, which corresponds to the debonded zone. The lower the stress level associated with an isocontour, the larger the debonding surface. Figure 3(b) shows the isocontour for a normal stress level of $\sigma_{nn}/\sigma^\infty = 0.4$ and, as expected, a larger debonding surface is obtained. Since the normal stress monotonically decreases from the free edges, it is likely that debonding initiates at the free edge and then propagates through the thickness along the fiber pole. The debonding thus evolves from a localized debonding in the free edge vicinity to the creation of a tunnel reaching the other free
<figure><img src="image_4.png" /><figcaption>**Figure 2** Analysis of the fiber-matrix interface stress with (a) Variation of the normalized normal stress ($\sigma_{nn}$) along the interface emphasizing the stress singularity for $-z = 3$ mm and (b) corresponding isocontours.</figcaption></figure>

edge debonding. This debonding process closely resembles what was experimentally observed by Martyniuk et al. (2013), for a single fiber sample solicited under tension.

<figure><img src="image_5.png" /><figcaption>**Figure 3** Debonding surface delimitation for stress isocontours depicted in Figure 2(b) for (a) a stress level $\sigma_{nn}/\sigma^{\infty} = 1.5$ so that the debonding remains close to the free edge and (b) $\sigma_{nn}/\sigma^{\infty} = 0.4$ resulting in the two debondings merging by tunneling between each others.</figcaption></figure>

The benefit of basing debonding shapes on normal stress isocontours is that no interface properties are a priori needed. Consequently, the determination of isocontours becomes significantly more efficient, as the crack shape remains independent of fracture properties. However, the fiber interface may experience both opening and shear modes from the pole to the equator of the fiber. The equivalent stress

$$ \sigma_{eq} = \sqrt{\sigma_{nn}^2 + \frac{1}{\mu^2}(\tau_{nt}^2 + \tau_{nz}^2)}, \quad (1) $$

where $\mu = \tau_c/\sigma_c$ describes the ratio between the shear strength $\tau_c$ and the tensile strength $\sigma_c$, is employed to account for the complex stress field and both phenomena. The quantity $\sigma_{nn}$ represents the normal tensile stress at the fiber-matrix interface, while $\tau_{nt}$ and $\tau_{nz}$ correspond to the in-plane and out-of-plane shear stresses respectively.

Similarly to the representation shown in Figure 2(a) focusing on the normal stress only, the variation in equivalent stress along the interface can be computed. Since both normal and equivalent stresses decrease monotonically from the free edges, an efficient approach to consider the equivalent stress using normal stress-based isocontours is to provide the equivalent stress values along each isocontour. This assumption results in a non-constant equivalent stress, and the
minimum stress value is retained, accordingly with the stress condition proposed by Leguillon (2002), recalled in Section 2.3.

Figure 4 shows the normalized equivalent stress as a function of the debonding surface for different values of $\mu$. The larger the value of $\mu$, the closer the variation in equivalent stress is to the variation in normal stress. Consequently, when $\mu$ is large, the equivalent stress isocontours revert to the normal stress isocontours. Isocontours based solely on normal stress therefore induce a smaller error. Reducing $\mu$ to 1 (e.g., $\sigma_c = \tau_c$) results in a slight deviation of the equivalent stress from the trends observed with larger $\mu$ values (e.g., 3 and larger). Using the variation of equivalent stress based on normal stress isocontours seems to be a suitable approach for enhancing the efficiency of the numerical model. Moreover, the disparity between equivalent stress and normal stress remains minimal, highlighting the insignificance of shear in calculating the equivalent stress in comparison to the normal stress component. As a result, the use of the equivalent stress may become unnecessary, and the identification of shear strength based on equivalent stress could hold little significance. Nonetheless, in the subsequent sections, equivalent stress condition is used to maintain consistency with previously developed 2D models (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023), where shear has a more pronounced influence.

<figure><img src="image_6.png" /><figcaption>**Figure 4** Minimum value of $\sigma_{eq}$ calculated along each isocontours based on $\sigma_{nn}$ variation along the interface (equivalent to isocontours based on $\sigma_{eq}$ with $\mu = \infty$).</figcaption></figure>

## 2.3 CC implementation

Once the potential initiation debonding shapes are determined, the CC can be implemented. In order to provide the initiation loading and the debonding size, the CC combines the two conditions

$$ \left\{ \begin{array}{l} \sigma_{eq}(\vec{x}, \sigma^{\infty}) \ge \sigma_c, \quad \forall \vec{x} \in \Gamma, \\ G_{\text{inc}}(S, \sigma^{\infty}) \ge \bar{G}_c(S). \end{array} \right. \qquad (2) $$

The equivalent stress $\sigma_{eq}$ has to exceed the tensile strength $\sigma_c$ at any location on the debonding area $\Gamma$ before initiation. Likewise, the minimum $\sigma_{eq}$ has to exceed $\sigma_c$ over the isocontour as the stress decreases monotonically. The quantity

$$ G_{\text{inc}}(S, \sigma^{\infty}) = \frac{W(0, \sigma^{\infty}) - W(S, \sigma^{\infty})}{S} \qquad (3) $$

refers to the incremental energy release rate (IERR) which is the variation of the elastic strain energy $W$ per unit crack surface $S$. $G_{\text{inc}}$ must be larger than the global critical ERR $\bar{G}_c$, which is defined in Equation (4), to fulfill the energy criterion.

### 2.3.1 IERR calculation

From a numerical point of view, cracks can be defined either based on the initial interface mesh, or after remeshing so that the new mesh topology coincides with the stress isocontours trajectory. The first method is easier to implement because it does not require remeshing, but also yields an
approximate, non-smooth crack shape depending on the mesh size. In this context, the IERRs obtained through both methods have to be compared.

In Figure 5(a), a comparison of the IERR obtained using two different methods is presented. For the first method, denoted ①, a converged mesh size of 15 µm is adopted. With the first

<figure><img src="image_7.png" /><figcaption>**Figure 5** Influence of the methodology adopted to describe the debonding shapes on (a) (a) IERR variation as a function of the debonding surface with the method ① and ② where DOF translates the number of degrees of freedom of the model. (b) Example of debonding surface based on both method ① and ② with relatively large mesh size to highlight the potential difference observed between both methods.</figcaption></figure>

method based on the original interface mesh, a constant mesh size along the interface is used as well as near the free edges and then gradually increased once the tunnel has started (see Figure 2(b)). The second method, involving remeshing and identified as ②, also undergoes a mesh convergence study. This validation is essential to ensure the validity of the comparison, given that this method yields a more accurate representation of the isocontours topology. Achieving a converged solution requires a similar 15 µm mesh size over the isocontours, see Section 2.1, and larger sizes are adopted elsewhere to be consistent with method ①. Both methods yield overall comparable results in terms of IERR. For small debonding surfaces, method ① yields a slightly lower IERR compared to method ②. This difference may emanate from the presence of a relatively small debonding surface near the free edges, making the surface difficult to be captured using a constant mesh size that may be too large. The remeshing method thus has an advantage because the mesh size is varied accordingly to the surface geometry. Contours of the debonding surfaces obtained using the two methods are depicted in Figure 5(b). Using a large mesh size effectively highlights the contrast between the surface contours derived from remeshing and the shape defined by the original mesh. Hence, the area encompassed by method ① might be larger for the same isocontour level, potentially explaining the observed shift in the IERR. For larger debonding surfaces, method ① aligns with method ② IERRs. Furthermore, the overall gain of the method ② over the method ① in terms of IERR is negligible, about 3% difference for converged solutions. This highlights the higher computational cost of the method ① to be converged, twice the DOFs number in the method ②, as it requires a finer mesh that contributes to an expanded numerical model.

The first method is easier to implement but requires a finer mesh to accurately capture the isocontours shapes. While the second method is more challenging to model but may involve larger mesh size to precisely describe the isocontour. As a result, the first method may involve a longer calculation time, whereas the second approach might require larger implementation time but is more computationally efficient. All in all, both methods remain relevant due to their negligible disparity; however, the second method could allow small debonding surfaces to be captured more accurately and reduce the number of DOF.
### 2.3.2 Critical ERR calculation

The global critical ERR $\overline{G}_c$ (Equation (4)) is determined by calculating the average of the local critical ERR $G_c$ over the entire debonding front:

$$ \overline{G}_c = \frac{1}{S} \int_{0}^{S} G_c(\psi(s))ds. \qquad (4) $$

The local critical ERR is calculated according to the formulation introduced by Hutchinson and Suo (1991):

$$ G_c(\psi(S)) = G_{IC} [1 + \tan^2[(1 - \lambda)\psi(S)]]. \qquad (5) $$

The parameter $\lambda$ is a function of the $G_{IC}$ to $G_{IIC}$ ratio, and more information about its influence can be found in Section 5.2 and in (Doitrand et al. 2024). The critical ERR is denoted $G_{IC}$ in opening mode (mode I) and $G_{IIC}$ in shear mode (mode II). The global and local critical ERRs depend on the mode mixity $\psi$ that is the ratio between shear and tensile stresses acting at the debonding front. This quantity is averaged at the interface, at a distance of 7.5 µm, corresponding to half the element size, beyond the crack front according to the method used in 2D by Mantič (2009) and also employed in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023).

Besides, a condition may be applied to the calculation of the mode mixity $\psi$ when the interface experiences compression in order to only account for the shear-induced propagation at this location:

$$ \psi(S) = \begin{cases} \arctan\left(\frac{|\tau(S)|}{\sigma_{nn}(S)}\right) & \text{if } \sigma_{nn}(S) \ge 0, \\ \frac{\pi}{2} & \text{if } \sigma_{nn}(S) < 0. \end{cases} \qquad (6) $$

In the above condition, the shear component $\tau(S)$ can be affiliated to i) in-plane shear $\tau_{nt}$ with respect to the free edge, ii) out-of-plane shear $\tau_{nz}$, iii) the maximum of the two shear stresses or iv) an average shear stress. Figure 6 shows the variation in both mode mixity $\psi$ (Figure 6(a)) and critical ERR $\overline{G}_c$ (Figure 6(b)) for the above calculation possibilities as a function of debonding surface.

<figure><img src="image_8.png" /><figcaption>**Figure 6** Comparison of different calculation possibilities for determining (a) mode mixity and (b) global critical ERR as a function of the debonding surface for the different calculation possibilities with the interface properties identified previously in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023): $G_{IC} = 3.7 \times 10^{-3}$ N/mm, $\lambda = 0.13$.</figcaption></figure>

Taking $\tau(S) = \max(\tau_{nt}, \tau_{nz})$ will increase the influence of shear on the calculation of the critical IERR since larger $\psi$ are considered. This results in an increase in the critical ERR, see
Figure 6(b), as the mode II critical ERR is often larger than the mode I one. Taking the out-of-plane shear into account increases the mode mixity for small debonding surface because the debonding is close to the free edges. In this region, the elastic mismatch between the fiber and the matrix induces a large out-of-plane shear stress compared to the normal stress. In addition, the in-plane shear acts mainly on the equators of the fiber interface, at a specific distance from the free edge. Therefore, taking into account the in-plane shear stress in the mode mixity calculation increases the critical ERR for a larger debonding surface. Finally, an average of the shear stresses $\tau(S) = 0.5(\tau_{nt} + \tau_{nz})$ will reduce the influence of shear by considering an intermediate value. It should be mentioned that opting for the maximum of the two shear stresses results in a more important assessment of the mode mixity, involving larger critical ERR. This is due to the initial evaluation of the mode mixity at multiple positions along the debonding front for each surface. Subsequently, these values are averaged by alternating the maximum shear stresses along the debonding front whether it is $\tau_{nt}$ or $\tau_{nz}$. This approach yields a higher value compared to considering in-plane or out-of-plane shear stresses independently. Consequently, the initiation process on a surface dominated by shear will either be dampened or intensified, depending on the selected mode mixity calculation. Interestingly, the variation in $\bar{G}_c$ shown in Figure 6(b) remains small since the mode mixity is close to 0 meaning that the interface experiences mostly mode I. Consequently, the selected mode mixity calculation has a negligible influence on the critical ERR determination in the present example. Furthermore, it is noteworthy that while the consideration of one critical shear ERR is included, there are potentially two critical ERRs that can be accounted for ($G_{IIC}^{nt}$ and $G_{IIC}^{nz}$). However, for the sake of simplicity, the subsequent analysis will focus on a single critical shear ERR, employing the maximum shear stresses for the mode mixity calculation.

### 2.3.3 ERR calculation

The CC provides the debonding initiation loading and surface. Two configurations can be encountered at initiation depending on the IERR and the critical ERR. On the one hand, $d(G_{inc}/\bar{G}_c)/dS = 0$ leading to a stable debonding after initiation such that the initiation debonding surface corresponds to the arrest surface. On the other hand, $d(G_{inc}/\bar{G}_c)/dS > 0$ leading to $G/G_c > 1$ and further propagation of the debonding without any increase in the remote load until $G/G_c = 1$ occurs. Thus, in the second configuration, the ERR is evaluated using

$$G = -\frac{dW}{dS} \qquad (7)$$

as the opposite of the elastic strain energy $W$ derivative with respect to the debonding surface $S$.

Consequently, debonding shape associated with unstable debonding propagation must be determined. Exploring all the shape possibilities starting from the debonding initiation one would be computationally time-consuming. As a result, for the sake of efficiency, in particular to perform an inverse identification, the debonding shapes based on stress isocontours are used.

# 3 Comparison between 2D and 3D models

## 3.1 Models and debonding opening

The use of a 2D numerical model instead of a 3D model enables a significant reduction in computational costs. Considering the problem of a single fiber subjected to remote tension can be modeled in 2D using either plane strain (PE) or plane stress (PS) assumptions. Two distinct 2D models are suitable for assessing debonding initiation and propagation. Figure 7(a) shows the front model which corresponds to a stiff circular cylindrical inclusion within a softer squared medium. This model, commonly employed for investigating such problems (Mantič 2009; García et al. 2015; Gentieu et al. 2019; Martin et al. 2020; Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023), allows monitoring the debonding angle ($\theta_d$) relative to the imposed displacement ($U$) or corresponding remote stress ($\sigma^\infty$). While it accurately describes the fields around the specimen mid-plane, it does not account for the singularity at the free edges caused by elastic property mismatch. Figure 7(b) shows the current 3D modeling of the problem which considers simultaneously the debonding angle, the debonding length ($l_d$) and the
debonding surface (S) as a function of the remote stress or the imposed displacement. Finally, Figure 7(c) shows the side model which provides information on the debonding length and offers a more accurate description of the free edge singularity compared to the front model.

<figure><img src="image_9.png" /><figcaption>**Figure 7** Schematic of (a) the 2D front model, (b) the 3D model and (c) the 2D side model. The 2D representations are depicted in color on the 3D scheme. The quantity $\delta_{nn}$ denotes the debonding normal opening extracted at the free edge, shared by all the models.</figcaption></figure>

It is worth noting that the only debonding characteristic shared by the three models is the maximum normal opening at the free edge, $\delta_{nn}$, see Figure 7.

Figure 8(a) shows the variation in normal opening as a function of debonding angle for the PS and PE front models and the 3D model. The two front 2D models exhibit similar relationship between the two quantities, presenting relatively large debonding openings for any fixed angle compared to 3D modeling. In contrast, the 3D model provides a small opening even for a complete debonding. The 3D model therefore rapidly provides large debonding angle at the free edge for fixed debonding opening, whereas the 2D models provide a more progressive increase. The previous observations could be attributed to the debonding shape that exhibits a large debonding angle associated to small debonding length. A more concave than convex trend could result in a better correspondence between the 2D and 3D models.

Figure 8(b) shows the variation in opening as a function of debonding length for the side and 3D models. For a small debonding length, the side and 3D models show good correspondence. Beyond a certain point ($l_d/h \approx 0.05$), corresponding to the initiation of the tunnel effect, the 3D opening stops increasing, while the 2D front models continue to increase. This tunnel effect leads to complete debonding across the sample thickness, resulting in a constant debonding length increase for larger debonding openings. Once the tunneling starts, the debonding propagation influence on the opening becomes negligible. This effect is not considered by the side model, which predicts a progressive increase in opening as a function of debonding length.

<figure><img src="image_10.png" /><figcaption>**Figure 8** Comparison of the crack opening ($\delta_{nn}$) as a function of (a) the debonding angle ($\theta_d$) for the 3D and 2D front models and (b) the normalized debonding length ($l_d/h$) for the 3D and 2D side models.</figcaption></figure>

As mentioned already, the debonding opening is the only parameter shared by the three
models. It will thus be used to compare the stress and energy fields between the 2D front, 2D side and 3D models, as well as a full CC comparison, which are detailed hereafter.

## 3.2 Stress comparison

Figure 9 shows the variation of equivalent stress ($\sigma_{eq}$) as a function of debonding opening ($\delta_{nn}$) for the front, side and 3D models.

<figure><img src="image_11.png" /><figcaption>**Figure 9** Equivalent stress ($\sigma_{eq}$) as a function of the debonding opening ($\delta_{nn}$) for a similar (a) imposed displacement and (b) remote stress.</figcaption></figure>

The subsequent results are presented for homogeneous fixed material and interface properties, specifically either plane strain or plane stress is used. Consequently, the remote stress varies according to models for a similar imposed displacement, as indicated in Table 2. These values are calculated as the total reaction forces divided by the cross-section. In the following, all the calculations are performed with a prescribed displacement, which magnitude is adapted to obtain a given remote stress. Another possibility would be to impose a specific remote stress as a boundary condition (not considered here). For the sake of simplicity, only the results for the side PE model are provided since the trends are similar for the side PS model.

<table><caption>**Table 2** Remote stress extrapolated for each model under an imposed displacement of $U = 0.03L$ mm.</caption><thead><tr><th>Model</th><th>3D</th><th>2D Front PE</th><th>2D Front PS</th><th>2D Side PE</th></tr></thead><tbody><tr><td>Remote stress ($\sigma^{\infty}$) [MPa]</td><td>73.2</td><td>86.7</td><td>72.7</td><td>104.3</td></tr></tbody></table>

As a result, differences in both stress and energy fields are obtained, whether a configuration with given displacement magnitude (Figure 9(a)) or remote stress level (Figure 9(b)) is considered.

Contrary to the convex trends of the front models, the side model varies in a concave manner. The side model effectively captures the stress singularity originating from the free edges, exhibiting a pronounced downward trend. Furthermore, the side model concave trend is similar to the 3D model one. Both front models also show a decreasing convex trend as a function of the debonding opening, but converge to a finite stress value for vanishing debonding opening. Consequently, the front models are not preferable for accurately replicating the actual stress trends observed in 3D. Nonetheless, for a similar imposed displacement, the front PE model exhibits good correspondence in terms of stress levels with the 3D model for $\delta_{nn} > 20$ µm, with a relative difference smaller than 15%. For a similar remote stress, the side model better matches the 3D stress levels, exhibiting a similar downward trend.

## 3.3 Energy comparison

Figure 10 shows the variation of IERR as well as the IERR to critical ERR ratio as a function of debonding opening for the four models, all characterized by identical $G_{IC}$ and $\lambda$.
<figure><img src="image_12.png" /><figcaption>**Figure 10** IERR as a function of the debonding opening for a similar (a) imposed displacement and (b) remote stress. IERR to critical ERR ratio as a function of the debonding opening for a similar (c) imposed displacement and (d) remote stress.</figcaption></figure>

All of the above quantities are shown for either the same applied displacement or remote stress, left or right column, respectively. For a similar imposed displacement (Figure 10(a)), the energy released per surface increment by the 3D and front models are close, implying that 2D front models are valid options for describing IERR in this loading configuration. The IERR obtained by the side model is larger than that of the 3D model. With a similar remote stress (Figure 10(b)), the side model more accurately coincides with the behaviors of both 3D and 2D front models. However, the difference with the front PE model becomes more pronounced, while the front PS model maintains good correspondence with the 3D model. In both loading cases, the IERR shows a monotonic increase with the debonding opening for the 2D models, while the 3D model exhibits a maximum. Actually, the front PS model also shows a maximum for a larger opening. However, neither the front nor the side models are able to accurately describe the maximum position of the 3D curve in terms of debonding opening. This maximum is particularly noteworthy given the assumption of a constant $G_c$. Such a maximum could therefore result in a stable growth of the debonding after initiation, see (Mantič 2009; Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023). However, as mixed modes are taken into account in the calculation of $G_c$, it may induce a maximum in the $G_{inc}/\bar{G}_c$ profile based on the stress state (i.e., $\psi$) acting at the interface.

In this regard, Figures 10(c) and 10(d) illustrate the influence of mode mixity on the IERR to critical ERR ratio. For both 3D and side models, a trend similar to the one depicted for $G_{inc}$ is obtained, suggesting that mode mixity remains relatively constant with changes in debonding size. Such finding is further explored in Figure 11 where the mode mixity at the crack front is displayed as a function of the debonding opening. Both the 3D and side models maintain a nearly constant mode mixity close to 0. This indicates that the interface predominantly experienced mode I loading, as described in Equation (6). Consequently, $G_c$ maintains an almost constant value throughout the debonding propagation, ensuring that the trends in IERR remain unaffected.
In contrast, the front models show an increase in mode mixity as a function of the debonding opening. The migration of the crack tip from the fiber pole to the fiber equator leads to a shift from tensile-induced propagation to shear-induced propagation of the debonding. Consequently, the value of $G_c$ rises as mode mixity increases. This increase in $G_c$ consequently influences the trends observed in $G_{inc}/\bar{G}_c$ (Figures 10(c) and 10(d)), which now emphasize a different maximum position compared to the IERR alone. This suggests that the front model might undergo stable debonding growth following initiation, but for a larger opening. However, even when considering variations in $\lambda$, the front models fail to accurately predict the 3D model maximum prediction. Notably, varying $\lambda$ in the 3D model does not influence the opening position determined by the maximum in $G_{inc}/\bar{G}_c$, as discussed in Section 5.2.

Only the side model does not exhibit a maximum in $G_{inc}/\bar{G}_c$. Once the debonding has initiated, crack propagation can be assessed by calculating the ERR ($G$) which can be derived from the IERR using

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

Consequently, in the side PE model where debonding is predominantly driven by mode I, resulting in a constant $G_c$, the IERR shows a monotonic increase. According to Equation (8) the ERR also follows a monotonically increasing trend, along with the ERR to critical ERR ratio. It means that unstable propagation takes place after crack initiation, and no arrest length can be deduced from the side model. This model is unsuitable for performing inverse identification based on arrest length and is therefore disregarded in the following.

<figure><img src="image_13.png" /><figcaption>Figure 11 Mode mixity at the crack tip for front PE, side PE and 3D models.</figcaption></figure>

## 3.4 Comparison of the CC solution

Figure 12 shows the CC solution obtained from the two front models and the 3D model for similar fixed interface fracture properties.

Both strategies devoted to the required displacements or the remote stresses to satisfy the CC solution are investigated. As a result, the 2D model that best matches the 3D model and the corresponding imposed conditions can be defined. As expected, increasing the strength or the critical ERR leads to an increase in the required displacement and remote stress required to fulfill the CC. For sufficiently low strengths, an insensitivity with respect to $\sigma_c$ is obtained, describing an energy-driven configuration where only the critical ERR drives the CC solution. On the one hand, the remote stress required to fulfill the CC is similar for the two 2D front models, where the solutions closely align (Figure 12(b)). However, the 2D front models slightly overestimate the solution of the 3D model. On the other hand, the displacement necessary to satisfy the CC for the two 2D models exhibits a more significant difference, with the disparity being more pronounced in the plane stress assumption, Figure 12(a). Notably, the difference with the 3D model is more prominent in this scenario. The plane strain model provides a more accurate representation of the 3D model solution, especially in cases involving small strengths and configurations with large critical ERRs. This observation is supported by the similarity in the
<figure><img src="image_14.png" /><figcaption>**Figure 12** Required (a) displacement and (b) remote stress to fulfill the CC as a function of a range of strengths and critical ERRs for $\mu = 1.6$.</figcaption></figure>

maximum values of $G_{inc}/G_c$ for the same imposed displacement, as depicted in Figure 10(c). These results highlight the significance of the 2D front PE model in accurately describing the 3D model solution, particularly when considering cases with imposed displacement. However, for the sake of consistency with the 2D identification performed in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023), remote stress solutions are compared in the sequel, where experimental measurements are used to identify interface fracture properties.

## 4 Experimental results

Dog bone samples with a single incorporated fiber are manufactured following the geometry shown in Figure 1. A full description of the manufacturing process can be found in (Livingston and Koohbor 2022). Uniaxial quasi-static tensile tests are carried out under constant cross-head speed at room temperature. Spray paint is applied on the front surface to enable DIC measurements. Debonding initiation and propagation are monitored in situ using a high-magnification camera synchronized with the load cell.

Figure 13(a) shows the $U_y$ displacement field with removed rigid body motion obtained via 2D DIC at the fiber vicinity, superimposed with the associated RAW image where the fiber location can be identified in-between the two displacement discontinuities by the black circle.

<figure><img src="image_15.png" /><figcaption>**Figure 13** Experimental characterization of the debonding process using the displacement jump across the fiber-matrix interface measured by DIC. (a) Displacement field $U_y$ obtained at the fiber vicinity superimposed with the RAW image and the mesh used for the DIC. Locations of the gauge and the fiber are indicated by the straight black line and the black circle, respectively. (b) Debonding opening ($\delta_{nn}$) extracted experimentally as a function of the remote loading ($\sigma^\infty$). (c) Displacement along the gauge for several loading levels.</figcaption></figure>
DIC is assessed using UFreckles (Réthoré 2018) with a converged mesh size of 35 pixels. The pixel size corresponds to 10 µm. The mesh can be visualized in Figure 13(a). Different elements and mesh topologies are also compared, such as mesh topology along the interface or the use of two separate meshes in the fiber and the matrix, but minor differences ranging from 0.2 to 1.5 µm are obtained on the displacement field depending on the remote stress. A 2.5 mm long virtual gauge aligned with the loading direction and passing through the fiber and the matrix at the top pole, as seen in Figure 13(a), is used to determine the $U_y$ displacement as a function of the applied remote load within both the fiber and matrix. The gauge size was chosen to embed several points of measurements (see the mesh in Figure 13(a)). Figure 13(c) shows the variation in displacement along the gauge for several remote loadings. The profile obtained for $\sigma^\infty = 23.3$ MPa is depicted with colors consistent with the associated displacement field shown in Figure 13(a). For sufficiently small loading, i.e., $\sigma^\infty < 4.2$ MPa, the displacement exhibits slight variations along the virtual gauge without any discontinuities at the fiber-matrix interface. Obviously, the displacement presents larger variations within the softer matrix compared to the stiffer fiber. From a certain loading level, between $\sigma^\infty = 4.2$ and $\sigma^\infty = 5.4$ MPa, displacement discontinuously develops at the fiber-matrix interface, which can be associated with the debonding initiation. From this loading level, the displacement remains almost constant within either the matrix or the fiber, allowing the debonding opening $\delta_{nn}$ to be assessed by computing the displacement jump, as depicted in Figure 13(c) with $\delta_{nn}$. It is noteworthy that the displacements are averaged in both phases to compute the debonding opening and the two extremes are used to compute the measurement uncertainties. This allows the remote loading at initiation and the corresponding debonding opening to be extracted. Therefore, opening variation can be assessed from this specific remote loading, see Figure 13(b). The variation in opening is similar to that obtained in (Girard, Koohbor, et al. 2023) measured on the same system under similar condition. Debonding initiation thus occurs at a remote stress ranging from 4.2 to 5.4 MPa. The obtained initiation opening results are summarized in Table 3.

<table><caption>**Table 3** Experimental observation of the initiation debonding opening using DIC measurements and a virtual gauge placed at fiber top pole.</caption><thead><tr><th></th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>Remote stress [MPa]</td><td>4.2</td><td>5.4</td></tr><tr><td>Debonding opening [µm]</td><td>2.1 ± 1.0</td><td>3.8 ± 0.8</td></tr></tbody></table>

Since debonding occurs between two optical acquisitions, the two opening values are used as the lower and upper bounds to evaluate the measurement uncertainties. In the subsequent sections, a two-step methodology is presented to deduce a range of appropriate interface properties. More specifically, inverse property identification is performed by numerically identifying properties to match both the remote stress and the debonding opening after initiation.

# 5 Inverse identification

## 5.1 Identification based on remote stress at initiation

Since linear elastic properties are considered for the fiber and the matrix, FE calculation with several debonding configurations at a single imposed displacement are required to solve the CC for any interface fracture properties. Consequently, CC becomes a highly efficient tool for inverse identification. The inverse identification process is similar to the approach used in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023).

The first step in the inverse identification is to identify a range of properties (i.e., $\sigma_c$, $\tau_c$, $G_{IC}$ and $G_{IIC}$) leading to the debonding initiation remote stress between 4.2 and 5.4 MPa. The initiation stress level observed slightly differs from the one obtained in the previous study (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023). This can be explained by the fact that no relevant debonding opening can be extracted from this sample so that a second similar one is used here. Although debonding could have occurred between 4.2 and 5.4 MPa, the smaller value is chosen as the target to be consistent with that observed using the
debonding angle in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023), i.e., 4.1 MPa.

Figure 14(a) shows range of $G_{IC}$ identified for a range of strengths leading to the initiation of debonding at a remote stress of 4.2 MPa. For sufficiently small strengths (< 6 MPa), a plateau is reached where constant $G_{IC}$ value is identified. In such a configuration, the energy condition drives the initiation where the solution is described by the maximum in the IERR to critical ERR ratio. As a result, tensile and shear strengths no longer influence the CC solution, since the stress condition is already fulfilled; constant critical ERRs are therefore sufficient to match the experimental remote stress. Conversely, for larger strengths, both stress and energy conditions drive debonding initiation. As tensile and shear strengths increase, decreasing critical ERRs are required to counterbalance the increasing stress condition. Figure 14(b) shows the $G_{IC}$ plateau values identified as a function of $\lambda$ where the latter depends on the $G_{IIIC}$ to $G_{IC}$ ratio. It shows that $\lambda$ has a negligible influence on the range of identified $G_{IC}$ values, given that the interface mainly experiences mode I loading, whereas $\lambda$ mostly denotes mode II influence. In fact, each range of identified $G_{IC}$ values, based on a fixed $\lambda$, enables the identification of the corresponding $G_{IIIC}$.

<figure><img src="image_16.png" /><figcaption>**Figure 14** Inverse identification of a (a) range of critical opening ERR as a function of both tensile and shear strengths leading to a debonding initiation at a remote stress of 4.2 MPa and (b) maximum critical ERR in mode I identified as a function of $\lambda$, corresponding to the plateau value.</figcaption></figure>

## 5.2 Properties constraint based on debonding opening at initiation

Once a range of properties is deduced from debonding initiation, a second step involving additional constraints can be undertaken. The arrest opening given by the CC solution is compared with the experimental debonding opening obtained after initiation (see Table 3), enabling the determination of upper and lower opening limits. Figure 15 shows the debonding opening $\delta_{nn}^{\max}$ corresponding to either the crack length or surface maximizing $G_{inc}/\bar{G}_c$ as a function of $\lambda$.

Since the $\bar{G}_c$ calculation relies on the mode mixity, the maximum location remains constant for the 3D due to the dominance of mode I configuration at the interface. As a result, the shear mode has a negligible influence on the calculation of $\bar{G}_c$, and consequently, the maximum location will correspond to that of $G_{inc}$. On the contrary, the front model undergoes a transition from mode I to mode II as the opening increases. The smaller the value of $\lambda$, the larger the value of $G_{IIIC}$, leading to a smaller arrest opening and making debonding propagation more challenging.

With a fixed remote stress of 4.2 MPa, which corresponds to the experimental remote stress observed at initiation, the corresponding debonding opening can be deduced for both models. When debonding initiation is solely driven by the energy condition, stable debonding occurs, resulting in the smallest possible arrest opening, which corresponds to the location of the maximum $G_{inc}/\bar{G}_c$. However, when both stress and energy conditions drive the initiation, larger arrest openings are observed. As a result, the openings shown in Figure 15 correspond to the minimum arrest openings that can be obtained with both models. The 3D model shows a good correspondence with experimental opening whereas the front model overestimates
<figure><img src="image_17.png" /><figcaption>**Figure 15** Debonding opening corresponding to the $G_{inc}/\bar{G}_c$ maximum for both front and 3D models as a function of $\lambda$ and compared with the experimental upper and lower bounds for $\sigma^\infty = 4.2$ MPa.</figcaption></figure>

the results for $\lambda$ larger than 0.6 since the 2D model overestimates $\delta_{nn}$ for a given angle $\theta_d$, see Figure 8(b). Therefore, the front model is able to perform inverse identification based on $\delta_{nn}$ only for small $\lambda$; such a configuration would imply large $G_{IIC}$ to restrain debonding opening at initiation since the arrest opening is primarily driven by shear in 2D. As a consequence, $\lambda = 0.13$ could still be identified in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023) based on the debonding angle which also results in debonding opening consistent with experiments. The 3D model does not allow for the identification of $\lambda$, as the latter does not influence the arrest opening. A meaningful identification of $G_{IIC}$ is no longer feasible using the 3D model and based on $\delta_{nn}$ solely. Possible identification of $G_{IIC}$ could be achieved based on debonding propagation located at the specimen center, where the 2D PE model is applicable. However, this would imply observing the debonding propagation on the fiber equator. Nonetheless, additional constraints could be applied to $G_{IIC}$ based on arrest openings obtained from further unstable growth.

Figure 16 shows the arrest debonding opening as a function of a range of strengths identified during the first identification step. Similarly to the first step, when dealing with small strengths, an energy-driven configuration prevails for stable debonding, i.e., $G = G_c$ at initiation. Hence, there is no additional unstable propagation of debonding, with the initiation opening depending only on constant ERRs over such a range of strengths. Consequently, constant arrest openings, equivalent to initiation openings, are obtained, leading to the formation of a plateau. With larger strengths, there is a possibility of unstable propagation of the debonding, resulting in additional opening growth until $G < G_c$. These configurations therefore induce a larger arrest opening. By comparing the numerical arrest opening with the experimental bounds, the range of identified strengths can be validated since the arrest openings remain within the two experimental bounds. No additional restrictions can be applied to the range of strengths. Similarly, since the complete range of strengths matches the experimental observations, the identified range of critical ERR can also be approved.

### 5.3 Comparison with 2D properties identification

Table 4 compares the properties identified with the 3D model with the ones obtained with the 2D Front PE model in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023) to quantify the relevance of a 2D model in accurately identifying the interface properties. Even though, the stress components experienced by the interface in both models are different, the 2D model is consistent with the 3D obtained properties regarding the critical ERR in mode I. Similar properties are obtained, with the 3D ones encompassing the 2D ones. The 2D model slightly underestimates the tensile strengths as it does not account for the singularity at the free edge. However, due to the fact that the 3D model involves less shear at the interface compared to the 2D model, differences in shear properties are likely and raise questions about the relevance of such identification. On the one hand, the choice of an equivalent stress in 3D might be justifiable since it closely aligns with a normal stress criterion. The shear strengths deduced from this stress
<figure><img src="image_18.png" /><figcaption>**Figure 16** Arrest debonding opening as a function of range of strengths with lower and upper bound obtained from the experimental observation for $\lambda = 0.13$.</figcaption></figure>

<table><caption>**Table 4** Comparison of the interface properties obtained by inverse identification from the 3D model and previously in (Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023) using a 2D front model under plane strain assumption based on the debonding angle instead of the debonding opening.</caption><thead><tr><th>Interface properties</th><th>3D</th><th>2D front PE</th></tr></thead><tbody><tr><td>σ<sub>c</sub> [MPa]</td><td>4 - 11</td><td>1 - 5</td></tr><tr><td>τ<sub>c</sub> [MPa]</td><td>/</td><td>1 - 10</td></tr><tr><td>G<sub>IC</sub> [N/mm]</td><td>0.0022 - 0.0047</td><td>0.0037</td></tr><tr><td>G<sub>IIIC</sub> [N/mm]</td><td>/</td><td>0.09</td></tr></tbody></table>

condition might be erroneous because regardless of the chosen properties, their influence on the stress condition is negligible, so that comparison between the two models is thus no longer feasible. On the other hand, the identification of $\lambda$ is challenging since it does not drive the arrest opening. Consequently, it is not possible to identify any critical ERR in mode II due to the unknown value of $\lambda$. It would require analyzing further debonding propagation with larger debonding angle in the specimen middle plane where shear has larger influence. It is worth mentioning that debonding shapes are based on stress isocontours while the energy isocontours could also be suitable. Such energy isocontours may yield debonding shapes experiencing more shear and allow for a better capture of both $\tau_c$ and $G_{IIIC}$. Furthermore, these new debonding shapes may better describe debonding angles observed experimentally with a variation of $\delta_{nn}$. However, the determination of energy isocontours is challenging to implement and will be developed in subsequent studies.

# 6 Conclusion

Initiation of fiber-matrix debonding was assessed using 3D and 2D models and the CC on a single transverse fiber configuration. The different concluding remarks are summarized below:

1. Crack path determination: shear mode has a negligible influence over the potential crack path so that deriving debonding shapes based on normal stress isocontours is a relevant solution since it does not rely on any fracture properties;

2. Model comparison: 2D front models are unable to capture stress singularities at the free edge, unlike 2D side model. The IERR to critical ERR ratio of the side model exhibits a monotonic increase, rendering this model unsuitable for predicting a crack arrest after initiation and consequently inappropriate for inverse identification based on the crack arrest configuration. When comparing the CC solution, agreement is met between the 3D solution and the 2D front PE solution with a similar imposed displacement;

3. Experimental results: debonding initiates at a remote stress ranging between 4.2 and 5.4 MPa for an opening comprised between 2.1 and 3.8 µm;

4. Inverse identification: the properties identified using the 2D front model in (Girard, Doitrand,
Koohbor, Rinaldi, Godin, Long, Bikard, and Trouillet-Fonti 2023) align with those in 3D concerning tensile strength and mode I critical ERR. This makes the front PE model relevant for efficiency purposes in identifying the interface fracture properties even though the 2D front model yields a slight underestimation of the tensile strength. The shear strength can no longer be identified using the 3D model since crack initiation mainly occurs under opening mode, unlike in 2D front models where both opening and shear modes are involved. Similarly, it is not possible to identify the critical shear ERR because the influence of $\lambda$ is negligible.

Overall, the numerical approach can be extended to more realistic composite configurations, including those with fiber diameters of just a few microns. However, for smaller fibers, a scanning electron microscope would be necessary to accurately characterize the debonding process. In addition, neighboring fibers could also be considered to mimic the real microstructure of composites. They would influence the crack path and the intensity of the stress and energy fields, as observed in (Girard et al. 2024). As a perspective, new debonding paths could be incorporated into the CC implementation to account for energy considerations. Energy-based isocontours might offer a more consistent representation of the observed debonding angle at the free edge, and provide a more pronounced shear influence during debonding initiation. Additionally, potential unstable debonding propagation after initiation is assessed based on stress isocontours, the optimal debonding shape could be determined by identifying the debonding shape that minimizes the imposed loading, among all possible initiation debonding shapes.

## References

Ageorges, C., K. Friedrich, T. Schüller, and B. Lauke (1999). Single-fibre Broutman test: fibre-matrix interface transverse debonding. *Composites Part A: Applied Science and Manufacturing 30*(12):1423-1434. [DOI].

AhmadvashAghbash, S., I. Verpoest, Y. Swolfs, and M. Mehdikhani (2023). Methods and models for fibre-matrix interface characterisation in fibre-reinforced polymers: a review. *International Materials Reviews* 68(8):1245-1319. [DOI].

Broutman, L. J. (1966). Glass-resin joint strength and their effect on failure mechanisms in reinforced plastics. *Polymer Engineering & Science* 6(3):263-272. [DOI].

Carrere, N., A. Doitrand, E. Martin, and D. Leguillon (2021). 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:102979. [DOI], [OA].

Dève, H. E. and S. Schmauder (1992). Role of interface properties on the toughness of brittle matrix composites reinforced with ductile fibers. *Journal of Materials Research* 7(11):3132-3138. [DOI].

Doitrand, A., T. Duminy, H. Girard, and X. Chen (2024). A review of the coupled criterion. *Journal of Theoretical, Computational and Applied Mechanics*, 1-37. [DOI], [OA].

Doitrand, A., C. Fagiano, N. Carrère, V. Chiaruttini, and M. Hirsekorn (2017). Damage onset modeling in woven composites based on a coupled stress and energy criterion. *Engineering Fracture Mechanics* 169:189-200. [DOI], [HAL].

Doitrand, A. and D. Leguillon (2018a). 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:175-182. [DOI], [OA].

Doitrand, A. and D. Leguillon (2018b). 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:69-76. [DOI], [HAL].

Duminy, T., R. Henry, J. Adrien, A. Doitrand, and S. Meille (2023). Anisotropic fracture in nacre-like alumina. *Theoretical and Applied Fracture Mechanics* 123:103710. [DOI], [OA].

García, I., B. Carter, A. Ingraffea, and V. Mantič (2016). A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics. *Composites Part B: Engineering* 95:475-487. [DOI].

García, I., V. Mantič, and E. Graciani (2015). Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds? *European Journal of Mechanics - A/Solids* 53:75-88. [DOI].
García, I., M. Paggi, and V. Mantič (2014). 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:96-110. [DOI].

Gentieu, T., J. Jumel, A. Catapano, and J. Broughton (2019). Size effect in particle debonding: Comparisons between finite fracture mechanics and cohesive zone model. *Journal of Composite Materials* 53(14):1941-1954. [DOI], [OA].

Girard, H., A. Doitrand, B. Koohbor, R. Rinaldi, N. Godin, D. Long, and J. Bikard (2024). Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination. *Journal of the Mechanics and Physics of Solids* 183:105498. [DOI], [OA].

Girard, H., A. Doitrand, B. Koohbor, R. Rinaldi, N. Godin, D. Long, J. Bikard, and L. Trouillet-Fonti (2023). Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties. *Engineering Fracture Mechanics* 286:109254. [DOI], [OA].

Girard, H., B. Koohbor, A. Doitrand, and R. Livingston (2023). Experimental characterization of in-plane debonding at fiber-matrix interface using single glass macro fiber samples. *Composites Part A: Applied Science and Manufacturing* 171:107573. [DOI], [HAL].

Gundel, D., B. Majumdar, and D. Miracle (1995). Evaluation of the transverse response of fiber-reinforced composites using a cross-shaped sample geometry. *Scripta Metallurgica et Materialia* 33(12):2057-2065. [DOI].

Hutchinson, J. and Z. Suo (1991). Mixed Mode Cracking in Layered Materials. *Advances in Applied Mechanics*. Ed. by J. W. Hutchinson and T. Y. Wu. Vol. 29. Elsevier, pp 63-191. [DOI].

Johnson, A. C., S. A. Hayes, and F. R. Jones (2012). The role of matrix cracks and fibre/matrix debonding on the stress transfer between fibre and matrix in a single fibre fragmentation test. *Composites Part A: Applied Science and Manufacturing* 43(1):65-72.

Koyanagi, J., P. D. Shah, S. Kimura, S. K. Ha, and H. Kawada (2009). Mixed-mode interfacial debonding simulation in single-fiber composite under a transverse load. *Journal of Solid Mechanics and Materials Engineering* 3(5):796-806. [DOI].

Kushch, V., S. Shmegera, P. Brøndsted, and L. Mishnaevsky (2011). Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading. *International Journal of Engineering Science* 49(1):17-29. [DOI].

Lee, D., J.-H. Kim, S. J. Lee, M. Kim, and D.-J. Kwon (2024). Microdroplet pull-out testing: Significance of fiber fracture results. *Polymer Testing* 141:108631. [DOI], [OA].

Leguillon, D. (2014). An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case. *Theoretical and Applied Fracture Mechanics* 74:7-17. [DOI], [HAL].

Leguillon, D. (2002). Strength or toughness? A criterion for crack onset at a notch. *European Journal of Mechanics - A/Solids* 21(1):61-72. [DOI], [HAL].

Livingston, R. and B. Koohbor (2022). Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements. *Composites Part C: Open Access* 7:100229. [DOI], [OA].

Mandell, J., J. Chen, and F. McGarry (1980). A microdebonding test for in situ assessment of fibre/matrix bond strength in composite materials. *International Journal of Adhesion and Adhesives* 1(1):40-44. [DOI].

Mantič, V. (2009). 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):1287-1304. [DOI], [OA].

Martin, E., D. Leguillon, A. Catapano, and N. Carrère (2020). Prediction of interfacial debonding between stiff spherical particles and a soft matrix with the coupled criterion. *Theoretical and Applied Fracture Mechanics* 109:102749. [DOI], [OA].

Martyniuk, K., B. F. Sørensen, P. Modregger, and E. M. Lauridsen (2013). 3D in situ observations of glass fibre/matrix interfacial debonding. *Composites Part A: Applied Science and Manufacturing* 55:63-73. [DOI].

Meer, F. P. van der, S. Raijmaekers, and I. B. Rocha (2019). Interpreting the single fiber fragmentation test with numerical simulations. *Composites Part A: Applied Science and Manufacturing* 118:259-266. [DOI], [OA].

Meurs, P., B. Schrauwen, P. Schreurs, and T. Peijs (1998). Determination of the interfacial normal strength using single fibre model composites. *Composites Part A: Applied Science and Manufacturing* 29(9-10):1027-1034. [DOI].
Muñoz-Reja, M., L. Távara, V. Mantič, and P. Cornetti (2016). 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. *Procedia Structural Integrity 2*:2022-2029. [DOI], [OA].

Ogihara, S., Y. Sakamoto, and J. Koyanagi (2009). Evaluation of Interfacial Tensile Strength in Glass Fiber/Epoxy Resin Interface using the Cruciform Specimen Method. *Journal of Solid Mechanics and Materials Engineering* 3(9):1071-1080. [DOI], [OA].

Park, S.-J., M.-K. Seo, and J.-R. Lee (2006). Roles of interfaces between carbon fibers and epoxy matrix on interlaminar fracture toughness of composites. *Composite Interfaces* 13(2-3):249-267. [DOI].

Pitkethly, M., J. Favre, U. Gaur, J. Jakubowski, S. Mudrich, D. Caldwell, L. Drzal, M. Nardin, H. Wagner, L. Di Landro, A. Hampe, J. Armistead, M. Desaeger, and I. Verpoest (1993). A round-robin programme on interfacial test methods. *Composites Science and Technology* 48(1-4):205-214. [DOI].

Réthoré, J. (2018). *UFreckles*. Version 2.0. [DOI].

Sandino, C., E. Correa, and F. París (2016). Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load. *Engineering Fracture Mechanics* 168:58-75. [DOI], [OA].

Damage in composite materials (2012). *Damage and Failure of Composite Materials*. Ed. by C. V. Singh and R. Talreja. Cambridge University Press, pp 36-56. [DOI].

Tandon, G. P. (2000). Evaluation of interfacial normal strength in a SCS-0/epoxy composite with cruciform specimens. *Composites Science and Technology* 60(12-13):2281-2295. [DOI].

Tandon, G. P., R. Y. Kim, and V. T. Bechel (2002). Fiber-Matrix Interfacial Failure Characterization Using a Cruciform-Shaped Specimen. *Journal of Composite Materials* 36(23):2667-2691. [DOI].

Totten, K. R., B. Kutub, and L. A. Carlsson (2016). In situ determination of the fiber-matrix interface tensile strength. *Journal of Composite Materials* 50(5):589-599. [DOI].

Tripathi, D. and F. R. Jones (1998). Single fibre fragmentation test for assessing adhesion in fibre reinforced composites. *Journal of Materials Science* 33(1):1-16. [DOI].

Velasco, M., E. Correa, and F. París (2020). Interaction between fibres in the transverse damage in composites. *Engineering Fracture Mechanics* 239:107273. [DOI].

Zarges, J.-C., C. Kaufhold, M. Feldmann, and H.-P. Heim (2018). Single fiber pull-out test of regenerated cellulose fibers in polypropylene: An energetic evaluation. *Composites Part A: Applied Science and Manufacturing* 105:19-27. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** H. Girard: conceptualization, methodology, writing, investigation; A. Doitrand: conceptualization, methodology, writing, investigation; B. Koohbor: investigation, validation, resources; R. Rinaldi: conceptualization, methodology, writing, investigation; N. Godin: supervision, writing, validation; J. Bikard: funding acquisition, supervision

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.