# Fiber-matrix interface debonding and transverse cracking in macro fiber composites

Behrad Koohbor, Zaynab Hazaveh, Aurélien Doitrand, Hugo Girard

► To cite this version:

Behrad Koohbor, Zaynab Hazaveh, Aurélien Doitrand, Hugo Girard. Fiber-matrix interface debonding and transverse cracking in macro fiber composites. Engineering Fracture Mechanics, 2025, 325, pp.111345. 10.1016/j.engfracmech.2025.111345. hal-05115375

HAL Id: hal-05115375

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

Submitted on 16 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fiber-Matrix Interface Debonding and Transverse Cracking in Macro Fiber Composites

Behrad Koohbor<sup>a,c,*</sup>, Zaynab Hazaveh<sup>a</sup>, Aurélien Doitrand<sup>b</sup>, Hugo Girard<sup>b</sup>

a. Department of Mechanical Engineering, Rowan University, 201 Mullica Hill Rd., Glassboro, NJ 08028 USA

b. Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, Villeurbanne, UMR5510, 69621, France

c. Advanced Materials and Manufacturing Institute (AMMI), Rowan University, Glassboro, NJ, 08028 USA

* Corresponding Author, Email: koohbor@rowan.edu, Phone: +1-856-256-5328

## Abstract

Fiber-matrix interface debonding is a precursor to transverse matrix cracking at the mesoscopic scale in fiber composites. The mechanisms controlling fiber-matrix interface debonding and subsequent transverse crack formation have been explored primarily by computational methods with limited experimental verification. This study aims to establish an experimental approach for characterizing fiber-matrix interface debonding in model macro fiber specimens that replicate realistic microstructures. The primary goal is to measure strain fields at individual fiber-matrix interfaces using optical digital image correlation (DIC) and link these measurements to the initiation and propagation of transverse cracks. Macro fiber composite specimens are fabricated by embedding dozens of randomly distributed glass macro fibers (1 mm dia.) in an epoxy matrix. These specimens are then subjected to controlled transverse loading, and their local strain fields are monitored and quantified with high-magnification optical DIC. The experimentally obtained kinematic fields are first used to connect global and local deformation responses and to investigate the mechanisms governing matrix cracking between the fibers. The experimental data are then used to set up and validate a modeling framework created based on cohesive zone and phase field formulations to investigate fiber-matrix interface debond initiation and matrix cracking, respectively. The experimental protocols described here provide a practical approach for characterizing deformation and failure at the fiber-matrix interface and tracking their evolution into larger transverse cracks. Complementary simulation studies highlight the significance of boundary conditions and the uncertainty in the fiber-matrix interface fracture properties in realistic and reliable predictions of debonding kinetics and matrix crack formation. The presented approach is transferable to smaller length scales to enable the quantitative assessment of the effects of geometric and morphological factors, such as inter-fiber distance and angle, on transverse crack formation in fiber composites.

**Keywords:** Fiber-matrix debonding; Damage mechanics; Transverse cracking; Digital image correlation
<table><thead><tr><th colspan="2">Nomenclature</th></tr></thead><tbody><tr><td>CT</td><td>Computed Tomography</td></tr><tr><td>DIC</td><td>Digital Image Correlation</td></tr><tr><td>ERR</td><td>Energy Release Rate</td></tr><tr><td>RVE</td><td>Representative Volume Element</td></tr><tr><td>SEM</td><td>Scanning Electron Microscopy</td></tr><tr><td>d</td><td>Phase Field Damage Variable</td></tr><tr><td>E</td><td>Elastic Modulus</td></tr><tr><td>E<sub>f</sub></td><td>Elastic Modulus of Fiber</td></tr><tr><td>E<sub>m</sub></td><td>Elastic Modulus of Matrix</td></tr><tr><td>G<sub>c</sub></td><td>Critical Energy Release Rate for Matrix</td></tr><tr><td>G<sub>Ic</sub></td><td>Critical Energy Release Rate in Mode I Fracture</td></tr><tr><td>G<sub>IIc</sub></td><td>Critical Energy Release Rate in Mode II Fracture</td></tr><tr><td>L<sub>sub</sub></td><td>DIC Subset size</td></tr><tr><td>L<sub>st</sub></td><td>DIC Subset spacing,</td></tr><tr><td>L<sub>window</sub></td><td>DIC strain window size</td></tr><tr><td>L<sub>VSG</sub></td><td>Virtual strain gauge size</td></tr><tr><td>n<sub>d</sub></td><td>Number of Debonded Fibers</td></tr><tr><td>n<sub>f</sub></td><td>Total Number of Fibers</td></tr><tr><td>δ<sub>n</sub></td><td>Fiber-Matrix Interface Opening Displacement</td></tr><tr><td>ε</td><td>Global Strain</td></tr><tr><td>ε<sub>yy</sub></td><td>In-plane Axial Strain</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>σ</td><td>Remote (Global) Tensile Stress</td></tr><tr><td>σ<sub>c</sub></td><td>Interface Normal Strength</td></tr><tr><td>σ<sub>f</sub></td><td>Matrix Tensile Strength</td></tr><tr><td>τ<sub>c</sub></td><td>Interface Shear Strength</td></tr></tbody></table>

# 1. Introduction

Most damage types in composites initiate at nano or sub-micrometric length scales, rapidly developing into larger macroscale failure modes. An example of such multiscale failure mechanisms is transverse and matrix cracking in unidirectional and cross-ply composites. These damage types initiate at the fiber-matrix interface. An interface debond then propagates
circumferentially around the fiber, eventually kinking into the matrix. The coalescence of matrix microcracks between neighboring fibers forms larger cracks that propagate through the lamina thickness, significantly reducing the mechanical load-bearing capacity of the lamina (e.g., stiffness reduction) and leading to other damage types like delamination [1, 2]. Because *in-situ* and real-time visualization and characterization of this type of damage are extremely difficult—if not impossible—most research works on fiber-matrix debonding and transverse cracking rely on numerical modeling with limited experimental validation or *ex-situ* experiments where material behavior is observed and quantified *before and after* failure, without capturing the damage progression *during* the process.

