# Numerical estimation and optimization of the fracture toughness in short-fibre ceramic matrix composites

Sara Jimenez Alfaro, Hyung-Jun Chang, Dominique Leguillon, Sébastien

Denneulin

► To cite this version:

Sara Jimenez Alfaro, Hyung-Jun Chang, Dominique Leguillon, Sébastien Denneulin. Numerical estimation and optimization of the fracture toughness in short-fibre ceramic matrix composites. 2025. hal-05048399v2

HAL Id: hal-05048399

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

Preprint submitted on 1 Sep 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
# Numerical estimation and optimization of the fracture toughness in short-fibre ceramic matrix composites

Jimenez-Alfaro, Sara<fnref n="1" />, Chang, Hyung-Jun<fnref n="2" />, Leguillon, Dominique<fnref n="3" />, and Denneulin, Sebastien<fnref n="4" />

Achieving greater engine efficiency by increasing gas temperatures is a key challenge set by the European Commission for the aerospace industry by 2050. However, traditional metal alloys may fail under such extreme conditions. Ceramic matrix composites (CMCs) have emerged as a promising alternative, particularly short-fibre CMCs, which offer enhanced performance in components with complex geometries. This paper presents a novel computational framework to predict the fracture toughness of short-fibre CMCs based on the Coupled Criterion approach. The influence of various composite parameters has been systematically analysed and validated against existing experimental data, demonstrating the model's reliability and potential as a robust design tool for next-generation aerospace applications.

Keywords coupled criterion, ceramic matrix composites, fracture toughness, aeroengines

## 1 Introduction

In 2019, the European Commission introduced the European Green Deal, setting ambitious targets to achieve net-zero greenhouse gas emissions by 2050 (European Commission 2019). In response, the aerospace industry has focused on reducing aircraft weight and improving engine efficiency. A key strategy to achieve these goals involves increasing the operating gas temperature from approximately 1700°C to a target of more than 2000°C, whereas current metallic alloys used for the gas turbine engine cannot resist a temperature higher than 1300°C as indicated in Yasmin Legerstee 2022. This limitation can be partially addressed through thermal barrier coatings and cooling techniques, although it cannot be completely overcome. This substantial temperature rise presents both challenges and opportunities for material development and engineering innovation.

To achieve this objective, a transition from conventional metal alloys in aircraft engines is essential, as these materials would be unable to perform effectively at the proposed temperatures, even with the application of advanced thermal barrier coatings. Ceramics present a compelling alternative, offering numerous advantages such as low density and exceptional refractory properties, while maintaining mechanical strength at elevated temperatures. Moreover, ceramic matrix composites (CMCs) offer a distinct advantage by enhancing fracture toughness while retaining high thermal stability, low density, and chemical inertness, as explained in Donald et al. 1976. Consequently, they were identified two decades ago as a key material system for improving the thrust-to-weight ratio of high-performance aircraft engines, as indicated in Ohnabe et al. 1999. Their applications extend beyond aviation, encompassing diverse industrial sectors such as telecommunications, modern electronics, photonics, multifunctional sensors, thermal engines and military technologies, see Dhanasekar et al. 2022. In addition, it was established that the fracture toughness of the CMCs can be enhanced by introducing weak interfaces between the fibre and the matrix. Then Martin, Peters, et al. 1998; Martin, Leguillon, and Lacroix 2001 considered crack deflection mechanism at the interface between fibre and matrix as a principal toughening mechanism of the CMCs. This mechanism of microscale cracking in the CMCs was further investigated by many other researchers. (See Oriol Gavalda-Diaz et al. 2025 for a recent critical review) While many of these applications predominantly employ long-fibre CMCs, short-fibre

<footnote n="1">¹ Departamento de Mecánica de Medios Continuos y Teoría de Estructuras, Universidad de Sevilla, Camino de los Descubrimientos, 41092, Sevilla, España</footnote>

<footnote n="2">² DST Safran Tech, Rue des Jeunes Bois - Châteaufort, 78722 Magny Les Hameaux Cedex, France</footnote>

<footnote n="3">³ Institut Jean Le Rond d'Alembert, Sorbonne University, 75005, Paris, France</footnote>

<footnote n="4">⁴ Safran Ceramics, Rue de Touban, 33185, Le Haillan, France</footnote>
CMCs offer an optimal solution for components with complex geometries, such as high-pressure turbine blades. In 2022, a short-fibre CMC was characterised both numerically and experimentally in Yasmin Legerstee 2022 as a candidate for the next generation of turbine blades. While the study estimated nearly all the mechanical properties of this material, the fracture toughness of the composite was not numerically analysed, despite its critical importance.

On the other hand, the Coupled Criterion (CC) is a fracture mechanics approach used to predict crack nucleation and initial propagation. In recent years, it has gained recognition for its application to ceramic materials. This method has been applied to describe the concept of tensile strength in polycrystalline ceramics, see Martin, Leguillon, Sevecek, et al. 2018, introducing a novel definition of “inherent” tensile strength applied in numerical models, as described in Leguillon, Martin, Sevecek, et al. 2018. Moreover, Jiménez-Alfaro and Leguillon 2021 explained how the CC approach has showcased its versatility in capturing fracture behavior across multiple scales, with noteworthy applications in nacre-like alumina materials, see Doitrand et al. 2020; Duminy et al. 2023. The CC has also been effectively applied in evaluating the influence of notch radius on fracture toughness in Picard et al. 2006. Significant progress has been achieved in applying the CC approach to layered ceramics, see Ševeček et al. 2016, addressing key phenomena such as internal compressive stresses in Leguillon, Martin, Ševeček, et al. 2015, edge and tunneling crack formation in Papšík, Ševeček, Hofer, et al. 2023, and thermal shock damage in Papšík, Ševeček, Schlacher, et al. 2024, among others. Recently, Jiménez-Alfaro and Leguillon 2022 developed a novel computational framework, integrating the Coupled Criterion with the Matched Asymptotics (MA) approach to predict the apparent fracture toughness of glass matrices reinforced with alumina platelets. This innovative model was extended in Jiménez-Alfaro and Leguillon 2023 to account for the effects of residual stresses and composite volume fraction, offering a more comprehensive understanding of fracture behavior in these materials.

Other experimental and theoretical studies have been applied for explaining the crack path in ceramic matrix composites. Martin, Peters, et al. 1998 applied an energy criterion for predicting a crack deflection at the fibre-matrix interface. Barsoum et al. 1992 investigated the relationship between microstructural factors and matrix cracking initiation stress in fibre-reinforced ceramic matrix composites. However, the exact process by which these evolve into a macroscopic, fibre-bridged crack remained unresolved.

In this context, this paper aims to apply the Coupled Criterion to estimate the apparent fracture toughness of short-fibre CMCs and to present the predominant toughening mechanisms as a function of the geometrical parameters in the composite. A key aspect of the computational framework developed in this study is its ability to account for the influence of fibre interactions on crack propagation. The paper is organized as follows: Section 2 provides a comprehensive description of the Coupled Criterion. Section 3 details the material properties and the case study applied to illustrate the computational methodology. The main results, including a comparison of numerical predictions with existing experimental data, are presented and discussed in Section 4. Finally, Section 5 summarizes the principal conclusions drawn from this study.

## 2 The Coupled Criterion

The Coupled Criterion (CC) model enables the prediction of crack nucleation in brittle materials. This criterion is supported by prior experiments on transverse cracking in cross-ply laminates, see Weißgräber et al. 2016. According to the CC model first introduced in Leguillon 2002, two necessary conditions, when met together, are sufficient to predict crack nucleation, as fracture occurs when both are simultaneously satisfied: (i) an energy condition, based on an energy balance, and (ii) a stress condition. In Jiménez-Alfaro and Leguillon 2023 the two conditions of the CC are summarized considering a structure with N different materials for a 2D case. First, the stress condition is written as a function of the tensile stress $\sigma(s)$ prior to any crack extension, i.e., the stress component in the normal direction $n$ highlighted in Fig. 1

$$ \bar{\sigma} = \frac{\sigma(s)}{\sigma_c(s)} \ge 1, \text{ for } 0 \le s \le \delta l, \qquad (1) $$
where s denotes the coordinate along the expected crack path, see Figure 1, and $\delta l$ is the newly created crack length. The function $\sigma_c(s)$ describes the tensile strength along the expected crack path $0 \le s \le \delta l$,

$$ \sigma_c(s) = \sum_{i=1}^{N} \sigma_c^{(i)} [H(s - \delta l^{(i-1)}) - H(s - \delta l^{(i)})], \qquad (2) $$

being $H$ the Heaviside function and $\delta l^{(0)} = 0$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Expected crack path in a specimen with a heterogeneous structure made of two materials A and B.</figcaption></figure>

The energy condition for a newly created crack length $\delta l$ is written in a heterogeneous structure as

$$ G_{\text{inc}} = -\frac{\Delta\Pi_{\text{p}}(\delta l)}{\delta l}, \qquad (3) $$

and therefore

$$ \frac{G_{\text{inc}}}{\bar{G}_{c}(\delta l)} \ge 1, \qquad (4) $$

where $\Delta\Pi_{\text{p}}(\delta l)$ is the difference in potential energy between the initial state prior to any crack extension and the final state when a crack of length $\delta l$ is generated $\Delta\Pi_{\text{p}}(\delta l) = \Pi_{\text{p}}(\delta l) - \Pi_{\text{p}}(0)$. The function $\bar{G}_{c}(\delta l)$ describes the critical energy release rate along the expected crack path,

$$ \bar{G}_{c}(\delta l) = \frac{\sum_{i=1}^{N} G_{c}^{(i)} \delta l^{(i)}}{\delta l}. \qquad (5) $$

being $\delta l^{(i)}$ the crack length associated with material (i). Hence, the minimum admissible load for which both Equations (1) and (4) are satisfied determines the crack onset. This, in turn, defines the nucleation length, $\delta l_c$, which represents the initial crack jump upon nucleation, within the framework of Finite Fracture Mechanics.

## 3 A computational framework to estimate the fracture toughness in CMCs

This section presents the computational tool employed to numerically assess the fracture toughness of short fibre ceramic matrix composites (CMCs). Section 3.1 provides an overview of the material properties used, alongside the experimental determination of fracture toughness. Section 3.2 outlines the applied methodology, which utilises the Matched Asymptotics technique
to derive the solution in the vicinity of the crack tip under Mode I loading conditions. This solution is then linked to the composite's fracture toughness through the Coupled Criterion. The proposed tool significantly reduces the computational cost associated with such problems, where fine mesh resolution is necessary to accurately capture the behaviour in the region surrounding the small-scale constituents of the matrix (the short fibres).

## 3.1 The material

This paper examines two ceramic matrix composites (CMCs) reinforced with short fibres, both of which have been numerically and experimentally characterised in Yasmin Legerstee 2022. The first CMC (named in the text as **Material 1**) features a matrix composed of TiSi₂ and MoSi₂. Notably, MoSi₂ undergoes a striking transformation from brittle behaviour at low temperatures to ductile behaviour at elevated temperatures (around 900°C), making it an ideal candidate for high-pressure turbine blade applications. Complementing this, TiSi₂ is strategically incorporated to counteract the oxidation of MoSi₂ in air — a well-known degradation phenomenon termed "pest". The composite is further reinforced with silicon carbide (SiC) Hi-Nicalon fibres, selected for their outstanding high-temperature performance and their ability to overcome the oxidation limitations typical of conventional SiC fibres. The second CMC (denoted in the text as **Material 2**) retains the same matrix and fibre composition but introduces a Pyrocarbon (PyC) interphase. This addition significantly enhances the composite's overall resilience and mechanical performance.