Modeling studies on transverse cracking mechanisms in composites vastly outnumber experimental reports. Among these studies, Guillén-Hernández et al. [3] developed a combined phase-field and cohesive zone model to study the failure of a single fiber system, focusing on fiber-matrix debonding and crack kinking. Simulations were then used to investigate the influence of neighboring fibers and crack coalescence between debonded fibers. In another work, Guillén-Hernández et al. [4] employed similar strategies to model fiber-matrix debonding and matrix cracking in realistic RVEs representing actual composite microstructures. Despite using realistic microstructures and effective modeling of early-stage damage, the conclusiveness of their findings was challenged by idealized material properties and boundary conditions, making this modeling framework only qualitatively applicable to practical conditions. Combined phase-field and cohesive zone models have also been used to predict complex crack paths in ceramic-matrix composites [5]. Several other notable studies have developed micromechanics models to examine transverse cracking in cross-ply composites [6] and thin-ply laminates [7], as well as UD composites under compression [8].
Unlike the abundance of modeling studies, experimental works on fiber-matrix debonding and transverse cracking are scarce. *In-situ* experimental studies that enable characterization of interface debonding and inter-fiber cracking have only become feasible in recent years owing to advancements in electron microscopy and non-contact full-field measurements, such as digital image correlation (DIC) [9]. A pioneering study by Tracy *et al.* [10] examined multiscale damage in continuous fiber composites, from single-fiber to laminate scales. Full-field strain data obtained and analyzed in this work emphasized the need for multiscale analyses to bridge damage initiation at the single-fiber level to laminate-scale failure. Mehdikhani *et al.* [11] combined in-SEM testing with DIC to validate finite element models at single-fiber and RVE scales, showing that the detection of small-scale strain concentrations requires enhanced DIC patterns and optimized image correlation parameters. To address the challenges pointed out in these works, Montgomery *et al.* [12] developed a new speckle pattern production technique that proved effective in measuring strain fields near single glass fibers tested inside an SEM chamber. This method was later applied to studying fiber-matrix debonding and matrix cracking in composites with realistic microstructures [13]. Despite the valuable information gained from these SEM-based *in-situ* measurements, challenges remain with testing polymer composite specimens inside an SEM chamber. Issues such as the need for a well-dispersed nanometric speckle pattern, the intermittent loading-imaging process, and image distortion due to beam-matter interactions (e.g., charging) can negatively impact the accuracy of these measurements [14]. Some of these challenges also affect alternative X-ray Computed Tomography (CT) approaches used for measuring fiber-matrix interface failure [15].

To address some of the aforementioned limitations, the application of macro fibers (100 to 1000 µm in diameter) combined with optical measurement methods has gained considerable attention.
Combining large-diameter fibers with optical techniques eliminates the need for intermittent
loading, allowing continuous data collection without the potential complications caused by high-
energy electron or X-ray beams. This approach also simplifies the application of DIC-suitable
patterns. Studies by Tabiai *et al.* [16-18] demonstrated the effectiveness of macro fiber composite
testing with optical full-field measurements in revealing mechanisms of fiber-matrix interface
debonding. These studies also suggested that damage behavior in fiber bundles mirrors that of
isolated fiber-matrix interfaces in single macro fiber specimens, showing similar debonding
initiation and propagation mechanisms. Furthermore, studies conducted on macro fibers with
diameters ranging from a few hundred micrometers to millimeters were shown to provide results
that were independent of fiber dimension, especially when using DIC to extract local stain
information from regions in the matrix that are far enough from the fiber-matrix interface [16, 17].
Recent works by Nikforooz *et al.* [19] and Verschatse *et al.* [20] further highlighted the application
of optical methods in the real-time *in-situ* measurement of interface debonding in short-fiber and
single long fiber composites, respectively.

The present authors have developed experimental techniques to further demonstrate the
effectiveness of macro fiber composites in understanding transverse cracking and fiber-matrix
interface failures. Full-field measurements on single and double glass macro fiber specimens were
used to analyze local vs. global strain fields before and during interface debonding [21]. Using
these techniques, strain fields were measured across various scales and utilized to identify interface
fracture properties based on coupled-criterion fracture models [22, 23]. These models were then
employed to correlate matrix crack formation with inter-fiber distance and angles [24]. Attempting
to provide a more realistic three-dimensional view of the interface debonding behavior and to
provide a sufficiently accurate substitute for X-ray CT analyses [15, 25], Uddin *et al.* [26]
developed a dual-vision setup that enabled simultaneous characterization of in-plane debonding and through-thickness crack tunneling in single glass macro fiber specimens. Tunneled crack shape predictions were recently performed by using experimental results extracted from the dual-vision experimental setup that calibrated coupled fracture stress-energy criteria, providing a unique insight into the micromechanics and failure behaviors at the fiber-matrix interface in fiber composites [27].

Despite the recent advancements in understanding the fiber-matrix interface debonding from both experimental and simulation perspectives, there are knowledge gaps on how to bridge the failure response between single fibers and fiber bundles. In particular, the representativeness and reliability of damage mechanism identification in single-fiber specimens and their relevance to actual composites are yet to be fully understood. A recent study by Chatziathanasiou et al. [28] addressed this question through a series of X-ray synchrotron holo-tomography measurements on single carbon and glass fiber specimens as well as their multiple fiber counterparts. Their findings suggest certain mechanistic differences between single- and multi-fiber composites pertaining to interface damage and matrix cracking. While they provided no concrete justification for such differences, Chatziathanasiou et al. [28] hypothesized that the discrepancies between the damage mechanics in single- and multi-fiber specimens could originate from matrix size effects, fiber-matrix interphase formations, fiber volume fraction effects, thermal expansion differences, or energy release upon fracture. The latter hypothesis was justified by estimating the strain energy that was shown to be several times higher in single-fiber composites than their multi-fiber counterparts. This significantly higher strain energy must be dissipated by substantial damage developed in the form of more pronounced matrix cracking and interfacial debonding in single-fiber compared with multi-fiber specimens.
The present study aims to address some of the unanswered questions relevant to the applicability
of macro fiber specimens in studying the transverse cracking behavior of fiber composites.
Specifically, this work provides additional insights into the representativeness of single macro fiber
composites (comprehensively discussed in previous publications) in damage characterization and
modeling in multi-fiber composites. To this end, a hybrid experimental-computational approach is
developed wherein multi glass macro fiber specimens are first fabricated and mechanically tested.
DIC is used to measure the local displacement and strain fields developed in the specimens, also
allowing the determination of interface debonding evolution with global tensile loading. In
parallel, a coupled cohesive zone-phase field modeling framework is developed based on realistic
DIC-measured boundary conditions and interface fracture properties identified previously in our
single-fiber specimen tests. The experimental protocols developed here provide a practical
approach for quantifying the kinetics of damage at the fiber-matrix interface and the progression
of local interface debonding into larger transverse cracks. Complementary finite element
simulations further highlight the roles of boundary conditions and the scatter in interface fracture
properties on transverse cracking behavior in fiber composites.

**2. Experimental Methods**

**2.1. Materials and Specimen Preparation**

Tensile test specimens were fabricated by embedding 39 individual glass rods (hereafter referred
to as macro fibers) in a thermoset epoxy. Quartz glass macro fibers (1 mm diameter,
McMaster Carr, USA) with nominal Young's modulus of 72 GPa were embedded in a translucent
room temperature curable epoxy matrix (105/209 epoxy, West System, MI, US) with the aid of
custom-designed 3D printed jigs. As shown in **Figure 1a**, the jigs were designed with holes slightly larger than 1 mm to allow for free insertion of the macro fibers in a silicone mold. The dimensions of the silicone mold cavity were chosen in accordance with ASTM D638 with a *ca.* 13 mm gauge width. Before insertion into the mold, glass macro fibers were cleaned with isopropanol and spray coated with a silane glass treatment solution (AP115 Silane Glass Treatment, 3M, MN, USA). To do so, the fibers were cut to dimensions and then sprayed with the silane solution. The coated fibers were kept in ambient condition for five minutes, allowing for the solution to settle on the surface. The excess solution was then gently wiped from the surface of the fibers with a clean microfiber cloth. After wiping off the sprayed solution, a hydrophobic surface was created that led to stronger bonds with the cured epoxy resin.

After resin casting, specimens were degassed and then cured at room temperature inside a chemical fume hood for 48 hours. The room temperature curing ensured minimal effects of any potential residual thermal stresses in the specimen. After curing, the specimens were extracted from silicone molds and left in room conditions for an additional 12 hours to ensure full curing on all surfaces and through the thickness. The 3D printed jig was then removed and the extra length of the glass macro fibers outside of the cured specimens was carefully trimmed using a diamond glass cutter followed by gentle polishing with 1000 and 2000 grit sandpapers. The rear surface of the specimen (i.e., the surface in contact with the bottom mold cavity) was also lightly polished with fine sandpapers to ensure that fiber surfaces were flush with the surface of the epoxy matrix. **Figure 1b** illustrates the as-polished surface of a representative specimen, showing fiber locations. The coordinates of the fiber centers were identified using the circle detection tool in the commercially available image processing program ImageJ (NIH, USA). The fiber coordinates were used to characterize inter-fiber distance and angles, as well as the input to create finite element models.
**Figure 1c** and **Figure 1d** show the high-magnification view of the location of the fibers and the corresponding digitized fiber map in a representative specimen, respectively.

After marking the location of the fibers, the polished surface of the specimen was coated with a thin layer of white matte paint and then a black speckle pattern for DIC measurements. The speckle-patterned specimen was affixed in a Shimadzu 10 kN universal test frame and subjected to uniaxial tensile force applied at a constant rate of 1 mm/min. Tensile loading continued until complete failure of the specimens. During mechanical tests, a 5-megapixel camera (Basler Ace USB3, Germany) fitted with a 100 mm macro lens (Tokina AT-X M100 PRO D, Japan) was used to capture high magnification images from the speckle patterned area on the specimen surface. The image acquisition rate was set to 1 Hz, synced with the load data sampling rate. Images captured during testing were analyzed in the commercially available DIC program Vic-2D (Correlated Solutions, Inc., SC, USA) with image correlation parameters optimized for maximum signal-to-noise ratio. Image correlation parameters used in this work were similar to those reported in our previous studies [26] and are listed in **Table 1**.
<figure><img src="image_2.png" /><figcaption>**Figure 1-** (a) Schematic of the specimen fabrication process. (b) An as-prepared specimen with 39 glass macro fibers with a magnified view of the macro fiber cross sections shown in (c). The red markers in (b) and (c) were used to overlay the DIC region of interest on the specimen. (d) The digitized map of macro fibers reproduced using the circle detection tool in ImageJ. Tensile load was applied in y-direction.</figcaption></figure>

<table><caption>**Table 1-** Imaging and digital image correlation parameters</caption><tr><td>Digital magnification (µm/pixels)</td><td>10.5</td></tr><tr><td>Subset size, L<sub>sub</sub> (pixels)</td><td>21×21</td></tr><tr><td>Subset spacing, L<sub>st</sub> (pixels)</td><td>1</td></tr><tr><td>DIC strain window size, L<sub>window</sub> (data points)</td><td>5×5</td></tr><tr><td>Virtual strain gauge size*, L<sub>VSG</sub> (µm)</td><td>262</td></tr><tr><td>Strain metric</td><td>Gaussian Weighted, Lagrangian large strain</td></tr><tr><td>Strain noise floor</td><td>4.8×10<sup>-4</sup></td></tr><tr><td>Strain bias</td><td>7.3×10<sup>-5</sup></td></tr></table>

<footnote>* Determined as L<sub>VSG</sub> = (L<sub>window</sub>-1) L<sub>st</sub> + L<sub>sub</sub> [29]</footnote>
## 2.2. Specimen Types

The specimen fabrication process described above allows the design of macro glass fiber composites with various fiber numbers and distributions. The present report focuses on three separate specimen designs and conditions. First, triplicate tensile test specimens with no glass macro fibers were fabricated using the process described in **Sec. 2.1**. These neat epoxy specimens were used to characterize the stress-strain response of the epoxy resin to be used as input for FE models, as well as to assess the interface debonding-induced progressive damage discussed in later sections.

Two types of tensile test specimens embedding 39 individual glass macro fibers were designed and fabricated. The rationale for using 39 fibers was (1) based on the tensile specimen dimensions that limited the inclusion of higher fiber numbers, and (2) to ensure that the total number of fibers in the area of interest is higher than those previously identified as the minimum numbers required for a statistically equivalent representative volume element in fiber composites. To further elaborate on the latter point, it should be noted that previous studies by Swaminathan *et al.* [30] and Koohbor *et al.* [13] suggest that the smallest volume element that can represent the macroscale response of a fiber composite must incorporate at least 26 randomly distributed fibers. Consistent fiber numbers to represent the statistically representative volume elements in fluid flow in fibrous and porous media were reported by Geoffre *et al.* [31, 32] and Trias *et al.* [33]. Therefore, we chose the abovementioned number to ensure that RVE scale measurements herein represent the macroscopic behavior in a fiber composite.

The first 39-fiber specimen (hereafter referred to as RVE-1) was fabricated considering the actual microstructures reported in Koohbor *et al.* [13]. As shown in **Figure 2a**, the fiber locations revealed by the SEM image of the specimen discussed in Ref. [13] were tracked and used to create its macro
fiber replica (Figure 2b). Note that due to the limitations imposed by the design and fabrication of the 3D printed jigs, fiber packing densities in our macro fiber specimens were significantly lower than those of a realistic composite. For example, the fiber volume fraction of RVE-1 was approximated as 20.5%, determined as the ratio of the total surface of the 39 fibers and the rectangular area encompassing them. RVE-1 was considered the main test specimen in this work for image correlation and debonding study purposes. In addition, experimental data obtained from testing of this specimen were used to calibrate the finite element model replicas.

The second 39-fiber specimen (labeled RVE-2) was fabricated using the same procedure discussed above but with a different distribution of fibers. This specimen was designed primarily as a test condition to verify the accuracy of the finite element model established based on RVE-1 datasets. A magnified photograph of the area of interest in RVE-2 and its digitized map of macro fibers are shown in Figure 2c and Figure 2d, respectively. The fiber volume fraction in RVE-2 was approximated as 17.7%. All specimens described above were made with 7 mm thickness. As discussed above, the experimental results and their associated discussions will focus mainly on RVE-1.
<figure><img src="image_3.png" /><figcaption>**Figure 2-** (a) SEM image of an actual cross-ply composite used as a guideline for designing the macro fiber specimen, RVE-1. (b) Image of the as-prepared specimen RVE-1. (c) Image of the as-prepared specimen RVE-2 with its corresponding digitized map of macro fibers shown in (d). Panel (a) is reproduced from Ref. [13] with permission (Elsevier, 2022).</figcaption></figure>

**2.3. Fiber Distribution in Macro Fiber Specimens**