The two materials are characterised both numerically, using statistical methods, and experimentally. In both cases, the fibre volume fraction is approximately $V_f \approx 12\%$. The fibres exhibit a length range from a minimum of $l_{min} = 24.41 \ \mu\text{m}$ to a maximum of $l_{max} = 2272 \ \mu\text{m}$, while their diameters vary between $t_{min} = 8.500 \ \mu\text{m}$ and $t_{max} = 53.30 \ \mu\text{m}$. Based on statistical distributions, the average fibre length and diameter are determined to be $l_{avg} = 152.0 \ \mu\text{m}$ and $t_{avg} = 16.67 \ \mu\text{m}$, respectively.

Table 1 presents the mechanical properties of the constituents, as reported in Yasmin Legerstee 2022. The material properties of the composite are obtained applying the Voigt rule, and depend on the volume fraction. Notably, for the interphase, while the values of $\sigma_c$ and $K_c$ are experimentally determined for the PyC material, the Young's modulus provided in Yasmin Legerstee 2022 corresponds to that of BN, due to confidentiality reasons. For the sake of consistency, this same value is adopted in the present analysis. Notice that the interface between the interphase, the fibre and the matrix is considered strong.

<table><caption>**Table 1:** Mechanical properties of the constituents of the composite material under study, as reported in Yasmin Legerstee 2022.</caption><thead><tr><th>Constituents</th><th>E [GPa]</th><th>σ<sub>c</sub> [MPa]</th><th>K<sub>c</sub> [MPa m<sup>1/2</sup>]</th></tr></thead><tbody><tr><td>Fibre</td><td>420</td><td>1222</td><td>1.74</td></tr><tr><td>Interphase</td><td>20</td><td>107</td><td>0.35</td></tr><tr><td>Matrix</td><td>359</td><td>200</td><td>2.42</td></tr></tbody></table>

It is worth noting that in Yasmin Legerstee 2022, a study on the sensitivity of the elastic properties in the three spatial directions is conducted, concluding that the differences are not significant, and thus considering anisotropy is not necessary for this material.

The fracture toughness of the composite is estimated in Yasmin Legerstee 2022 using an analytical approximation based on Linear Elastic Fracture Mechanics (LEFM) and experimentally determined through Digital Image Correlation (DIC). In the first one, a shape factor $f(\frac{a}{W})$ is used,

$$K_{\text{LEFM}} = f\left(\frac{a}{W}\right) \sigma_c \sqrt{\pi a}, \qquad (6)$$

where a and W are the notch and the specimen width, respectively (see Figure 2). The tensile strength $\sigma_c$ of the composite is obtained considering the critical force of the test $P_c$:

$$\sigma_c = \frac{3P_c(L_e - L_i)}{2eW^2} \qquad (7)$$
being *e* the thickness of the specimen, and *L*<sub>e</sub>, *L*<sub>i</sub> the position of the boundary conditions (see Figure 2). The shape factor *f* (a/W) is defined as

$$f\left(\frac{a}{W}\right) = \frac{1.1215}{\sqrt{\left(1 - \frac{a}{W}\right)^3}} \left[ \frac{5}{8} - \frac{5}{12}\frac{a}{W} + \frac{1}{8}\left(\frac{a}{W}\right)^2 + 5\left(\frac{a}{W}\right)^2\left(1 - \frac{a}{W}\right)^6 + \frac{a}{W}\exp\left(-\frac{6.1342\frac{a}{W}}{1 - \frac{a}{W}}\right) \right], \quad (8)$$

The second methodology is based on the well-known Williams's expansion in Williams 1952, that defines the elastic solution $u(r, \theta) = K_I r^{0.5} u_I(\theta)$ in the neighbourhood of the Mode I crack tip. The DIC technique is applied to obtain the displacement field at the moment of failure, and therefore $K_c = K_I$. The results are presented in Table 2.

<table><caption>**Table 2:** Fracture toughness obtained in the composite, as reported in Yasmin Legerstee 2022.</caption><thead><tr><th>Fracture Toughness</th><th>K<sub>c</sub><sup>LEFM</sup> [MPa m<sup>1/2</sup>]</th><th>K<sub>c</sub><sup>DIC</sup> [MPa m<sup>1/2</sup>]</th></tr></thead><tbody><tr><td>Material 1</td><td>2.80 - 3.45</td><td>3.50 - 4.40</td></tr><tr><td>Material 2</td><td>3.87 - 4.25</td><td>5.20 - 6.32</td></tr></tbody></table>

Additionally, the fibre orientation study presented in Yasmin Legerstee 2022 reveals that most fibres are aligned within the same plane. This alignment enables the definition of an equivalent 2D problem, akin to the approaches described in Jiménez-Alfaro and Leguillon 2022; Jiménez-Alfaro and Leguillon 2023.

## 3.2 Model description to estimate the fracture toughness

The problem addressed in this paper, as illustrated in Figure 2, is inspired by the computational tool developed in Jiménez-Alfaro and Leguillon 2022; Jiménez-Alfaro and Leguillon 2023. It consists of a 3-point bending test with a pre-existing crack, $\Gamma_c$, subjected to Mode I loading conditions. At the crack tip, a fibre is positioned perpendicularly to the crack. The fracture toughness is estimated by studying the crack propagation of $\Gamma_c$, considering two possible scenarios: the crack may either penetrate the fibre or deviate along the interface if the first material described in Section 3.1 is considered, or follow the interphase path if the second material is analysed.

To determine this, the CC approach is applied, which requires the elastic solution $U(x_1, x_2)$ in the vicinity of the crack tip. An estimation of this solution is obtained using the quasi-analytical Matched Asymptotics (MA) technique, which is extensively described in Jiménez-Alfaro, García, et al. 2025. This approach is particularly useful in minimizing numerical errors that may arise when employing very small elements near perturbations. Additionally, it helps reduce the computational complexity of problems involving features such as platelet reinforcements embedded in brittle matrices Jiménez-Alfaro and Leguillon 2022; Jiménez-Alfaro and Leguillon 2023.

In the MAE, the solution $U(x_1, x_2)$ in the vicinity of the tip of $\Gamma_c$ can be approximated by an inner and an outer expansion. In the outer expansion

$$U(x_1, x_2) = U_0(x_1, x_2) + \text{small perturbation}, \quad (9)$$

where $U_0(x_1, x_2)$ is the solution of the same elastic problem as the one indicated in Figure 2, but assuming that the influence of the short fibre located at the tip of $\Gamma_c$ is negligible. Hence, the solution can be described by the William's expansion explained in Williams 1952 that in polar coordinates $(r, \theta)$ is written as

$$U_0(x_1, x_2) = U_0(0, 0) + K_I \sqrt{r} u_I(\theta) + \dots \quad (10)$$

being $K_I$ the stress intensity factor and $u_I(\theta)$ the opening shape function for a crack under Mode I, see Jiménez-Alfaro, García, et al. 2025 for more details.

In the inner expansion, the solution is expressed through a change of variable, where the primary dimensions of the domain are non-dimensionalised by one of the characteristic lengths of the problem. This ensures that the dimensions of the specimen are significantly larger than the
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Schematic representation of the problem being solved: a 3-point bending test with a pre-existing crack, $\Gamma_c$, interacting with a perpendicular short fibre.</figcaption></figure>

selected length scale. In this study, two main parameters are considered: the length ($l$) and the diameter ($t$) of the fibre. Hence, the change of variable can be $x_i = ly_i$ or $x_i = ty_i$, respectively. If the problem depicted in Figure 2 is modified by this change of variable, an unbounded domain, $\Omega_{in}$, is obtained, under the condition that the dimensions $l$ and $t$ are much smaller than the specimen's overall dimensions. In this paper, two possible configurations of the inner domain are considered: the one proposed in Jiménez-Alfaro and Leguillon 2022, shown in Figure 3, and a new configuration that incorporates the effect of the interaction between fibres, as illustrated in Figure 4.

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Scheme of the inner domain considering a single fibre. Initial undeformed shape.</figcaption></figure>

In Figure 3, the parameter $\hat{e}$ represents the thickness of the interphase, $e$, normalised by the characteristic dimension of the fibre used in the non-dimensionalisation of the inner domain, either $l$ or $t$. The three geometric parameters $l$, $e$, and $t$ are related through a shape factor, $f_e = \frac{lt}{e}$. It is important to note that if the first material is considered, there is no interphase in the composite material, and thus $e = 0$ (see Section 3.1). In the proposed model, it is assumed that around the fibre there exists a region of a certain size, free of fibres, and composed solely of the matrix. The size of this region is defined by the fibre volume fraction, $V_f$, as explained in Jiménez-Alfaro and Leguillon 2023. For instance, when the change of variable is performed using the fibre diameter, $t$, the dimensions of this region are

$$ \hat{a}_t = \frac{t}{l\sqrt{V_f}}, \quad \hat{a}_l = \frac{1}{\sqrt{V_f}} \qquad (11) $$
The model shown in Figure 4 considers the effect of fibres in proximity to the preexisting crack $\Gamma_c$, enabling the study of potential interactions between fibres by modifying the separation distance. Two new parameters are introduced: $\hat{h}_l$ and $\hat{h}_t$, which represent the horizontal and vertical distances between fibres, respectively. The geometric parameters $\hat{a}_l$ and $\hat{a}_t$ in Figure 4 are defined as

$$ \hat{a}_l = 3\hat{l} + 2\hat{h}_l + 6\hat{e}, \quad \hat{a}_t = 3\hat{t} + 2\hat{h}_t + 6\hat{e}, \qquad (12) $$

<figure><img src="image_5.png" /><figcaption>Figure 4: Scheme of the inner domain considering multiple fibres. Initial undeformed shape.</figcaption></figure>

Some restrictions are naturally deduced from the model

$$ \hat{a}_l > 0 \quad \text{and} \quad \hat{a}_t > 0, \qquad (13) $$

$$ \text{if } \hat{h}_l\hat{h}_t > 0 \quad \text{then} \quad \hat{h}_l, \hat{h}_t > 0, \qquad (14) $$

$$ 5\hat{l}\hat{t} = V_f \hat{a}_l \hat{a}_t, \qquad (15) $$

$$ \hat{a}_l - 2\hat{l} - 6\hat{e} > 0 \quad \text{and} \quad \hat{a}_t - 2\hat{t} - 6\hat{e} > 0. \qquad (16) $$

The first restriction, as shown in Equation (13), ensures that the dimensions $\hat{a}_l$ and $\hat{a}_t$ are always positive. The second restriction, expressed in Equation (14), stipulates that one fibre cannot penetrate another in the model. Equation (15) simply reiterates that the definition of the fibre volume fraction must be satisfied. Finally, the last condition, given in Equation (16), is applied to prevent fibres from overlapping in the geometry. It is important to note that these restrictions alone are not sufficient to yield a unique solution for the geometry of the problem, as $\hat{t}$ and $\hat{l}$ remain free parameters. Consequently, we select these two parameters in such a way that they minimise the surface area of the multi-fibre model, which is defined as

$$ S(\hat{h}_l, \hat{h}_t) = \hat{a}_l \hat{a}_t = (3\hat{l} + 2\hat{h}_l + 6\hat{e}). \qquad (17) $$

Although we use this hypothesis, at the end of this section we will change it to study the influence of $\hat{h}_l$ and $\hat{h}_t$.

As an example, the development of the inner expansion is written below considering the length of the fibre as the characteristic length in the adimensionalization, but a similar development can be obtained with the diameter of the fibre. The inner expansion is therefore defined as