Inter-fiber spacing and angles were determined for the macro fiber specimens by identifying the coordinates of the fiber centers. The closest neighbors for a specific fiber were first identified visually from the digitized fiber location maps shown earlier in **Figure 1d** and **Figure 2d** for RVE-1 and RVE-2 specimens, respectively. The center-to-center distances for a specific fiber relative to its closest neighbors were then calculated using trigonometric formulae. The actual closest neighbor was then selected based on the shortest inter-fiber (center-to-center) distance. Once known, the angular orientation between pairs of closest neighbors was determined as an angle
between 0 and 90°, as schematically shown in **Figure 3a**. Inter-fiber distances and angles for the two specimens are shown in **Figure 3b-e**. As depicted by the histograms, the inter-fiber distances cover a range of 1 to 2 mm for RVE-1 (specimen with higher fiber volume fraction) and from 1 to 2.8 mm for RVE-2 (specimen with lower fiber volume fraction), noting that a center-to-center fiber distance approaching 1 indicates two fibers in very close proximity of each other. Similarly, the inter-fiber angles for the two specimens cover a range of 0 to 90° for the two specimens with no particular trend. An inter-fiber angle of 0° indicates that the two closest neighbors are aligned perpendicular to the direction of the applied tensile force. Conversely, a 90° angle indicates that the two closest neighbors are aligned with the direction of the tensile load. The fiber coordinates used to quantify fiber distribution were used as input for creating the finite element models, as well.

<figure><img src="image_4.png" /><figcaption>**Figure 3-** (a) Visual description of the inter-fiber distance and angle used in this work. (b) Center-to-center distance and (c) inter-fiber angles between the closest fiber pairs in RVE-1. (d) Center-to-center distance and (e) inter-fiber angles between the closest fiber pairs in RVE-2.</figcaption></figure>
### 3. Finite Element Simulations

A two-dimensional finite element modeling framework of multiple fiber-matrix debonding initiation and propagation as well as matrix cracking was set up in ABAQUS to replicate the tested specimen geometries. The actual fiber locations and radii measured experimentally were considered. Linear 4-node elements were used to discretize the multi-fiber specimens under the plane strain assumption. Validity and accuracy of the 2D plane strain assumption in lieu of a 3D model have been comprehensively studied in a recent article by the authors [34]. Instead of modeling the whole tensile specimens, a rectangular sub-model corresponding to the specimen width encompassing all the fibers was selected. Boundary conditions were applied on the top and bottom edges of the model, whereas the lateral edges were set free to replicate free edges in the experiments. Vertical and horizontal displacement boundary conditions were applied on all the mesh nodes lying on the sub-model bottom and top edges. These boundary conditions were directly extracted from the displacement fields measured by DIC in order to account for possible asymmetry due, for instance, to slight specimen misalignment in the tensile grips. The epoxy matrix was modeled as a linear elastic material with $E = 3.04$ GPa Young's modulus and $v=0.41$ Poisson's ratio. Fiber-matrix debonding was modeled using cohesive zone elements with a bilinear traction-separation profile [35-37]. The tensile and shear strengths, as well as the critical energy release rate (ERR), were identified based on an inverse identification approach on single-fiber specimens reported in Girard et al. [23]. These parameters were identified to reproduce the measured debonding angle variation as a function of the applied loading in single-fiber specimens. The identified interface strengths were $\sigma_c = 3.5$ MPa and $\tau_c = 4.5$ MPa, and the interface critical ERR was defined by $G_{Ic} = 11$ J/m$^2$ and $G_{IIc}/G_{Ic} = 3$. To account for the variability in the measured interface properties [27], we also studied a configuration where the critical ERR at each interface
followed a normal distribution with a $G_{ic}$= 11 J/m² mean and 1.5 J/m² standard deviation. The
application of such variable fracture properties was justified by the results of a recent study
wherein the inverse identification of fracture properties on three separate specimens made in
nominally similar conditions yielded slightly different $G_c$ values. This level of uncertainty
indicated that the interface fracture properties were not a single set of values, but rather a scattered
distribution [27].

The mesh size was chosen so that the process zone before debonding initiation and during
propagation contained at least 10 integration points. In addition, a phase-field fracture model was
used to describe matrix cracking based on the formulation proposed by Doitrand *et al.* [38] with
$G_c$= 400 J/m² for matrix critical ERR and $\sigma_f$= 60 MPa tensile strength, assuming a Rankine
strength surface. An AT1 length-free formulation was adopted in which the regularization length
was locally adjusted so as to follow the input Rankine strength surface based on the previously
determined correlation between the phase-field internal length and the material tensile strength
[39, 40].

The model contained around 50,000 degrees of freedom. Visual examples of the finite element
meshes are shown in **Figure 4**.
<figure><img src="image_5.png" /><figcaption>**Figure 4-** Finite element model replicas for (a) RVE-1 and (b) RVE-2 specimens. For RVE-1, displacement boundary conditions were applied with values directly obtained from experimental DIC measurements. Idealized boundary conditions (perfectly uniaxial loading) were used for RVE-2. (c) Meshing of the fibers and the matrix.</figcaption></figure>

# 4. Results and Discussion

## 4.1. Tensile Response of the Neat Epoxy

**Figure 5a** shows the tensile stress-strain response of the fully-cured neat epoxy. The neat epoxy used in this work shows a true failure strain slightly higher than 0.03. The tensile mechanical properties of the epoxy indicated a Young's modulus of 3.04±0.11 GPa and a tensile strength of 59.5±2.0 MPa. Poisson's ratio of the epoxy was determined to be 0.41±0.03, characterized as the ratio between DIC-measured transverse and axial strain components [21]. Variation of stiffness (=dσ/dε) with strain was calculated by applying a forward difference operator on the global stress-strain data. The stiffness curve in **Figure 5b** shows a linearly decreasing trend. This curve will be used in later sections to compare the stiffness variation in the epoxy with and without fibers and
to highlight the role of fiber-matrix debonding as a progressive damage mechanism. Note that the
oscillations shown in the stiffness-strain data at small strains are the direct artifact of the numerical
differentiation, which significantly magnifies the noise level in the raw (i.e., stress-strain) data.

<figure><img src="image_6.png" /><figcaption>**Figure 5-** (a) True stress-strain and (b) stiffness-strain curves for the epoxy resin. Scatter bars indicate experimental variability (standard deviation) among three independent measurements. Data in this figure represent the global stress-strain response of the epoxy.</figcaption></figure>

4.2. Global vs. Local Strain Fields in Macro Fiber Specimen RVE-1

The global (macroscale) tensile stress-strain response of the macro fiber specimens was studied.

**Figure 6a** shows the stress-strain response for RVE-1 with its corresponding stiffness-strain curve in **Figure 6b**. The stress values in **Figure 6a** were calculated as the ratio of the measured force and the original cross-sectional area of the tensile test specimen. Global strain values were measured by DIC with a vertical optical extensometer overlaid on the DIC area of interest, directly measuring the elongation of this area, as shown schematically by the inset in **Figure 6a**.
<figure><img src="image_7.png" /><figcaption>**Figure 6-** (a) Macroscale (global) stress-strain curve for RVE-1. The figure inset shows the dimensions of the optical extensometer used for global strain measurement. (b) Stiffness-strain curve for RVE-1 and the neat epoxy. (c) Evolution of full-field in-plane axial strain component, εyy, at different stages of deformation. Circled numbers in (c) correspond to those in (b).</figcaption></figure>