$$ U(x_1, x_2) = U(ly_1, ly_2) = F_0(l)V_0(y_1, y_2) + F_1(l)V_1(y_1, y_2) + \dots \qquad (18) $$

The terms of Equation (18), and consequently the elastic solution in the vicinity of the pre-existing crack tip, $\Gamma_c$, are derived by relating the outer and inner expansions in a common region where
both approximations accurately represent the actual solution. In that region, the so-called matching conditions are applied,

$$F_0(l)V_0(y_1, y_2) \approx U(0, 0), \quad \text{when } \rho \to \infty, \qquad (19)$$

$$F_1(l)V_1(y_1, y_2) \approx K_I l^{\lambda} \rho^{\lambda} u(\theta), \quad \text{when } \rho \to \infty, \qquad (20)$$

where $\rho = \sqrt{y_1^2 + y_2^2}$ and $\approx$ means "behave like". Hence, it can be set

$$F_0(l) = 1, \quad V_0(y_1, y_2) \approx U(0, 0), \quad \text{when } \rho \to \infty, \qquad (21)$$

$$F_1(l) = K_I l^{\lambda}, \quad V_1(y_1, y_2) \approx \rho^{\lambda} u(\theta), \quad \text{when } \rho \to \infty, \qquad (22)$$

The expression in Equation (18) is a good approximation of the real solution $U(x_1, x_2)$ in the vicinity of the crack tip. To obtain the solution $V_1(y_1, y_2)$, the finite element method is applied in a circular domain whose dimensions are much larger (in this case, $R = 400 \ \mu m$) than those of the internal domain described in Figures 3 and 4. The element size in the vicinity of the pre-existing crack tip $\Gamma_c$ is sufficiently small to ensure numerical convergence, with a size of 0.008 $\ \mu m$. On the boundary of the domain, a prescribed displacement is applied as a boundary condition, as shown in Figure 5.

<figure><img src="image_6.png" /><figcaption>Figure 5: Scheme of the numerical problem.</figcaption></figure>

Once the elastic solution is obtained, the two CC conditions, summarised in Equations (1) and (5), are applied, yielding

$$K_I(s) \geq \frac{\sqrt{l} \sigma_c(\hat{s})}{\sigma(\hat{s})}, \qquad (23)$$

$$K_I(s) \geq \sqrt{\frac{G_c(\hat{s}l)}{-\frac{\Delta\Pi_p}{\delta l}}}. \qquad (24)$$

The apparent fracture toughness of the composite, denoted as $K_{IC}^{app}$, corresponds to the minimum value of $K_I(s)$ that satisfies both conditions of the Coupled Criterion, as summarised in Equations (23) and (24). In Equation 23 the tensile stress $\sigma(\hat{s})$ is obtained as

$$\sigma(\hat{s}) = \frac{1}{\sqrt{l}} \nabla_y V_1(y_1, y_2) \qquad (25)$$

if the inner domain is built considering the length of the fibre $l$ or

$$\sigma(\hat{s}) = \frac{1}{\sqrt{t}} \nabla_y V_1(y_1, y_2) \qquad (26)$$

if the inner domain is built considering the diameter of the fibre $t$. The term $\nabla_y$ is the gradient operator with respect to the variables $y_1, y_2$. As explained in Jiménez-Alfaro and Leguillon 2022
the increment in potential energy $\Delta\Pi_p$ is the one directly obtained from the inner problem (solution $V_1(y_1, y_2)$).

Finally, note that the problem is solved under the assumption of plane strain elasticity to estimate fracture toughness. While this simplification offers practical advantages, it may not fully capture the complexities of the real material. Given that the fibres have a round cross-section, their thickness in the perpendicular direction is considerably smaller than their length, introducing limitations to the model's representativeness. Consequently, certain phenomena, such as crack propagation around the fibres, cannot be accurately represented. Nevertheless, despite these constraints, the model provides a reliable first estimate of fracture toughness that aligns well with experimental data, as will be demonstrated in this section.

# 4 Results and discussions

This section presents the key results obtained from the application of the numerical tool described in Section 3.2 considering the material outlined in Section 3.1. The effect of various geometric parameters of the composite is analysed. Section 4.1 explores the influence of fibre size, proposing a configuration that optimises the fracture toughness of the composite. Section 4.2 investigates the impact of fibre volume fraction, as well as the changes in fracture properties that may arise from including or excluding an intermediate interface between the fibre and the matrix. Additionally, the effect of fibre spacing on fracture toughness is examined. Finally, Section 4.3 compares the two models proposed in Section 3.2, the single-fibre model and the multi-fibre model, with results compared to existing experimental data in the literature in Section 4.4.

## 4.1 The influence of the size of the fibre

In this section, we examine the influence of two variables: the length $l$ and the diameter $t$ of the fibre. The analysis is conducted for **Material 2**, as described in Section 3.1, using the multi-fibres model detailed in Figure 4. Only two possible paths for the propagation of the pre-existing crack, $\Gamma_c$, are analysed, as outlined in Section 3.2: crack penetration into the fibre, and lateral deflection through the interphase. Examples of the two crack paths are represented below in Figures 10 and 12. These are referred to as toughening mechanisms, as they are expected to enhance the fracture toughness of the composite compared to that of the matrix alone. The predominant mechanism — i.e., the one that governs crack propagation — will be the one associated with the lowest calculated value of the apparent fracture toughness $K_{\text{IC}}^{\text{app}}$.

Figure 6 shows the evolution of the apparent fracture toughness, $K_{\text{IC}}^{\text{app}}$, with respect to $t$, where the chosen range $t = 3 - 54 \ \mu\text{m}$ is based on the data provided in Section 3.1. For very thin fibres, it can be observed that penetration governs the failure, whereas for thicker fibres, the predominant mechanism is deflection. This suggests that when the fibre is thicker, the interphase acts as a mechanical fuse, causing the crack to deflect through the interphase. Note that a lower $K_{\text{IC}}^{\text{app}}$ is directly associated with a lower failure load, and therefore, the dominant failure mechanism is the one with the lowest $K_{\text{IC}}^{\text{app}}$. The transition from penetration to deflection occurs for $t$ in the range [15 – 25] $\mu\text{m}$. This explains the experimental observations of fibre deflection when the interphase is present, as the fibre diameter is sufficiently large to favour this mechanism (see Yasmin Legerstee 2022). Furthermore, a larger fibre diameter results in a higher fracture toughness for the composite, with the maximum enhancement of $K_{\text{IC}}^{\text{app}} = 1.38 K_{\text{IC}}^{\text{m}}$ occurring at $t = 54 \ \mu\text{m}$ (the maximum diameter in the range considered), corresponding to $K_{\text{IC}}^{\text{app}} = 3.34 \ \text{MPa m}^{0.5}$.

On the other hand, Figures 7 and 8 show the evolution of $K_{\text{IC}}^{\text{app}}$ with respect to the fibre length, $l$, for two different diameters: $t = t_{\text{avg}} = 16.67 \ \mu\text{m}$ and $t = 25 \ \mu\text{m}$. The range selected for $l$ is based on the data provided in Section 3.1. First, it can be observed that $l$ does not change the predominant mechanism, which remains penetration for $t = t_{\text{avg}} = 16.67 \ \mu\text{m}$ and deflection for $t = 25 \ \mu\text{m}$. Therefore, it can be concluded that the predominant mechanism mainly depends on the fibre diameter. However, $l$ does influence the apparent fracture toughness, with higher values observed when $l \approx t$. Moreover, Figures 7 and 8 show that for very long fibres, the influence of $l$ on $K_{\text{IC}}^{\text{app}}$ is significantly reduced. Interestingly, in the case of $t = 25 \ \mu\text{m}$, the apparent fracture toughness for penetration and deflection mechanisms becomes quite similar when considering
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Effect of the fibre diameter *t* on the apparent fracture toughness of the composite $K_{IC}^{app}$ for a fibre volume fraction $V_f$ = 10% and a fibre length *l* = 152 µm.</figcaption></figure>

very short fibres. It is worth noting that the slight non-monotonic variation of $K_{IC}^{app}$ is due to the
reduced scale of the plot and remains below 3%, likely caused by minor numerical fluctuations.

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Effect of the fibre length *l* on the apparent fracture toughness of the composite $K_{IC}^{app}$ for a fibre volume fraction $V_f$ = 10% and a fibre diameter *t* = 16.67 µm.</figcaption></figure>

In conclusion, the fibre diameter plays a major role in determining the predominant toughening mechanism observed in the composite, although both *l* and *t* influence $K_{IC}^{app}$, which increases when *l* is close to *t*. Consequently, square fibres represent the most effective configuration for increasing $K_{IC}^{app}$. When *l* ≈ *t*, they become more robust and resistant to penetration, resulting in the highest values of $K_{IC}^{app}$. For this reason, Figure 9 represents the evolution of $K_{IC}^{app}$ with respect to the fibre size for *t* = *l*. First, it can be observed that $K_{IC}^{app}$ increases with increasing *l*. Moreover, for small fibres, the predominant mechanism is penetration, whereas for larger fibres, the crack may deflect through the interphase. It is also evident that $K_{IC}^{app}$ becomes independent of *t* for very small fibres in the deflection mechanism. This is because the energy condition of the Coupled Criterion (CC), which does not depend on the fibre size, governs the failure. An example is shown in Figures 11(a) and 12(a), where the CC point is governed by the minimum of the energy condition for a small fibre and the cross point between the stress and the energy condition for a large fibre.
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Effect of the length of the fibre *l* on the apparent fracture toughness of the composite $K_{IC}^{app}$ for a volume fraction of fibres $V_f = 10\%$ and a thickness of the fibre $t = 25$ µm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Effect of the fibre length and diameter on the apparent fracture toughness of the composite $K_{IC}^{app}$ for a fibre volume fraction $V_f = 10\%$ when $l = t$.</figcaption></figure>

To illustrate the case of deflection for small fibre lengths, an example of the application of the CC in these cases is shown in Figure 10, where the curves of the energy and stress conditions are presented, along with the deformed shape after failure. Note that the crack increment (referred to as the nucleation length in the framework of Finite Fracture Mechanics) is shorter than the length of the interphase, explaining why the crack is arrested before reaching the fibre corner.

If $t > 25$ µm in Figure 9, the apparent fracture toughness of the composite, $K_{IC}^{app}$, begins to increase with growing $t$ in both the deflection and penetration cases. Both situations correspond to conventional cases, where the CC point is the intersection of the energy condition and the stress condition, as shown in Figures 11 and 12 for deflection and penetration, respectively. Notice that, due to the small newly created crack length, we have represented only the vicinity of one fibre in the deformed shape, although the multi-fibre model has been applied.

## 4.2 The influence of other design parameters

In this section, the influence of other design parameters is studied, such as the fibre volume fraction, the interface effect, and the distance between fibres. The **Material 2** and the multi-fibre
<figure><img src="image_11.png" /><figcaption>**Figure 10:** Application of the Coupled Criterion for the case of deflection in the optimal configuration $l = t$, for the volume fraction of fibres $V_f = 10\%$, and $l = t = 8$ µm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Application of the Coupled Criterion for the case of deflection in the optimal configuration $l = t$, for a fibre volume fraction $V_f = 10\%$, and $l = t = 100$ µm.</figcaption></figure>

model are used unless otherwise stated. Throughout this section, a fibre with $t = l = 17$ µm is considered as the object of study. First, Figure 13 illustrates the influence of $V_f$ on $K_{IC}^{app}$. The values $V_f = 5, 10,$ and $20\%$ correspond to the range obtained statistically in Yasmin Legerstee 2022. In addition, the value $V_f = 40\%$ has been proposed to study the effect of saturating the composite with more fibres, although it is purely a numerical model, as it is unknown whether such a high fibre volume fraction is feasible in practice. For $V_f = 5\% - 20\%$, there is little variation in $K_{IC}^{app}$, with a maximum increase with respect to $K_{IC}^m$ of 35% for penetration and 49% for deflection. However, when considering a higher fibre volume fraction, such as $V_f = 40\%$, the fracture toughness of the composite increases 46% for penetration and 52% for deflection with respect to $K_{IC}^m$.

The reason of the high volume fraction effect can be explained by the Figure 14 illustrating the geometry of the multi-fibre model for different values of $V_f$. As the volume fraction increases, the distance between the fibres decreases, that can make extra interactions between fibres allowing higher fracture toughness.

The influence of the interphase in the composite is studied and depicted in Figure 15, where the apparent fracture toughness of the composite, $K_{IC}^{app}$, is shown as a function of the fibre diameter $t$ for both penetration and deflection when **Material 1** is considered, i.e., when a strong interface between the matrix and the fibre is present. It can be observed that $K_{IC}^{app}$ in the case
<figure><img src="image_13.png" /><figcaption>**Figure 12:** Application of the Coupled Criterion for the case of penetration in the optimal configuration $l = t$, for a fibre volume fraction $V_f = 10\%$, and $l = t = 100$ µm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Effect of the fibre volume fraction $V_f$ on the apparent fracture toughness of the composite $K_{IC}^{app}$ for $t = l = 16.67$ µm.</figcaption></figure>

of deflection increases significantly compared to the case where the interphase is included, making penetration the predominant mechanism, regardless of the fibre diameter. This explains the experimental observation that when the interphase is removed, penetration becomes the dominant damage mechanism, regardless of fibre size (see Yasmin Legerstee 2022).

Figure 16 shows $K_{IC}^{app}$ for the cases with (Material 2) and without interphase (Material 1). For small-diameter fibres, where the predominant mechanism is penetration, there is no difference between the interphase and non-interphase cases. However, for larger fibre diameters, $K_{IC}^{app}$ is greater in the non-interphase case because the interphase causes a crack propagation by deflection easier. Note that this observation appears to contradict experimental measurements from Yasmin Legerstee 2022, as shown in Table 2. This is because the measurement of $K_{IC}^{app}$ in experiments is based on the total failure of the structure, which is related to crack propagation until complete failure occurs. When the interphase between the fibre and matrix is included, the crack meanders and delays the moment of total failure, thereby increasing the overall load. This effect cannot be captured by our numerical tool, which only measures the interaction between the crack and a single fibre as a necessary criterion for initiation of the crack propagation.

## 4.3 Multiple fibres vs single fibre model

In Figures 17 and 19, $K_{IC}^{app}$ is shown as a function of $t$ for both the multi-fibre and single-fibre models, for two different values of $V_f = 10\%$ and $40\%$, considering only the deflection case (similar conclusions are obtained for the penetration case, but these have been omitted for simplicity). For $V_f = 40\%$, the maximum difference reaches 4.15%, which is greater than for $V_f = 10\%$, where the
<figure><img src="image_15.png" /><figcaption>**Figure 14:** Geometries of the multi-fibre model for different fibre volume fraction values.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15:** Effect of the fibre diameter *t* on the apparent fracture toughness of the composite, $K_{IC}^{app}$, for a fibre volume fraction of $V_f$ = 10% and a fibre length of *l* = 152 µm, considering no interphase.</figcaption></figure>

maximum difference is 1.5%. In other words, the larger the volume fraction, the greater the difference between the two models, although the difference remains relatively small. As observed in Fig. 17, the distinction between the multi fibres and single fibre models becomes significant for smaller fibre sizes, where the influence of neighbouring fibres - accounted for only in the multi-fibre model - is increasingly more pronounced.

In general, the difference between the multi-fibre model and the single-fibre model will be greater when the neighbouring fibres have a more significant effect on the central fibre, which is the one directly interacting with the crack. For the same $V_f$, a smaller distance between fibres in the multi-fibre model ($h_l$ and $h_t$) will increase the influence of the neighbouring fibres, thus amplifying the difference between the two models. This is clearly illustrated in Figure 18, for $V_f$ = 10% and $l$ = 152 µm, where $h_l$ and $h_t$ are reduced for smaller values of $t$. The maximum difference between the models is observed in this region of smaller $t$.

However, this is not observed in Figure 19, where it can be seen that the difference between the two models increases with growing $t$ for $V_f$ = 40%. Looking at Figure 20, where $h_l$ and $h_t$ are represented, it is observed that although the distance $h_l$ increases as $t$ increases, the distance $h_t$ remains almost constant except for the last point where it is slightly higher. This varies aspect ratio between the two distances representing variations of fibre stacking. Additionally, $h_t$ is zero. For this reason, the difference between the two models may be motivated by the distance between the fibres and their stacking format of fibres, which will be more pronounced if we consider more fibres.

Finally, in Figure 21, the influence of the distance $h_l$ is analysed. Previously, we used
<figure><img src="image_17.png" /><figcaption>**Figure 16:** Comparison of $K_{IC}^{app}$ with respect to $t$ with and without the interphase for $V_f = 10\%$ and $l = t$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17:** Comparison between the single fibre and the multi-fibre model in the deflection case, for $l = l_{avg} = 152 \ \mu m$. Case of $V_f = 10\%$.</figcaption></figure>

the minimisation of surface area as a criterion to calculate $h_l$ and $h_t$ in the heterogeneous microstructure of the inner domain. However, this criterion lacks a physical explanation and may not be representative. For this reason, the influence of the distance between fibres is studied by selecting a wide range of $h_l$ values and defining $h_t$ so that the constraints in Equations (13)–(16) are satisfied. It is observed that for negative values of this distance, which would imply that fibres are interleaved in the microstructure, there is a jump in the fracture toughness. The improvement in fracture toughness due to an interleaved fibre structure is not surprising. For instance, in construction, brick structures often utilise a technique known as "stretcher bond", where the bricks are interleaved to prevent the joints from aligning, thereby enhancing the strength of the wall. It is important to note that this effect can only be observed when considering a multi-fibre model and not with the previously presented single fibre model.

Two examples of the geometries obtained for $h_l > 0$ and $h_l > 0$ are shown in Figures 22 and 23. The parameter $h_t$ is obtained to fulfill the geometrical restrictions in Equations (13)–(16).
<figure><img src="image_19.png" /><figcaption>Figure 18: Evolution of distances $h_l$ and $h_t$ for $l = l_{avg} = 152$ µm. Case of $V_f = 10\%$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison between the single fibre and the multi-fibre model in the deflection case, for $l = l_{avg} = 152$ µm. Case of $V_f = 40\%$.</figcaption></figure>

## 4.4 Comparison to experiments

Finally, it is challenging to validate this model through a precise comparison with the experimental results shown in Table 2, as the specimens contain fibres of various sizes and aspect ratios. Nevertheless, the range numerically estimated for $K_{IC}^{app}$, considering the average $l_{avg}$, $V_f = 10\%$, and $t = 3 - 54$ µm, is

$$K_{IC}^{app} = [2.34 - 4.50] \text{ MPa m}^{0.5}. \qquad (27)$$

with the interphase,

$$K_{IC}^{app} = [2.34 - 3.34] \text{ MPa m}^{0.5} \qquad (28)$$

without the interphase. The results are similar to the range obtained experimentally without an interphase, $K_{IC}^{app} = [2.80 - 4.40] \text{ MPa m}^{0.5}$, although it differs more from the range obtained in Yasmin Legerstee 2022 with the interphase, $K_{IC}^{app} = [3.87 - 6.32] \text{ MPa m}^{0.5}$. This discrepancy may be attributed to several factors, including the fact that dynamic effects of the interface were not considered here in our numerical approach, and some of the mechanical properties used for the interphase were those of BN instead of PyC, which is the material used in the experiments. As
<figure><img src="image_21.png" /><figcaption>**Figure 20:** Evolution of distances $h_l$ and $h_t$ for $l = l_{avg} = 152$ µm. Case of $V_f = 40\%$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21:** Influence of the distance between fibres $h_l$ on the apparent fracture toughness of the composite $K_{IC}^{app}$ for $V_f = 10\%$ and $l = t = 17$ µm, considering the interphase (Material 2).</figcaption></figure>

mentioned in Section 3.1, the properties of BN are the ones mentioned in Yasmin Legerstee 2022 to describe the interphase for confidential reasons. Moreover, 3D effects and stochastic features (fibre orientation, materials properties, defects) may lead to quantitative discrepancies. It is important to highlight that only crack initiation from the pre-existing crack is considered here, not the crack propagation across a distance larger to crack initiation length.

However, qualitative conclusions obtained in experiments are also observed using this computational tool. When considering the interphase, a very small diameter of the fibre results in lower $K_{IC}^{app}$, and the failure mechanism is penetration. However, when fibres have larger $t$, they are more robust and difficult to penetrate, resulting in higher $K_{IC}^{app}$, and the failure mechanism is deflection, allowing the interface to act as a mechanical fuse. This would explain why experimental observations show fibre deflection when an interphase is present, since the diameter used is sufficiently large to favor this mechanism, according to our numerical estimations. In case of the system without interface, penetration becomes the dominant damage mechanism, regardless of fibre size. This numerical result is also confirmed experimentally.

Note that, unlike the penetration mechanism, the deflection mechanism causes a much longer crack propagation path, which requires additional energy for the increasing path to failure.
<figure><img src="image_23.png" /><figcaption>**Figure 22:** Fibre distribution in the inner domain with $h_l > 0$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 23:** Fibre distribution in the inner domain with $h_l < 0$</figcaption></figure>

Consequently, this must be taken into account for an accurate estimation of $K_{IC}^{\text{app}}$ based on the deflection mechanism. However, our current numerical model is not yet able to consider this effect. This could be the main reason for the small discrepancy observed without an *interphase* and the larger discrepancy with an interface between the numerical estimations and experimental results.

# 5 Conclusions

In this paper, a computational tool for estimating fracture toughness in short-fibre ceramic matrix composites has been developed and applied to an industrial composite material. The following key findings are highlighted:

* The fracture toughness estimated by the computational tool was thoroughly validated by comparing them with experimental data, both qualitatively and quantitatively.

* The fibre volume fraction effect on the $K_{IC}^{\text{app}}$ can be negligible when lower values up to $V_f = 20\%$, are considered, while a very high value of volume fraction, such as $V_f = 40\%$, increases the $K_{IC}^{\text{app}}$ significantly. It is because of the effect of interactions between fibres.

* The predominant failure mechanism is primarily influenced by the presence of the interphase and the fibre size. When the interphase is present, smaller fibre diameters result in lower $K_{IC}^{\text{app}}$, with penetration being the dominant failure mechanism. In contrast, larger fibre diameters promote fibre deflection, allowing the interphase to act as a mechanical fuse, which is same as experimental observation. Conversely, when the interphase is removed, penetration becomes the dominant damage mechanism, regardless of fibre size, as observed experimentally. Notice that this computational tool could be used to determine the properties of the interphase that govern the competition between penetration and deflection.

* There is not a significant difference between the multi-fibre and the single fibre model if there is no change of fibre stacking represented by the ratio between the two distances between fibres, $h_l$ and $h_t$. That is why, for general estimations, the use of the single fibre model is proposed as it has a lower computational complexity. For the same level of refinement, it requires approximately half of the number of elements.