The stress-strain response in **Figure 6a** exhibits a concave nonlinear curve with a failure strain of approximately 0.01, i.e., a three-fold decrease compared with the failure strain of the neat epoxy. The corresponding stiffness-strain curve shows initial high stiffness values that rapidly decrease with strain. The rate of decrease of stiffness with strain reduces at global strains greater than 5×10⁻⁴. The stiffness-strain curve of the neat epoxy (previously discussed in **Sec. 4.1**) is also shown in **Figure 6b** to facilitate a direct comparison of stiffness evolution in macro fiber composite and neat epoxy specimens. Despite the initial higher stiffness of the composite specimen which
originates from the inclusion of rigid fibers, the stiffness reduction in the composite specimen is more substantial than that of the neat epoxy. The significant stiffness reduction in the composite specimen is correlated with the development of fiber-matrix debonding in the specimen. The contour maps in **Figure 6c** represent the evolution of the axial strain component, $\epsilon_{yy}$, in RVE-1 at increasing global strains. The global strains associated with these contour maps are marked by circled numbers in **Figure 6b**. The development of fiber-matrix interface debonding reveals itself in the form of narrow high-strain bands located at the top or bottom poles of the fibers. The location of these high-strain bands depends on the local strain and stress fields relative to the nearest neighboring fiber. Development of fiber-matrix debonding progresses steadily until the complete failure of the specimen, depicted by the last contour map in **Figure 6c**. The strain release following the final failure leads to a reduction in the intensity and heterogeneity of the axial strain field.

### 4.3. Fiber-Matrix Interface Debonding and Progressive Damage Evolution in RVE-1

The continuous development of fiber-matrix interface debonding at different locations is regarded as a progressive damage mechanism that reduces the stiffness of the macro fiber composite specimens. Given the resolution of the experimental measurements conducted here, it is impractical to accurately visualize the progression of fiber-matrix debonding at different stages of the loading. Furthermore, the interface crack openings due to debonding are asymmetric and often too small to be detected visually. However, our previous studies on single macro fiber specimens made with the same fiber and matrix materials indicated that upon debonding initiation, areas in the vicinity of the fiber-matrix interface will experience a local opening strain of ~0.002 [22, 26]. This critical strain was used herein as a threshold to detect and manually count the number of debonded fibers at different global strains.
**Figure 7a** shows examples of the debonded fibers indicated by binary color maps created based on the abovementioned 0.002 strain criterion. The debonded interfaces are clearly shown with red bands developed at the top and/or bottom poles of the fibers at increasing global strains. The visualization approach utilized to generate **Figure 7a** has practical significance as it allows the number of debonded interfaces to be counted and correlated with the global deformation and stiffness response in the specimen. **Figure 7b** shows the evolution of the number of debonded interfaces with global strain for the specimen RVE-1. Note that the scatter bars denote a ±1 experimental variability due to human error as a possible source of uncertainty. The number of debonded interfaces increases with the global strain in a piecewise trilinear trend, with zero debonded interfaces detected at global strains less than ~0.002. The majority of interfaces debond at global strains ranging from 0.002 to 0.006. The remaining intact interfaces will debond at global strains greater than 0.006. Finally, as the global strain approaches 0.01, the specimen fails completely with the crack path shown in **Figure 7a**. The final crack observed for RVE-1 indicates a path that initiates between the left edge of the specimen and fiber No. 1 (Numbering scheme shown earlier in **Figure 1d**). The crack bridges between fibers No. 1 and 2, and then between 2 and 3. A final inter-fiber crack bridging takes place between fibers 3 and 4. Having reached almost half of the width of the specimen, the crack will then propagate toward the right edge of the specimen at a straight line with no interactions with the remainder of the fibers. It is noteworthy to mention that the use of a brittle epoxy and the measurement resolution in this work did not allow for the detection of the matrix cracking as an intermediate stage between fiber-matrix interface debonding and transverse crack formation. Evidence for matrix microcracks originating from a debonded interface in similar macro fiber systems was visualized and discussed in our previous works [21].
<figure><img src="image_8.png" /><figcaption>**Figure 7-** (a) binary contour maps showing the evolution of damage at different fiber locations in RVE-1. The debonded fibers are indicated by yellow star markers. The final crack path is highlighted by a thick solid white line. (b) Variation of the number of debonded fibers, $n_d$, with global strain determined visually.</figcaption></figure>

Comparing the stiffness curve discussed earlier (Figure 6b) and the debonded interface numbers presented in Figure 7b, there appears to be a slight discrepancy between the two data sets. The stiffness values show a significant drop in the range of 0 to 0.0005 global strain which is likely due to debonding initiation. However, the first interface debonding detected visually and reported in Figure 7b corresponds to a ~0.002 global strain. This discrepancy leads us to assume that the visual detection of the debonded interfaces would not necessarily be the most accurate method of assessing the damage progression in the macro fiber specimens, at least when it is conducted on only one surface of the specimen. Note that, as was also shown in our previous work [26], debonding can initiate from one or both exposed surfaces (front and rear), rapidly tunneling along the interface to the other side of the specimen. Here, in-plane debonding may have occurred on all or a few fibers from the other side of the specimen where visuals were unavailable. Therefore, the
debonding initiation on the rear side of the specimen is completely neglected, possibly contributing to the discrepancies discussed earlier. One potential solution to mitigate such a discrepancy is using a synchronized three-camera system to view damage detection on the front, rear, and side of the specimen simultaneously. While this idea remains an interesting topic of future research, for the purposes of correlating global (macro) and local (micro/meso) scale mechanics and failure, we will continue to use the data extracted via the single front-view camera hereafter.

The number of debonded interfaces in **Figure 7b** can be further analyzed to establish correlations between the stiffness reduction and the progression of damage in the macro fiber specimens.

**Figure 8** shows the variation in stiffness with respect to the volume fraction of debonded interfaces (i.e., the ratio between the failed interfaces, $n_d$, to the total number of fibers, $n_f=39$). The stiffness values are also normalized with the initial Young's modulus of the neat epoxy to enable a direct comparison of our data with those reported in previous studies. Finally, the data presented in **Figure 8** are truncated to include only the values after 0.002 global strain to ensure that the possible sources of discrepancy discussed above are excluded from the analyses. This figure compares the stiffness-debonding data obtained in this work with the numerical results reported by Teng [41].

The data reported in the work of Teng [41] represent composites with a fiber-to-matrix moduli ratio, $E_f/E_m=20$, which is comparable to the results of the present study (i.e., $72/3.04=23.7$).

Furthermore, as depicted in **Figure 8**, the simulation results reported by Teng [41] include fiber volume fractions of 30%, 50%, and 70%. Regardless of the fiber volume fraction, simulation results predict a monotonic reduction in stiffness with the volume fraction of debonded fibers, with higher fiber volume fractions showing a steeper decrease at low debonded fiber fractions. A decreasing trend is observed in our experimental measurements, noting that the fiber volume fraction in RVE-1 is approximately 20%, leading to a more modest decrease as expected by the
simulation curve trends. In addition, the simulation results of Teng [41] indicate that at a debonding
volume fraction greater than 30%, the stiffness reduction rate is moderated, minimizing the effect
of fiber volume fraction.

<figure><img src="image_9.png" /><figcaption>**Figure 8-** Normalized stiffness versus debonded fiber volume fraction. Solid lines represent the simulation data obtained by Teng (2007). Hollow circles denote the experimental data obtained in this work for the RVE-1 specimen with a 20.5% fiber volume fraction.</figcaption></figure>

It is necessary to reiterate here that the fiber volume fractions considered in this work are
excessively lower than those in real structural composites. For a given debonded fiber volume
fraction, more fiber-matrix debonding will occur per unit volume in a composite with higher
volume fractions. Consequently, while a higher volume fraction composite is expected to show
higher initial stiffness, the rate of stiffness decrease is expected to be higher for higher volume
fraction composites [41]. Despite the low fiber volume fractions considered in the present study,
we believe the experimental methodologies developed here can be transferred to real composites
with much higher fiber volume fractions. In addition, the dilute concentration of fibers in the
samples tested here provides wider matrix regions around the fibers, allowing us to better visualize
the possible interactions between a debonded fiber and its surrounding matrix.

4.4. *Mechanical and Failure Behaviors of specimen RVE-2*

While the data obtained from RVE-1 was used primarily to uncover some of the fundamental
mechanical and failure characteristics of the macro fiber specimens as well as to calibrate the finite
element models, RVE-2 datasets were also obtained and are briefly discussed here for
completeness. **Figure 9** depicts the key results obtained for this specimen. The global stress-strain
and stiffness variations of RVE-2 are shown in **Figure 9a** and **Figure 9b**, respectively. The general
shapes and trends observed for RVE-2 resemble those discussed in the case of RVE-1.
Nevertheless, the tensile strength and initial stiffness of RVE-2 appear to be slightly lower than
those of RVE-1 possibly due to its lower fiber volume fraction. The global stress-strain response
of RVE-2 is also exemplified by a nonlinear curve with a decreasing stiffness caused by the
progressive development of damage in the form of fiber-matrix interface debonding. The evolution
of such fiber-matrix debonding at different global strains is evidenced in **Figure 9c** by the
formation of high-strain bands surrounding fibers at different locations and at increasing global
strains. Compared with the previous case, RVE-2 shows a more intense development of local strain
bands at the right edge of the specimen, which may be either due to slight misalignment of the
specimens in the tensile grips or the higher fiber packing in these regions. Regardless, fiber-matrix
debonding and the resultant microcrack formation and coalescence appear to be more severe in the
vicinity of fiber numbers 15, 39, and 21 (see **Figure 2d** for the numbering scheme). As depicted
in **Figure 9d**, the final transverse crack for RVE-2 follows the shortest path connecting fibers 21-39-17-14-...-38-13.

<figure><img src="image_10.png" /><figcaption>**Figure 9-** (a) global stress-strain and (b) stiffness-strain curves for RVE-2. The stiffness curve of the neat epoxy is included in (b) for reference. (c) Evolution of full-field in-plane axial strain component, $\varepsilon_{yy}$, at different global strains with circled numbers in (c) corresponding to those in (b). (d) Post-failure image of the RVE-2 showing the transverse crack path.</figcaption></figure>
4.5. Fiber-Matrix Debonding and Crack Path– Comparing Experiments and Simulations

The experimental observations were compared with the results of the numerical simulations.

**Figure 10** compares the experimental stress-strain curve with its numerically predicted counterpart for RVE-1. The model prediction correctly replicates the initial part of the stress-strain curve up to around 0.4% strain. Until this point, the nonlinearity is mainly caused by the fiber-matrix interface debonding. Therefore, it appears that the influence of the fiber-matrix debonding on the stress-strain curve is well captured by the model. However, at larger global strains, the nonlinearity becomes more pronounced in the experimental stress-strain curve. This observation can be explained by the matrix nonlinearities that are activated during actual testing but not considered in the numerical model, recalling our assumption of a linear elastic mechanical behavior for the epoxy.

<figure><img src="image_11.png" /><figcaption>**Figure 10- Comparing stress-strain curves obtained experimentally and numerically for RVE-1.**</figcaption></figure>
Next, the variation of the number of debonded interfaces as a function of the imposed global strain was determined. As shown in **Figure 11**, the assumption of constant values of interface strengths and critical interface ERR do not accurately replicate the debonding process assessed experimentally, because all fiber-matrix interfaces will suddenly debond within a 0.2% strain range following the first debonding initiation. In such conditions, the model predictions indicate a steep increase in the number of debonded interfaces, saturating at a global strain of ~0.5%.

<figure><img src="image_12.png" /><figcaption>**Figure 11**- Number of debonded fibers vs. global strain, comparing experimental and simulation data. Simulation data include those obtained with both constant and variable (normal distribution) fiber-matrix interface $G_c$ values.</figcaption></figure>

A more realistic description of the overall debonding process is achieved by considering the scattering and uncertainty quantification in $G_c$. We assumed a constant value for the interface critical energy release rate along the circumference of a single fiber-matrix interface. The scatter in properties is implemented by assigning different critical energy release rates to different fiber-matrix interfaces. Note that the range of variation of $G_c$ provided by the chosen normal distribution
(see **Sec. 3**) was consistent with the range of critical ERRs identified on several single-fiber
specimens [27]. The scatter in the fracture properties could be explained by the fiber surface
conditions and the variations in adhesion characteristics along the fiber-matrix interface, also
leading to slight variations from one fiber to the other. These unavoidable sources of discrepancy
highlight the importance of variability in fracture properties when considering crack initiation at
the fiber scale in a multi-fiber composite. These findings are consistent with previous works related
to crack initiation [42, 43]. For instance, Vereecke et al. [42] found that considering failure
property variability was essential to replicate the transverse crack density variation as a function
of the applied loading in composite laminates at the mesoscopic scale. In Doitrand et al. [43], the
crack density in thin films was found to be too abrupt compared to experimental observations when
no variability was considered in the model. Consistent with the previous reports, a direct
comparison between the present experimental and simulation data indicates that realistic and
reliable prediction of fiber-matrix debonding requires the inclusion of variations associated with
the uncertainty quantification of interface fracture properties.

As a complement to the overall debonding kinetics, a more in-depth analysis of fiber-matrix
debonding was performed. **Figure 12a** shows the experimental and model-predicted variations of
the interface opening displacement, $\delta_n$, at all fiber-matrix interfaces over a 1 mm gauge length as
a function of the applied global strain. Details of the procedure followed to extract the experimental
debonded interface opening are discussed in our previous publications [22]. Overall, similar trends
and consistent values for the interface openings are obtained, indicating that the general debonding
trends are correctly replicated by the numerical model here. This analysis can be further refined
by comparing the opening displacement variation for any given fiber in the model, as shown by
several examples in **Figure 12b-d**. While the opening displacements are accurately captured by
model predictions (see e.g., fibers No. 1 and No. 36 in RVE-1), there are also some discrepancies
in the opening displacement variation as a function of global strain for other fibers such as fiber
No. 13. In the latter fiber, the numerical prediction overestimates the experimentally measured
values. This difference is likely due to the simplifying assumptions in the model including, for
instance, the implementation of a normal distribution for critical ERR at the interface.

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

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

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

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

**Figure 12-** Comparing experimental and simulation results of interface opening displacement, $\delta_n$, as a function of global strain for (a) all 39 fibers in RVE-1, and individual fibers (b) No.1, (c) No. 13, and (d) No. 36. Experimental values were measured with a 1 mm optical gauge placed across the fiber-matrix interface.
An analysis of the final crack paths was performed for RVE-1 comparing DIC-measured boundary conditions with idealized (perfectly uniaxial) boundary conditions. **Figure 13** shows the transverse crack paths comparing these two conditions. To better visualize the differences between the two conditions, the displacement fields before matrix cracking are also illustrated. As depicted by the DIC-measured conditions, the experiments do not correspond to a perfectly uniaxial tensile loading parallel to the specimen axis. Instead, the displacement fields are slightly misoriented, possibly due to the specimen positioning in the grips. It is clearly evident that even the slightest variations in the boundary conditions can significantly affect the final matrix crack path. Furthermore, the application of realistic DIC-measured boundary conditions leads to a crack path that better resembles those observed experimentally in RVE-1, signifying the role of boundary conditions on the transverse crack path determination. Such observations also motivate the need to carefully monitor the experimental tests with DIC in order to provide the most realistic boundary conditions possible as an input for numerical simulations.