* Square fibres represent the most effective configuration for increasing $K_{IC}^{\text{app}}$. When $l \approx t$, they become more robust and resistant to penetration, resulting in the highest values of $K_{IC}^{\text{app}}$. It is important to note that the spatial distribution of fibres, and therefore the
statistical probability of the crack encountering a fibre, has not been included in this analysis. The optimal configuration could change if this effect is considered. Moreover, an additional condition should be considered to ensure that the crack does not propagate exclusively through the matrix. This can be addressed by imposing that the probability of a crack encountering a fibre is greater than 50%, i.e., $3^l/\hat{a}_l > 0.5$.

* Additionally, the intercalation of fibres also has an effect on $K_{IC}^{app}$, which can increase approximately 4%.

* Finally, no initial stress field is considered, which is likely to originate from the processing, e.g., residual stresses. This could be included in the model as a future work. An example is shown in Jiménez-Alfaro and Leguillon 2023.

# 6 Appendix

Details about the numerical implementation

## Acknowledgements

S. Jimenez-Alfaro acknowledges the Iberdrola Foundation under the Marie Skłodowska-Curie Grant Agreement No 101034297. She also acknowledges the support of the Ministerio de Ciencia e Innovación of Spain for his postdoctoral fellowship under reference JDC2023-050492-I.

## References

Barsoum, M. W., P. Kangutkar, and A. Wang (1992). "Matrix crack initiation in ceramic matrix composites Part I: Experiments and test results". *Composites science and technology* 44.3, pp. 257-269

Dhanasekar, S., A. T. Ganesan, T. L. Rani, V. K. Vinjamuri, M. Nageswara Rao, E. Shankar, Dharamvir, P. S. Kumar, and W. Misganaw Golie (2022). "A Comprehensive Study of Ceramic Matrix Composites for Space Applications". *Advances in Materials Science and Engineering* 2022, p. 6160591

Doitrand, A., R. Henry, H. Saad, S. Deville, and S. Meille (2020). "Determination of interface fracture properties by micro- and macro-scale experiments in nacre-like alumina". *Journal of the Mechanics and Physics of Solids* 145, p. 104143

Donald, I. W. and P. W. Mcmi Llan (1976). "Ceramic-matrix composites". *Journal of Materials Science* 11, pp. 949-972

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

European Commission (2019). *The European Green Deal*. COM

Jiménez-Alfaro, S. and D. Leguillon (Dec. 2022). "Modelling of glass matrix composites by the Coupled Criterion and the Matched Asymptotics approach. The role of a single platelet". *Theoretical and Applied Fracture Mechanics* 122, p. 103650

Jiménez-Alfaro, S. and D. Leguillon (Dec. 2023). "Modelling of glass matrix composites by the Coupled Criterion and the Matched Asymptotic Approach. The effect of residual stresses and volume fraction". *Theoretical and Applied Fracture Mechanics* 128, p. 104112

Jiménez-Alfaro, S., I. G. García, and A. Doitrand (2025). "Review of the matched asymptotic approach of the coupled criterion". *Comptes Rendus - Mecanique* 353.G1, pp. 339-357

Jiménez-Alfaro, S. and D. Leguillon (2021). "Finite fracture Mechanics at the micro-scale. Application to bending tests of micro cantilever beams". *Engineering Fracture Mechanics* 258, p. 108012

Leguillon, D., E. Martin, O. Ševeček, and R. Bermejo (2015). "Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses". *European Journal of Mechanics, A/Solids* 54, pp. 94-104

Leguillon, D. (2002). "Strength or toughness? A criterion for crack onset at a notch". *European Journal of Mechanics A/Solids* 21, pp. 61-72

Leguillon, D., E. Martin, O. Sevecek, and R. Bermejo (2018). "What is the tensile strength of a ceramic to be used in numerical models for predicting crack initiation?" *International Journal of Fracture* 212.1, pp. 89-103
Martin, E., D. Leguillon, and C. Lacroix (2001). "A revisited criterion for crack deflection at an interface in a brittle bimaterial". *Composites Science and Technology* 61, pp. 1671-1679

Martin, E., D. Leguillon, O. Sevecek, and R. Bermejo (2018). "Understanding the tensile strength of ceramics in the presence of small critical flaws". *Engineering Fracture Mechanics* 201, pp. 167-175

Martin, E., P. W. Peters, D. Leguillon, and J. M. Quenisset (1998). "Conditions for matrix crack deflection at an interface in ceramic matrix composites". *Materials Science and Engineering: A* 250.2, pp. 291-302

Ohnabe, H., S. Masaki, M. Onozuka, K. Miyahara, and T. Sasa (1999). "Potential application of ceramic matrix composites to aero-engine components". *Composites Part A: Applied Science and Manufacturing* 30.4, pp. 489-496

Oriol Gavalda-Diaz, O., E. Saiz, J. Chevalier, and F. Bouville (2025). "Toughening of ceramics and ceramic composites through microstructure engineering: A review". *International Materials Review* 70.1, pp. 3-30

Papšík, R., O. Ševeček, A. K. Hofer, I. Kraleva, J. Kreith, and R. Bermejo (2023). "Prediction of edge and tunnelling crack formation in layered ceramics using a stress-energy fracture criterion". *Journal of the European Ceramic Society* 43.7, pp. 2928-2934

Papšík, R., O. Ševeček, J. Schlacher, and R. Bermejo (2024). "Prediction of thermal shock induced cracking in multi-material ceramics using a stress-energy criterion". *Engineering Fracture Mechanics* 303, p. 110121

Picard, D., D. Leguillon, and C. Putot (2006). "A method to estimate the influence of the notch-root radius on the fracture toughness measurement of ceramics". *Journal of the European Ceramic Society* 26.8, pp. 1421-1427

Ševeček, O., M. Kotoul, D. Leguillon, E. Martin, and R. Bermejo (2016). "Modelling of edge crack formation and propagation in ceramic laminates using the stress-energy coupled criterion". *Engineering Fracture Mechanics* 167, pp. 45-55

Weißgraeber, P., D. Leguillon, and W. Becker (Jan. 2016). "A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers". *Archive of Applied Mechanics* 86.1-2, pp. 375-401

Williams, M. L. (Dec. 1952). "Stress Singularities Resulting From Various Boundary Conditions in Angular Corners of Plates in Extension". *Journal of Applied Mechanics* 19.4, pp. 526-528

Yasmin Legerstee (2022). *Caractérisation et modélisation du comportement thermomécanique de composites à matrice céramique renforcés par des fibres courtes*. Université de Bordeaux
# Notes on revision made to manuscript JTCAM-15621

## Response to the reviewers

### Reviewer 1

The paper presents a numerical study about the influence of a short fiber on the propagation of an existing crack in a ceramic composite, based on a 2D modeling including either one isolated fiber, or a fiber with some surrounding fibers, the fibers being represented by rectangles. The influence of several geometry parameter on the apparent propagation SIF is investigated based on the coupled criterion (CC), including either penetration of the crack through the fiber or deflection at the fiber interface. The presence of an interphase is also considered. Some interesting results are obtained about the most favorable fracture mechanism that may occur, the influence of the volume fraction, or the influence of other surrounding fibers on the critical SIF at crack propagation.

We appreciate the positive opinion of the reviewer and the constructive feedback provided. In this document, dark blue reflects changes made to the manuscript, while light blue indicates our responses.

Some shortcomings are also to be improved before publication of the manuscript:

**Comment 1.1** — The presentation of the different fracture scenario could be better explained, with an illustration of the configurations before and after the small crack increment.

**Reply:** Thank you for your comment. To clarify this aspect, Figures 3 and 4 have been changed to better represent the initial undeformed shape. A comment has been included in the caption of both figures.

**Comment 1.2** — Regarding the presentation of the CC implementation. Details of the MA expansions are missing (it could be added as an appendix or in the main text), as well as details on the numerical implementation and mesh convergence (could be added in appendix too)

**Reply:** Thank you for the comment. More details have been included in Section 3.2.

"...As an example, the development of the inner expansion is written below considering the length of the fibre as the characteristic length in the adimensionalization, but a similar development can be obtained with the diameter of the fibre. The inner expansion is therefore defined as

$$U(x_1, x_2) = U(ly_1, ly_2) = F_0(l)V_0(y_1, y_2) + F_1(l)V_1(y_1, y_2) + \dots \quad (1)$$

The terms of Equation (1), and consequently the elastic solution in the vicinity of the pre-existing crack tip, $\Gamma_c$, are derived by relating the outer and inner expansions in a common region where both approximations accurately represent the actual solution. In that region, the so-called matching conditions are applied,

$$F_0(l)V_0(y_1, y_2) \approx U(0, 0), \quad \text{when } \rho \to \infty, \quad (2)$$

$$F_1(l)V_1(y_1, y_2) \approx K_1 l^{\lambda} \rho^{\lambda} u(\theta), \quad \text{when } \rho \to \infty, \quad (3)$$

where $\rho = \sqrt{y_1^2 + y_2^2}$ and $\approx$ means "behave like". Hence, it can be set

$$F_0(l) = 1, \quad V_0(y_1, y_2) \approx U(0, 0), \quad \text{when } \rho \to \infty, \quad (4)$$

$$F_1(l) = K_1 l^{\lambda}, \quad V_1(y_1, y_2) \approx \rho^{\lambda} u(\theta), \quad \text{when } \rho \to \infty, \quad (5)$$

The expression in Equation (1) is a good approximation of the real solution $U(x_1, x_2)$ in the vicinity of the crack tip. To obtain the solution, the finite element method is applied in a circular domain whose dimensions are much larger (in this case, $R = 400 \ \mu m$) than those of the internal domain described in Figures 3 and 4. The element size in the vicinity of the pre-existing crack tip $\Gamma_c$ is sufficiently small to ensure numerical convergence, with a size of 0.008 $\mu m$. On the boundary of the domain, a prescribed displacement is applied as a boundary condition, as shown in Figure 1..."

**Comment 1.3** — A comparison of the different fracture scenario is only given in terms of critical SIF, but unfortunately no further effort is provided to explain why one or the other mechanism is predominant.
<figure><img src="image_25.png" /><figcaption>**Figure 1.** Scheme of the numerical problem.</figcaption></figure>

However, it is the main advantage of the CC that it is possible to explain why the transition between the different mechanisms occurs by studying the stress and the energy criteria. As a consequence, it would be worth detailing more the CC solution before presenting the results in terms of critical SIF to evidence why a specific fracture mechanism occurs rather than the other.

**Reply:** Thank you for the comment. The Stress Intensity Factor (SIF) is crucial in determining the dominant failure mechanism. This is because the SIF is directly proportional to the applied load. Therefore, the lowest SIF corresponds to the mechanism that becomes active at the lowest failure load, and is thus considered the dominant one. A comment regarding this has been included in Section 4.1: "Note that a lower $K_{IC}^{app}$ is directly associated with a lower failure load, and therefore, the dominant failure mechanism is the one with the lowest $K_{IC}^{app}$..."

**Comment 1.4** — The proposed model lacks some simple physical considerations, for instance a main result is that square fibers should be used for the best improvement of the composite toughness. However, with square fibers and low fiber volume fractions, it is unlikely that a crack encounters a fiber, it would rather propagate into the matrix between the fiber, thus not resulting in any toughening. The analysis should include such considerations to provide more physically based results.