Finally, **Figure 14** shows the model prediction of the crack path in RVE-2 based on idealized (perfectly uniaxial) boundary conditions. Compared with the experimental results shown in **Figure 9d**, the model prediction of the crack path in RVE-2 better resembles the experimental observations, especially on the left-hand side of the specimen. Consistent observations between experiments and modeling results seen on the left-hand side of both of the specimens could be due to the presence of non-uniform grip pressures applied on the tested specimens. Furthermore, while the application of realistic boundary conditions and the variability of interface fracture properties lead to more representative and realistic modeling output, there appear to be other considerations required to further converge model predictions with experimental observations. The inclusion of matrix materials nonlinearity and three-dimensional stress-strain states [34] are identified as two
major contributing factors. Last but not least, the evolution of interface debonding into transverse
cracks is accompanied by complex damage developments along the fiber axes (e.g., crack
tunneling, as revealed and discussed in previous studies [15, 26, 27]). The inclusion of the above
parameters could lead to significant improvements in the accuracy and reliability of modeling
frameworks for the study of transverse cracks in fiber composites. Nevertheless, as shown by the
results of the present study, mechanism identification and quantitative data collected accurately
from single fiber tests can indeed be a reliable and useful first step in understanding damage
characteristics in multi-fiber composites, at least for composites with low fiber volume contents.
The representativeness of such an approach is further justified by the acceptable predictions of
crack initiation and final crack path despite the use of many simplifying assumptions.
<figure><img src="image_17.png" /><figcaption>**Figure 13-** Displacement fields (left) and matrix crack paths (right) for RVE-1 with (a) DIC-measured boundary conditions and (b) idealized displacement boundary conditions applied as a constant vertical displacement at the top edge of the domain.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 14-** Distribution of phase field damage variable indicating transverse crack path in RVE-2. Experimental observation (dotted line) and model predictions are shown. Modeling results are based on idealized displacement boundary conditions.</figcaption></figure>

**5. Conclusions**

A hybrid experimental-modeling approach was implemented to explore the development of fiber-matrix interface debonding and its evolution to transverse cracking in model macro fiber composite specimens. Experimental full-field DIC measurements revealed different stages associated with the initiation and progression of fiber-matrix interface debonding, leading to transverse cracks. Experimental quantification of fiber-matrix interface debonding suggested a nonlinear damage progression trend, leading to a significant decrease in stiffness as a function of damage evolution. The stiffness reduction was found to be most substantial at the early stages of interface debonding.

The experimental results were replicated by a modeling framework that combined cohesive zone and phase field models to simulate fiber-matrix interface debonding and matrix cracking. At the onset, the model was set up based on constant (unique) fiber-matrix interface fracture properties identified previously from single macro fiber test results. It was shown that if a unique set of fracture properties is considered for all the interfaces, the predicted debonding kinetics is too steep
compared to those measured experimentally. The predicted debonding kinetics converged to the experimentally measured values when a normal distribution (i.e., variability) of the interface critical ERR was included in the model. This observation highlighted the importance of uncertainty and scattering in the interface fracture properties, which was also consistent with previous observations where slightly different ERR values were identified in the testing of single fiber specimens.

Despite the simplifying assumptions and the implementation of fracture property variability only through the critical ERR, the debonding kinetics were overall well captured by the model. In particular, simulations captured the variation of the number of debonded fibers and debonding opening variation with increasing loading, and also the nonlinearity in the stress-strain curve due to debonding for small strain levels. The modeling results also indicated that it is crucial to consider realistic boundary conditions to accurately capture the final crack path. Idealized boundary conditions provided acceptable crack path predictions for the sample with lower fiber volume fraction with longer fiber-to-fiber distances. Further improvements in the modeling must include the nonlinearity in the stress-strain response of the matrix material, implement 3D geometries, and consider variability in fracture parameters beyond the critical ERR.

## Acknowledgments

Z.H. gratefully acknowledges the financial support provided by the New Jersey Economic Development Authority (NJEDA) through the Wind Institute Fellowship. B.K. acknowledges the financial support provided by the Advanced Materials & Manufacturing Institute at Rowan University.
## Data availability

Data will be made available on request.

## References

[1] Zhuang L, Talreja R, Varna J. Transverse crack formation in unidirectional composites by linking of fibre/matrix debond cracks. *Compos Part A: Appl Sci Manufact* 2018; 107: 294-303.

[2] Talreja R, Yalvac S, Yats LD, Wetters DG. Transverse cracking and stiffness reduction in cross ply laminates of different matrix toughness. *J Compos Mater* 1992; 26(11): 1644-1663.

[3] Guillén-Hernández T, García IG, Reinoso J, Paggi M. A micromechanical analysis of inter-fiber failure in long reinforced composites based on the phase field approach of fracture combined with the cohesive zone model. *Int J Fract* 2019; 220: 181-203.

[4] Guillén-Hernández T, Quintana-Corominas A, García IG, Reinoso J, Paggi M, Turón A. In-situ strength effects in long fibre reinforced composites: A micro-mechanical analysis using the phase field approach of fracture. *Theor Appl Fract Mech* 2020; 108: 102621.

[5] Carollo V, Reinoso J, Paggi M. Modeling complex crack paths in ceramic laminates: A novel variational framework combining the phase field method of fracture and the cohesive zone model. *J Eur Ceram Soc* 2018; 38(8): 2994-3003.

[6] Herráez M, Mora D, Naya F, Lopes CS, González C, Llorca J. Transverse cracking of cross-ply laminates: A computational micromechanics perspective. *Compos Sci Technol* 2015; 110: 196-204.

[7] Kohler S, Cugnoni J, Amacher R, Botsis J. Transverse cracking in the bulk and at the free edge of thin-ply composites: Experiments and multiscale modelling. *Compos Part A* 2019; 124: 105468.

[8] Arteiro A, Catalanotti G, Melro AR, Linde P, Camanho PP. Micro-mechanical analysis of the effect of ply thickness on the transverse compressive strength of polymer composites. *Compos Part A* 2015; 79: 127-137.

[9] Sutton MA, Li N, Joy DC, Reynolds AP, Li X. Scanning Electron Microscopy for Quantitative Small and Large Deformation Measurements Part I: SEM Imaging at Magnifications from 200 to 10,000. *Exp Mech* 2007; 47: 775-787.

[10] Tracy J, Daly S, Sevener K. Multiscale damage characterization in continuous fiber ceramic matrix composites using digital image correlation. *Exp Mech* 2015; 50: 5286-5299.