**Reply:** Thank you very much for the comment. We agree that the probability of the pre-existing crack not encountering a fiber has not been included in the study, which assumes that the crack has already intersected a fiber. A remark has been added in Section 5, the Conclusions: "It is important to note that the spatial distribution of fibers, and therefore the statistical probability of the crack encountering a fiber, has not been included in this analysis. The optimal configuration could change if this effect is considered...."

**Comment 1.5** — The modeling of crack deflection through the interphase may be questioned, regarding i) its location (within the interphase or at the fiber/interphase interface), ii) the fracture properties employed (which are those of the interphase) and iii) the choice to neglect the influence of mode mixity on the stress criterion and the fracture properties.

**Reply:** Thank you for your contribution. (i) The pre-existing crack reaches the interphase, which has been more clearly highlighted in Figures 3 and 4. The deflection of the crack occurs at the interface. (ii) The fracture properties of the interface considered are provided in Table 1. If these properties were to change, the results of the analysis would also be affected accordingly. (iii) In this case, crack propagation is not considered — only crack nucleation is addressed. To determine this nucleation, we use a technique that avoids explicitly including mixed-mode fracture criteria. This is achieved by relating the solution near the crack tip in the heterogeneous domains (highlighted in Figures 3 and 4) to the solution near the crack tip assuming no fibers are present — what we refer to as the outer solution. The latter corresponds to a Mode I crack.

Here are some comments/questions to help the authors improve the current manuscript:

**Comment 1.6** — The introduction could be improved, it is basically describing CMC first, then describing the CC, then saying that the CC is applied to CMC. It seems that a clear motivation lacks, as well as a description of the literature about CMC cracking which is not investigated at all, even if many experimental and numerical works focus on the damage mechanisms of these materials, at different scales. Please improve the introduction to better motivate this work and situate this work with respect
to existing literature about fracture analysis in CMCs. The reference Martin et al., 1998 is relevant to this work but it is not an application of the CC, as mentioned in the text

**Reply:** Thank you very much for your contribution. More literatures about fracture analysis in CMCs including the reference Martin et al., 1998 were introduced in introduction as below.

In addition, it was established that the fracture toughness of the CMCs can be enhanced by introducing weak interfaces between the fibre and the matrix. Then Martin, Peters, et al. 1998; Martin, Leguillon, and Lacroix 2001 considered crack deflection mechanism at the interface between fibre and matrix as a principal toughening mechanism of the CMCs. This mechanism of microscale cracking in the CMCs was further investigated by many other researchers. (See e Oriol Gavalda-Diaz et al. 2025 for a recent critical review)"...

Other experimental and theoretical studies have been applied for explaining the crack path in ceramic matrix composites. Martin, Peters, et al. 1998 applied an energy criterion for predicting a crack deflection at the fibre-matrix interface. Barsoum et al. 1992 investigated the relationship between microstructural factors and matrix cracking initiation stress in fibre-reinforced ceramic matrix composites. However, the exact process by which these evolve into a macroscopic, fibre-bridged crack remained unresolved."...

**Comment 1.7** — Equation 3 - I would recommend a different naming of the different quantities: The overline is generally used for a mean quantity, so it would be adapted for $\bar{G}_c$ whereas $G_{inc}$ keeps its usual definition so I would keep $G_{inc}$ (without overline) and put $\bar{G}_c$ (with overline)

**Reply:** Thank you for the remark. It has been changed:

"The energy condition for a newly created crack length $\delta l$ is written in a heterogeneous structure as

$$G_{inc} = -\frac{\Delta\Pi_p(\delta l)}{\delta l}, \qquad (6)$$

and therefore

$$\frac{G_{inc}}{\bar{G}_c(\delta l)} \geq 1, \qquad (7)$$

where $\Delta\Pi_p(\delta l)$ is the difference in potential energy between the initial state prior to any crack extension and the final state when a crack of length $\delta l$ is generated $\Delta\Pi_p(\delta l) = \Pi_p(\delta l) - \Pi_p(0)$. The function $\bar{G}_c(\delta l)$ describes the critical energy release rate along the expected crack path,

$$\bar{G}_c(\delta l) = \frac{\sum_{i=1}^{N} G_c^{(i)} \delta l^{(i)}}{\delta l}. \qquad (8)$$

being $\delta l^{(i)}$ the crack length associated with material (i)..."

**Comment 1.8** — When referring to “the tensile stress”, it would be worth providing an axis system in figure 1 and detail which stress component is indeed involved.

**Reply:** Thank you for your comment. An additional axis representing the normal direction (n) has been added together with the axis s, that represents the crack path. Figure 1 has been changed and a new comment is included in the text to clarify the direction of the tensile stress: "i.e., the stress component in the normal direction n highlighted in Fig. 1"...

**Comment 1.9** — In CMCs, and more generally in composite materials, the constituents may be anisotropic, which is disregarded in the present analysis. Could the authors comment on the influence of considering the material anisotropy on the results?

**Reply:** Thank you for your comment. In Yasmin Legerstee 2022, a study on the sensitivity of the elastic properties in the three spatial directions is conducted, concluding that the differences are not significant, and thus considering anisotropy is not necessary for this material. However, a comment has been included in Section 3.1: It is worth noting that in Yasmin Legerstee 2022, a study on the sensitivity of the elastic properties in the three spatial directions is conducted, concluding that the differences are not significant, and thus considering anisotropy is not necessary for this material.

**Comment 1.10** — Table 2: the authors provide two values for the critical stress intensity factor, name $K_{LEFM}$ or $K_{DIC}$, but no indications are given on how they were respectively determined and the difference between each other.
**Reply:** Thank you for your comment. The fracture toughness of the composite is estimated in Yasmin Legerstee 2022 using an analytical approximation based on Linear Elastic Fracture Mechanics (LEFM) and experimentally determined through Digital Image Correlation (DIC). Figure 2 has been changed and some comments have been included:

In the first one, a shape factor $f(\frac{a}{W})$ is used,

$$K_{\text{LEFM}} = f\left(\frac{a}{W}\right) \sigma_c \sqrt{\pi a}, \qquad (9)$$

where a and W are the notch and the specimen width, respectively (see Figure ??). The tensile strength $\sigma_c$ of the composite is obtained considering the critical force of the test $P_c$:

$$\sigma_c = \frac{3P_c(L_e - L_i)}{2eW^2} \qquad (10)$$

being e the thickness of the specimen, and $L_e, L_i$ the position of the boundary conditions (see Figure ??). The shape factor $f(\frac{a}{W})$ is defined as

$$f\left(\frac{a}{W}\right) = \frac{1.1215}{\sqrt{(1 - \frac{a}{W})^3}} \left[ \frac{5}{8} - \frac{5}{12}\frac{a}{W} + \frac{1}{8}\left(\frac{a}{W}\right)^2 + 5\left(\frac{a}{W}\right)^2\left(1 - \frac{a}{W}\right)^6 + \frac{a}{W}\exp\left(-\frac{6.1342\frac{a}{W}}{1 - \frac{a}{W}}\right) \right], \quad (11)$$

The second methodology is based on the well-known Williams's expansion in Williams 1952, that defines the elastic solution $u(r, \theta) = K_I r^{0.5} u_I(\theta)$ in the neighbourhood of the Mode I crack tip. The DIC technique is applied to obtain the displacement field at the moment of failure, and therefore $K_c = K_I$.

**Comment 1.11** — Figure 4 Interaction between several fibers. The parameter $h_t$ represents the distance between the middle and upper fiber, is it also the distance between the middle and lower fibers? It seems to be from equation (8) but it is misleading as the middle fiber does not seem to be vertically centered on Figure 4. Equation (10) does $h_t$, $h_t > 0$ means $h_t > 0$ and $h_t > 0$? why is there a “if” at the beginning of eq (10)?

**Reply:** Thank you for your remark. Yes, it is also the distance between the middle and lower fibers. The figure has been changed to represent it accordingly. There was a typo in the equation that has been fixed.

**Comment 1.12** — From Equation (14) to equation (15-16). The asymptotic expansions are completely overlooked. Does the same expansions hold whatever the damage scenario? For instance, regarding the choice of the stress criterion. It would be worth adding some details (in the text, or in appendix) about this point

**Reply:** Thank you for your contribution. Some additional comments about the point have been included. In Equation 23 the tensile stress $\sigma(\hat{s})$ is obtained as

$$\sigma(\hat{s}) = \frac{1}{\sqrt{l}} \nabla_y V_1(y_1, y_2) \qquad (12)$$

if the inner domain is built considering the length of the fibre l or

$$\sigma(\hat{s}) = \frac{1}{\sqrt{t}} \nabla_y V_1(y_1, y_2) \qquad (13)$$

if the inner domain is built considering the diameter of the fibre t. The term $\nabla_y$ is the gradient operator with respect to the variables $y_1, y_2$. As explained in Jimenez Alfaro and Leguillon 2022 the increment in potential energy $\Delta\Pi_p$ is the one directly obtained from the inner problem (solution $V_1(y_1, y_2)$).

**Comment 1.13** — The studied fracture scenarios are not clearly described. It seems from figures 3 and 4 that the initial crack goes through the interphase so that the initial crack tip is located on the Fiber? But the configurations used to compute $\delta\Pi$ are not well described in the sequel. Actually, in Section 4.1, it is given: “crack penetration into the interphase and the fibre, and lateral deflection through the interphase”. For the former, it is not consistent with figures 3 and 4. Does symmetric or asymmetric lateral deflection is considered? How the choice is made? Overall, it is not straightforward to clearly understand the configuration before and after the small crack extension. Please improve these descriptions with schemes and explanation in the text.
**Reply:** Thank you for the comment. Figures 3 and 4 have been changed to better describe the initial scenario. Only the asymmetric case is considered, based on the experimental observations in Legerstee (2022). Consequently, the sentence describing penetration has been addressed with the better description about the initial scenario so that it is more clear: "Examples of the two crack paths are represented below in Figures 10 and 12..."

**Comment 1.14** — Regarding the interface deflection mechanism, it is probably not the same stress component that is involved in the stress criterion? Interface cracking also means that shear stress may be involved. The formulation should be detailed according to the different studied damage scenario. Also, what properties are used for the interface? It can be expected that they are different from the Interphase, Matrix and fiber properties.

**Reply:** Thank you for your comment. It is the same formulation, since $\sigma$ denotes the tensile stress that is the one responsible for the crack opening in both the deflection and the penetration scenarii. The interface between the interphase, the matrix and the fibre is considered strong, the crack is expected to deflect inside the interphase or inside the fibre. Figures 3 and 4 are now better explained. A comment has been included when presenting Table 1: "Notice that the interface between the interphase, the fibre and the matrix is considered strong..."

**Comment 1.15** — No details are given on the “homogenized material”, what properties are taken for this material? What homogenization scheme is used to get these properties, and how is it ensured that the homogenized material has the same properties than the local arrangement of fiber/matrix in the inner domain? The transition between the platelet/matrix area and the surrounding homogeneous material induces spurious stress singularities. How are these stress singularity melt with? What is their influence on the local stress field and energy release due to crack propagation?

**Reply:** Thank you for your remark. The Voigt rule has been used to obtain the properties of the homogenized material. A comment has been included: The material properties of the composite are obtained applying the Voigt rule, and depend on the volume fraction.

Thank you for asking about possible spurious stress singularities. We apply a method that avoids explicitly including them. This is achieved by relating the solution near the crack tip in the heterogeneous domains to the solution near the crack tip in the so-called outer expansion, that considers an homogeneous material.