[11] Mehdikhani M, Aravand M, Sabuncuoglu B, Callens MG, Lomov SV, Gorbatikh L. Full-field strain measurements at the micro-scale in fiber-reinforced composites using digital image correlation. *Compos Struct* 2016; 140: 192-201.
[12] Montgomery CB, Koohbor B, Sottos N. A Robust Patterning Technique for Electron Microscopy-Based Digital Image Correlation at Sub-Micron Resolutions. *Exp Mech* 2019; 59: 1063-1073.

[13] Koohbor B, Montgomery CB, Sottos N. Identification of RVE length scale in fiber composites via combined optical and SEM digital image correlation. *Compos Sci Technol* 2022; 227: 109613.

[14] Kammers AD, Daly S. Digital Image Correlation under Scanning Electron Microscopy: Methodology and Validation. *Exp Mech* 2013; 53: 1743-1761.

[15] Martyniuk K, Sorensen BF, Modregger P, Lauridsen EM. 3D in situ observations of glass fibre/matrix interfacial debonding. *Compos Part A: Appl Sci Manufact* 2013; 55: 63-73.

[16] Tabiai I., Delorme R, Therriault D, Levesque M. In-situ Full Field Measurements During Inter-Facial Debonding in Single Fiber Composite Under Transverse Load. *Exp Mech* 2018; 58: 1451-1467.

[17] Tabiai I, Texier D, Bocher P, Therriault D, Levesque M. In-situ Full Field Out of Plane Displacement and Strain Measurements at the Micro-Scale in Single Reinforcement Composites under Transverse Load. *Exp Mech* 2020; 60: 359-377.

[18] Tabiai I, Tkachev G, Diehl P, Frey S, Ertl T, Therriault D, Levesque M. Hybrid image processing approach for autonomous crack area detection and tracking using local digital image correlation results applied to single-fiber interfacial debonding. *Eng Fract Mech* 2019; 216: 106485.

[19] Nikforooz M, Verschatse O, Daelemans L, De Clerck K, Van Paepegem W. Experimental method for in-situ real-time measurement of mixed mode fibre/matrix interface debonding and comparison of fibre sizings in single short fibre composites. *Compos Sci Technol* 2024; 249: 110488.

[20] Verschatse O, Nikforooz M, De Clerck K, Van Paepegem W, Daelemans L. Analysis of fiber-matrix interface debonding under normal loading using in-situ microscopy and cruciform-shaped single glass fiber/composite specimens. *Compos Sci Technol* 2025; 263: 111089.

[21] Livingston R, Koohbor B. Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements. *Compos Part C* 2022; 7: 100229.

[22] Girard H, Koohbor B, Doitrand A, Livingston R. Experimental characterization of in-plane debonding at fiber-matrix interface using single glass macro fiber samples. *Compos Part A* 2023; 171: 107573.

[23] Girard H, Doitrand A, Koohbor B, Rinaldi RG, Godin N, Long D, Bikard J, Trouillet-Fonti L. Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties. *Eng Fract Mech* 2023; 286: 109254.

[24] Girard H, Doitrand A, Koohbor A, Rinaldi RG, Godin N, Long D, Bikard J. Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination. *J Mech Phys Solids* 2024; 183: 105498.
[25] Takahashi K, Shoya R, Matsuo T, Sato W, Nakamura T, Takeuchi A, Uesugi M, Uesugi K. X-ray nanoimaging of a transversely embedded carbon fiber in epoxy matrix under static and cyclic loads. *Sci Report* 2022; 12: 8843.

[26] Uddin KZ, Girard H, Mennie NB, Doitrand A, Koohbor B. Simultaneous Measurement of Fiber-Matrix Interface Debonding and Tunneling Using a Dual-Vision Experimental Setup. *Exp Mech* 2024; 64: 1497-1511.

[27] H. Girard, A. Doitrand, M. Tosti Umemura, B. Koohbor, R. G. Rinaldi, N. Godin, J. Bikard Towards the determination of optimal 3D initiation crack shapes. HAL (preprint) 2025.

[28] Chatziathanasiou T, Lee Y, Villanova J, Stamati O, AhmadvashAghbash S, Fazlali B, Breite C, Sinclair I, Mavrogordato MN, Spearing SM, Mehdikhani M, Swolfs Y. Questioning the Representativeness of Damage Mechanisms in Single-Fiber Composites via In Situ Synchrotron X-Ray Holo-Tomography. *Small* 2025; 20(52): 2406168.

[29] International Digital Image Correlation Society, Jones EMC and Iadicola MA (Eds.) 2018. A Good Practices Guide for Digital Image Correlation. DOI: 10.32720/idics/gpg.ed1

[30] Swaminathan S, Ghosh S, Pagano NJ. Statistically Equivalent Representative Volume Elements for Unidirectional Composite Microstructures: Part I - Without Damage. *J Compos Mater* 2006; 40(7): 583-604.

[31] Geoffre A, Moulin N, Bruchon J, Drapier S. Reappraisal of Upscaling Descriptors for Transient Two-Phase Flows in Fibrous Media. *Transp Porous Media* 2023; 147: 345-374.

[32] Geoffre A, Ghestin M, Moulin N, Bruchon J, Drapier S. Bounding transverse permeability of fibrous media: a statistical study from random representative volume elements with consideration of fluid slip. *Int J Multiphas Flow* 2021; 143: 103751.

[33] Trias D, Costa J, Turon A, Hurtado JE. Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers. *Acta Mater* 2006; 54: 3471-3484.

[34] Girard H, Doitrand A, Koohbor B, Rinaldi RG, Godin N, Bikard J. Comparison between 2D and 3D fiber-matrix debonding simulation for inverse identification of interface fracture properties. *J Theor Comput Appl Mech* 2025 (in press).

[35] Mi Y, Crisfield MA, Davies GAO, Hellweg HB. Progressive Delamination Using Interface Elements. *J Compos Mater* 1998; 32(14):1246-1272.

[36] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. *Int J Numer Methods Eng* 2001;50: 1701-1736.

[37] Alfano G, Sacco E. Combining interface damage and friction in a cohesive-zone model. *Int J Numer Methods Eng* 2006; 68(5), 542-582.
[38] Doitrand A, Molnár G, Estevez R, Gravouil A. Strength-based regularization length in phase field fracture. *Theor Appl Fract Mech* 2023;124, 103728.

[39] Molnár G, Doitrand A, Estevez R, Gravouil A. Toughness or strength? Correlation between length scale and tensile strength in phase-field fractures based on the coupled criterion. *Theor Appl Fract Mech* 2020;109, 102736.

[40] Molnár G, Doitrand A, Jaccon A, Prabel B, Gravouil A. Thermodynamically consistent linear-gradient damage model in Abaqus. *Eng Fract Mech* 2022; 266, 108390.

[41] Teng H. Transverse stiffness properties of unidirectional fiber composites containing debonded fibers. *Compos Part A* 2007; 38: 682-690.

[42] Vereecke J, Bois C, Wahl JC, Briand T, Ballère L, Lavelle F. Explicit modelling of meso-scale damage in laminated composites – Comparison between finite fracture mechanics and cohesive zone model. *Compos Sci Technol* 2024; 253: 110640.

[43] Doitrand A, Rubeck S, Meille S, Chevalier J, Steyer P, Gallois-Garreignot S. Influence of debonding and substrate plasticity on thin film multicracking. *Theor Appl Fract Mech* 2024; 131: 104374.