**Comment 1.16** — Figure 5 shows the competition between deflection through the interphase or penetration through the interphase and fiber. Penetration becomes more favorable than deflection for thin fibers. This result is not intuitive as a size effect could be expected so that the thin the fiber, the lower the dissipated energy due to crack nucleation through it. It would be worth detailing the stress and energy criteria for a) a thin ($t < 10$) and b) a large ($t > 30 \mu m$) fiber in order to better understand why one or the other mechanisms takes over. In particular, it is essential to elaborate on the choice of the criteria: is mode mixity considered when studying interface deflection. If not, could the authors provide an estimate of how the results would be influenced by considering a mixed mode stress criterion on the one hand, and a mixed-mode dependent critical energy release rate on the other hand?

**Reply:** Thank you for the remark. The mixed mode was mentioned in the comment 1.5. In this case, crack propagation is not considered — only crack nucleation is addressed. To determine this nucleation, we use a technique that avoids explicitly including mixed-mode fracture criteria. This is achieved by relating the solution near the crack tip in the heterogeneous domains (highlighted in Figures 3 and 4) to the solution near the crack tip assuming no fibers are present — what we refer to as the outer solution. The latter corresponds to a Mode I crack. The energy and stress criteria are provided for a thin and a large fibre in Figures 10 and 11. The predominant mechanism is chosen as the one associated with the lowest $K_I$, that is related to the lower failure load ($K_I$ is proportional to the load in the William's expansion).

**Comment 1.17** — In figure 6 and 7, how do the authors explain the non-monotonic variation of the critical SIF as a function of the platelet length in the case of penetration? Is it related to a physical mechanism (to be explained) or does it come from a numerical error (to be estimated)?

**Reply:** Thank you for your comment. The three last SIFs in the penetration curve (where the non-monotonic variation is observed) in figure 7 are 1.21, 1.23 and 1.22. It means that probably there is a numerical error, since there is not a change in the physical mechanisms that conduct the failure. This error would be less 3%.
A comment has been included: "...It is worth noting that the slight non-monotonic variation of $K_{IC}^{app}$ is due to the reduced scale of the plot and remains below 3%, likely caused by minor numerical fluctuations...."

**Comment 1.18** — “Consequently, square fibres represent the most effective configuration for increasing $K_{IC}^{app}$ → 10% of square fibers, the probability that a crack encounters a square fiber is smaller than with elongated fibers. As a consequence, it is most likely that the cracks go around such a fiber and finally only propagates in the matrix, thus not beneficiating from the toughening mechanism. Could the author revisit this result by considering the additional criterion that the probability that a vertical crack encounters a fiber is larger than 0.5, e.g. expressing that $3^l/\hat{a}_l > 0.5$, so that if a crack propagates in the matrix, there is no longer a toughening mechanism

**Reply:** We would like to thank the reviewer for the valuable comment. We agree that the effectiveness of the toughening mechanism depends not only on fibre geometry but also on the probability of crack-fibre interaction. Following the suggestion, we have included an additional condition in the manuscript to account for this effect. Moreover, an additional condition should be considered to ensure that the crack does not propagate exclusively through the matrix. This can be addressed by imposing that the probability of a crack encountering a fibre is greater than 50%, i.e., $3^l/\hat{a}_l > 0.5$.

**Comment 1.19** — Figure 8: “It is also evident that $K_{IC}^{app}$ becomes independent of $t$ for very small fibres in the deflection mechanism. This is because the energy condition of the Coupled Criterion (CC), which does not depend on the fibre size, governs the failure” → it is not evident if a more detailed analysis on the stress and energy condition is not provided. Please show them for small ($l < 10$ µm) and large ($l > 50$ µm) fiber length

**Reply:** Thank you for your comment. Figures 9 and 10 show the energy condition for a small fibre $l < 10$ µm and a large fibre $l > 50$ µm. The minimum point of the energy condition is the CC point in the case of a small fibre, indicating that it is governed by the energy condition. For a large fibre, the CC point is the cross point between the two conditions. A small comment has been added: "..."It is also evident that $K_{IC}^{app}$ becomes independent of $t$ for very small fibres in the deflection mechanism. An example is shown in Figures 11 and 12, where the CC point is governed by the minimum of the energy condition for a small fibre and the cross point between the stress and the energy condition for a large fibre...."

**Comment 1.20** — From figure 10 it appears that the crack is located at the interface between the fiber and the interphase. This interface is expected not to have the same fracture properties as the interphase or the fiber. However such properties are not given, only those of the interphase are. Why do they authors consider this mechanism instead of cracking within the interphase ?

**Reply:** Thank you for your remark. There was a problem in the Figure. It has been corrected now. The crack is growing in the interphase.

**Comment 1.21** — Figure 9: the crack path is hardly visible.

**Reply:** Thank you for your remark. It has been checked so that now the crack path can be observed clearly.

**Comment 1.22** — Figures 9, 10 and 11: their appear to be some small “oscillations” in the stress variations. Could the author elaborate on the possible origin of these variations? Did they ensure that the mesh used are indeed converged?

**Reply:** Thank you for your comment. These oscillations correspond to a change in the position from the fibre, to the interphase to the matrix. The blue curve is not exactly the stress, it is the stress divided by $\sigma_c$, that depends on the region where the stress is evaluated.

**Comment 1.23** — Page 11: “However, when considering a higher fibre volume fraction, such as $V_f = 40\%$, the fracture toughness of the composite increases by a factor of 10 for penetration and 3 for deflection, reaching 8.9% for penetration and 5.37% for deflection.” What is this 10 and 3 factors? It seems misleading as for $V_f = 40\%$, $K_{IC}^{app}/K_{IC}^m \approx 1.4 - 1.5$

**Reply:** Thank you for the remark. The message was really confusing. The comment has been changed: $K_{IC}^m$ of 35% for penetration and 49% for deflection. However, when considering a higher fibre volume fraction, such as $V_f = 40\%$, the fracture toughness of the composite increases 46% for penetration and 52% for deflection with respect to $K_{IC}^m$.
**Comment 1.24** — Figure 14: Could you please provide the analysis for the stress and energy criteria in presence or not of the interphase, to better understand why the presence of the interphase results in a change of mechanisms depending of the fiber diameter?

**Reply:** Thank you for the comment. The evolution of the stress and the change in potential energy are the same for both cases—with or without an interphase—that is, they correspond to those shown in Figures 10 and 11. The only difference is that $G_c$ and $\sigma_c$ correspond to those of the interphase, as indicated in Table 1, whereas in the case without an interphase, they are those of the matrix, which are higher. Therefore, in the presence of the interphase, since the fracture properties are lower for the same amounts of stress and energy, failure occurs under a lower failure load, that is, for a lower $K_{IC}^{app}$.

**Comment 1.25** — Are figure 17 and 19 necessary? They are only briefly discussed and can easily be obtained from geometry considerations.

**Reply:** Thank you for the comment. We believe that both figures address two key questions: the difference between the two models, and the influence of the distance between fibers. For this reason, we think it is better to keep them in the paper.

**Comment 1.26** — Figure 20: The variation seems discontinuous when introducing an overlap between the fibers. How to explain this discontinuity? Is it only an effect of the overlapping or a combined effect of $h_l$ and $h_t$? could the author show the fiber configurations for $h_l < 0$, $h_l = 0$ and $h_l > 0$, for instance as insets in figure 20, as a matter of example?

**Reply:** Thank you for the remark. It is related to both $h_l$ and $h_t$ because when $h_l$ is changed the parameter $h_t$ is also changed to fulfill Equations (13)-(16). A comment has been included, showing the two configurations: Two examples of the geometries obtained for $h_l > 0$ and $h_l > 0$ are shown in Figures 2 and 3. The parameter $h_t$ is obtained to fulfill the geometrical restrictions in Equations (13)-(16).

<figure><img src="image_26.png" /><figcaption>**Figure 2.** Fibre distribution in the inner domain with $h_l > 0$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 3.** Fibre distribution in the inner domain with $h_l < 0$</figcaption></figure>

**Comment 1.27** — Comparison to experiments: While it is discussed in previous sections that the interphase tends to increase $K_{IC}^{app}$, the range of critical SIF given in section 4.4 contains larger values up to 4.5 MPa m$^{0.5}$ without the interphase than with it (up to 3.3 MPa m$^{0.5}$). Is there an inconsistency here?

**Reply:** Thank you for the remark. It was a typo, which is now corrected.
**Comment 1.28** — Page 16: “Note that, unlike the penetration mechanism, the deflection mechanism causes a change in the crack propagation path, which requires additional energy for the path alteration prior to failure.” It is not clear what the authors means by “additional energy for the path alteration prior to failure”.

**Reply:** Thank you for the remark. it is modified as "Note that, unlike the penetration mechanism, the deflection mechanism causes a much longer crack propagation path, which requires additional energy for the increasing path to failure."

Misc.:

**Comment 1.29** — The introduction starts with a focus on gas emission reduction through the use of CMC, is it really pertinent to the present work? Half of the abstract and a part of the introduction is dedicated to this, whereas in the end this work is not related to these aspects so I wonder if it makes sense to display such a general message.

**Reply:** Thank you for the remark. The abstract has been changed in order to focus in the real topic, although a comment is still inside, since it is the main motivation to investigate these materials. The introduction has been changed so that a different perspective can be now perceived.

**Comment 1.30** — The number of digits could be carefully chosen, page 4 for $l_{min}/l_{max}$

**Reply:** The issue about number of digits has been checked and modified if necessary. Thank you for your comment.

**Comment 1.31** — I would avoid using too many terms as “novel”, “innovative”, especially when it refers to a model that is close to previous works by the authors.

**Reply:** Thank you for the remark. This has been changed in the manuscript.

**Comment 1.32** — The proposed approach does not account for the actual composite microstructure – how to improve it and what would be the influence?

**Reply:** Thank you for the remark. The model is a simplified representation of the real microstructure. In a real situation stochastics, defects, different sizes of the fibres can be found. This means that in a real situation several experimental tests must be done to define $K_{IC}$. Our model could be improved considering a mixture with a stochastic tool to allocate the fibres, maybe making a database with all of the results. A comment has been added in the conclusions.

**Comment 1.33** — Page 4: “assuming that .... is assumed”

**Reply:** Thank you for the remark. It has been changed: "...assuming that the influence of the short fibre located at the tip of $\Gamma_c$ is negligible..."

**Comment 1.34** — Page 7: “it may not fully capture the complexities of the real model.” =¿ What is the real model?

**Reply:** Thank you for the remark. It is the real material. This has been changed as “it may not fully capture the complexities of the real material.”

**Comment 1.35** — The number of digits should be carefully chosen (are 2 digits needed when providing the increase in KIC for instance, to be checked through the whole manuscript

**Reply:** Thank you for the remark. This has been homogenized through the whole manuscript. Thank you for the comment.

**Comment 1.36** — Units should not be in italic ($\mu$m, $\mu$ is in italic, the siunitx package could be used using micron)

**Reply:** This has been changed. Thank you for the remark.
Reviewer 2

The authors investigate the resistance to failure of Ceramic matric Composites, of various microstructure. The constitutive phases are (short) fibers embedded in a matrix, which account for an interphase. All constituents are linear elastic and exhibit brittle failure (order of magnitude of $G_c = K_c^2/E$ of 7, 6 or 16 J/m² or the fiber, the interphase and the matrix resp. Tab. 1, while the material strength reported is 1222, 107, 359 MPa).

By considering a preexisting crack, subjected to mode I, of which tip is located at the fiber-matrix interface, the authors somehow extend the seminal problem presented in Evans, He Hutchinson, JACS 1989): here the local stress field are dependent on the microstructure, the fiber diameter and length, and in a second step the influence of the vicinity of other fibers.
The paper is certainly of interest for the JTCAM audience.
After reading, several, some comments and questions can be formulated, that the authors could considered.

We appreciate the reviewer's constructive feedback. In this document, dark blue reflects changes made to the manuscript, while light blue indicates our responses.

**About the problem statement**

**Comment 2.1** — the dimensions of the microstructures are normalized in a smart way. One wonders if the thickness of the interphase e, of which normalization $\hat{e}$ is mentioned, is kept constant for all configurations or is also dependent on the reference length (fiber length or thickness)?

**Reply:** Thank you for the comment. It is constant in all the configurations, since according to the material we are analysing this is not a design parameter, it is kept constant.

**Comment 2.2** — It the initial crack tip located exactly at the fiber/matrix or fiber/interphase? this is worth to mention. Also, additional crack extension originates from this location, thus make the total crack length continuous (as understood from Eqs 2-4). Why not considering a crack tip located to some distance from this interface as in Martin et al. 1998 (quoted in the references)

**Reply:** Thank you for the comment. The pre-existing crack reaches the interphase, as more clearly illustrated in Figures 3 and 4, where the crack deflection occurs at the interface. Introducing a crack located at an arbitrary distance from the fiber would significantly complicate the analysis, as it would require the introduction of an additional parameter when defining the transition to the inner domain. The configuration with a pre-existing crack interacting with the fiber, as shown in Figures 3 and 4, has already been used in Jimenez Alfaro and Leguillon (2022, 2023), and is therefore adopted again in the present study for consistency and comparability.

**Comment 2.3** — In Figs. 2, but also 3 and 4, the 'dark blue' region refers to the homogenized composite. The corresponding effective or homogenized elastic properties are not reported? This is necessary to mention.

**Reply:** Thank you for the remark. The Voigt rule has been used to obtain the properties of the homogenized material. A comment has been included: The material properties of the composite are obtained applying the Voigt rule, and depend on the volume fraction.

**Comment 2.4** — About the inner problem, which is similar to a small scale approach with a remote boundary layer along which dominant K-fields are prescribed, the expression of these remote boundary conditions are worth to be reported

**Reply:** Thank you for the remark. The expression has been indicated in Figure 5. A reference is provided after Equation (10), where the function $u_I(\theta)$, the one of the William's expansion for Mode I, is given.

**About the methodology and outcomes of the study**

**Comment 2.5** — The linear elastic heterogeneous material under consideration here is worth to be analysed with a coupled criterion methodology. Other approaches could be (phase field models for fracture, cohesive models) could (and have been) used. Isn't the force of the couple criterionto provide an concise analysis the onset of (finite) crack initiation? Here the dependence of the competition of two
mechanism with the features of the microstructures (could also be the elastic contrast)? That the analysis
is restricted to the prediction of the onset of crack initiation can also be reported.

**Reply:** Thank you for the remark. A comment has been added in the Section 4.4: It is important to highlight that only crack initiation from the pre-existing crack is considered here, not the crack propagation across a distance larger to crack initiation length

**Comment 2.6** — The present study is inspired by the recent Y. Legerstee's PhD manuscript, in which the toughness of the material under consideration here is measured. Because the material is heterogeneous, one invokes an “effective” or “apparent” toughness which is different to that predicted in the CC analysis presented here: this is worth to mention clearly that crack initiation is considered here, not the crack propagation across a distance larger to one ‘motif’ of the microstructure

**Reply:** Thank you for the remark. As in the previous answer, the same comment could address this issue as well. It is important to highlight that only crack initiation from the pre-existing crack is considered here, not the crack propagation across a distance larger to crack initiation length.

**Comment 2.7** — A description of the mesh, its size, again the boundary conditions used is worth to present in complement to Figs. 2, 3 for instance. The element size with respect to the Irwin length could be mentioned.

**Reply:** Thank you for your comment. A description has been included together with Figure 5.

**Comment 2.8** — No initial stress field is considered, which is likely to originate from the processing. Mentioning this is worth, or indicating why this is not considered.

**Reply:** Thank you for the remark. A comment has been added in the conclusions: Finally, no initial stress field is considered, which is likely to originate from the processing, e.g., residual stresses. This could be included in the model as a future work. An example is shown in Jimenez Alfaro and Leguillon 2023.

**Comment 2.9** — The authors present a 2D plane strain description of a complicated 3D microstructure (as reported in Y Legerstee's manuscript). The present study is likely to provide routes to improve the resistance to failure. Isn't a comparison to measures on a 3D heterogeneous composite reported by Legerstee to consider with caution? The results presented by the authors are interesting, as commented in the sequel.

**Reply:** Thank you for your comment. We fully agree that comparing our 2D plane strain model to 3D experimental results should be done with caution, as the simplifications inherent to 2D modeling cannot capture the full complexity of a real 3D heterogeneous microstructure. However, the aim of the present study is not to provide exact quantitative predictions, but rather to explore and understand the mechanisms governing crack initiation. We believe that, despite the simplifications, the 2D approach offers valuable qualitative insights and can help identify trends that may guide future improvements in failure resistance. This is in line with the motivation behind using simplified models as a first step before more detailed 3D simulations.

## About the results reported

**Comment 2.10** — The authors consider first the interaction of the preexisting crack with a single fiber: variable thickness (height) for a given length. Fig. 5 present a typical and interesting result on the competition between penetration solely versus deflection along the interface solely. Crack through the fiber is predicted as $t < 16 \ \mu m$ approximately, and deflection for a larger value or t, the fiber length being 152 $\mu m$. In Fig. 6, the above fiber diameter $t = 16.6 \ \mu m$ is considered and the fiber length varies. Fig 5 predict a comparable apparent critical toughness which is not found in Fig. 6, the reason for this is not clear or is worth to comment. Fig. 7 with $t = 25 \ \mu m$ is consistent with Fig. 5, however.

**Reply:** Thank you for the remark. We think that the problem is the scale of both graphics. In Figure 7 (in the new manuscript) the point associated with $l = 152 \ \mu m$ is around 1.25 for penetration and 1.28 for deflection. In Figure 6 these two points are the same, but they are closer because the scale of the graphic is bigger.

**Comment 2.11** — Figs. 6-7 suggest that a ratio thickness/length closer to one, results in a larger load level for crack initiation, so called a square fibers. This observation is not related to the cuboidal $\gamma'$ precipitates in Ni-based superalloys which are promoted to improve the hardening and resistance to creep,
no specifically to improve fracture. The comment is a bit misleading, one wonders if this is absolutely
necessary to mention.

**Reply:** Thank you for the remark. The comment was removed in order not to mislead. Consequently, the same comment in Conclusion section was also removed.

**Comment 2.12** — It is also interesting to indicate in Figs. 9-10 the extent of the crack initiation, and its ratio with the fiber edge. How further crack propagation would proceed? In the case of the interfacial crack, mode mixity is likely to be observed along the interface: this seems to be neglected, is there any reason or motivation for this? Would this influence noticeably the competition between penetration vs deflection?

**Reply:** Thank you for the remark. In this work crack nucleation is determined using a technique that avoids explicitly including mixed-mode fracture criteria. This is achieved by relating the solution near the crack tip in heterogeneous domains (highlighted in Figures 3 and 4) to the solution near the crack tip assuming no fibers are present — what we refer to as the outer solution. The latter corresponds to a Mode I crack.

**Comment 2.13** — In the results incorporating multiple fibers (from section 4.2), the transition between the “homogenized composite” region and the detailed one is likely to results in “transient” fields, of the order of one row. Thus, having a larger “window” with more fibers could reduce this effect and one could me more confident on the estimation of the local stress/strain fields, which govern the penetration/deflection competition: did the author investigate this? And how is this likely to influence their predictions.

**Reply:** Thank you for this interesting comment. We did not investigate this specific aspect in the present work. However, we would expect its influence to be similar to that observed when comparing multiple fibres to a single fibre configuration. As shown in Figure 17 of the revised manuscript, this difference is relatively small, which suggests that the influence of the proposed idea might also be limited.

**Comment 2.14** — The influence of the interphase is investigated is interesting and the “softer” interface considered here promotes deflection. This mechanism becomes preferential beyond a thickness larger to approximately 16microns. An analysis of this prediction could be in place here by comparing the surface energy dissipated. Intermediate interphase properties cold be used for illustration, trying to get insight on which parameter, Young modulus of interphase toughness, is likely to govern the competition between penetration and deflection.

**Reply:** Thank you for the remark. We have included a comment in the conclusions: Notice that this computational tool could be used to determine the properties of the interphase that govern the competition between penetration and deflection.

**Comment 2.15** — As reported by the author at the end of section 4.2, caution is in place when comparing these predictions to Lagerstee’s measurements.

**Reply:** Thank you for the remark. We agree with you and we mentioned explicitly and repeatedly the possible limits of our analysis.

**Comment 2.16** — This influence of the fibers relative position is illustrated with the distances $h_t$ and $h_l$. Isn't it relevant to get an estimation of the “magnitude od the initial stress singularity” like the one derived from the Irwin length? Calculating $l_{Irwin} Kc/\sigma c$ result in order of magnitude (well) smaller to the distances considered in the study: isn't this ratio between the dimensions considered here and the Irwin length responsible for the small differences between the multiple vs single fibers descriptions? A comment in addition to the observation in Fig. 18 could be in place here.

**Reply:** Thank you for your comment. The stress singularity here is $\lambda = 0.5$, since it is the one related to the outer solution, that considers an homogeneous material with a crack under Mode I loading conditions. Therefore, it might not have an influence in the difference between the multiple fibre and the single fibre model.

## About the comparison with Lagerstee's measurements

**Comment 2.17** — Discrepancies between a 2D plane strain and measurements on a 3D real materials, and related stochastic effects (fiber orientation, materials properties, defects) are likely to be observed.
One expects to capture qualitatively trends. Again, by mentioning from the beginning that the present work focuses on the conditions for crack initiation rather than the prediction of an effective toughness is worth indicating. The and 3D effects, stochastic aspects of fibers orientations, variability of the parameters are likely to originate a discrepancy quantitatively. The interest and force of the present study is to provide insight on the expected mechanics from ‘simplified’ configurations.

**Reply:** Thank you for the remark. A comment has been included in Section 4.4: Moreover, 3D effects and stochastic features (fiber orientation, materials properties, defects) may lead to quantitative discrepancies.

**Comment 2.18** — Last, the abstract could focus a bit more on the problem position addressed and main outcome, and reduce a little bit the social/industrial motivations.

**Reply:** Thank you for the remark. Some changes have been implemented:

Achieving greater engine efficiency by increasing gas temperatures is a key challenge set by the European Commission for the aerospace industry by 2050. Ceramic matrix composites (CMCs) have emerged as a promising alternative, particularly short-fiber CMCs, which offer enhanced performance in components with complex geometries. This paper presents a novel computational framework to predict the fracture toughness of short-fiber CMCs based on the Coupled Criterion together with the Matched Asymptotics Approach. The influence of various composite parameters has been systematically analysed and validated against existing experimental data. In particular, the influence of fibre diameter, and length, as well as the fibre distance and the volume fraction have been studied. The model has demonstrated reliability and potential as a robust design tool for next-generation aerospace applications.

In conclusion, the authors have carried out an interesting work, of interest for the JTCAM audience and for sure, within the scope of the journal. I have raised several comments which illustrate the interest for the paper and which are also thought to be considered by the authors to clarify (and improve somehow) their paper. I would recommend for publication after these issues carefully considered by the authors.

**Reply:** Thank you very much for the all the interesting comments.