# Determination of interfacial fracture energy of an environmental barrier coating on ceramic matrix composite substrate

Pierre Bertrand, Cédric Huchette, Thibaut Archer, Thomas Vandellos, François Hild

► To cite this version:

Pierre Bertrand, Cédric Huchette, Thibaut Archer, Thomas Vandellos, François Hild. Determination of interfacial fracture energy of an environmental barrier coating on ceramic matrix composite substrate. Journal of the European Ceramic Society, 2025, 45 (2), pp.116886. 10.1016/j.jeurceramsoc.2024.116886. hal-04760047

HAL Id: hal-04760047

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

Submitted on 30 Oct 2024

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
Original article

# Determination of interfacial fracture energy of an environmental barrier coating on ceramic matrix composite substrate

Pierre Bertrand<sup>a,b,c,*</sup>, Cédric Huchette<sup>b</sup>, Thibaut Archer<sup>b</sup>, Thomas Vandellos<sup>a</sup>, François Hild<sup>c</sup>

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

## ARTICLE INFO

### Keywords:

Ceramic Matrix Composite (CMC)

Environmental barrier coating

Interface fracture energy

Digital image correlation

Crack length

## ABSTRACT

The present work concerns the study of the interface fracture energy between a SiC/SiC Ceramic Matrix Composite and an environmental barrier coating. Four-point flexural tests with no precrack were conducted. These tests enable for the stable propagation of two interfacial cracks. They were carried out at room temperature and were instrumented with visible light cameras. This instrumentation allowed for the analysis of the tests thanks to digital image correlation as well as comparisons between experimental and numerical results to locate crack tips and to calculate the interface fracture energy using numerical methods based on linear elastic fracture mechanics. The limits of the method as well as the uncertainties associated with the crack length and the fracture energy were assessed.

## 1. Introduction

Ceramic matrix composites (CMCs) are currently considered as an attractive option to design hot section components of next generation aircraft engines [1]. Their thermomechanical and physical properties at high temperatures and their low density provide increased performance under such extreme environments. In particular, silicon carbide (SiC) fiber-reinforced SiC matrices have shown good stability up to 1300 °C [2]. However, in the presence of water vapor, oxidation and rapid volatilization of the SiC surface occur [3]. Furthermore, the higher the temperature, the faster the kinetics of degradation [4]. To mitigate oxidation effects and extend lifetime [5], an Environmental Barrier Coating (EBC) is required. One of the main concerns regarding CMC/EBC systems is the adhesion between both constituents. The loss of adhesion may lead to spallation of the EBC from the substrate. A good adhesion is therefore essential to extend the lifetime of the system.

To characterize the adhesion of bi-material interfaces, several types of tests exist including tension [6,7], compression [8,9], shear loading [10,11], indentation [12] and flexural [13–18] tests. Among these tests, four-point flexure appears adapted to characterize the adhesion of stiff bi-materials. It was used, for instance, to quantify the adhesion between an aluminum alloy and a polymer layer [13], metallic alloys and ceramic thermal barrier coating (TBC) [14–16] but also on CMC/EBC systems [17]. However, because of the small thickness of the coating [14–18], an additional plate stiffening layer on the coated surface is needed. This additional stiffening part increases the amount

of elastic energy stored in the system and prevents multiple transverse cracking of the TBC or EBC from occurring. The stiffening part must be precracked in order to ensure delamination growth at the interface between the substrate and the coating. The control of the notch depth to induce crack initiation at the interface is challenging. Moreover, delamination growth may be unstable leading to an underestimation of the fracture energy [14,18].

To avoid these limitations, a simplified flexural test is presented hereafter. The setup is based on a four-point flexural test with no stiffening counter plate and no machined notch. To ensure that delamination is located at the interface, the EBC layer is thicker than usual [14,18]. This new geometry was used in a three-point flexural configuration [19] where Digital Image Correlation (DIC) was utilized to estimate the interface crack length, which is key to compute the interfacial fracture energy with Irwin–Kies compliance formula [20]. In the same idea, four-point flexural tests will be exploited using full-field measurements via stereocorrelation. Stereocorrelation is an extension of DIC that registers images acquired by multiple (in this case two) cameras to analyze a common region of interest [21]. This technique allows in-plane and out-of-plane displacement fields to be measured. Full-field measurements have proven their usefulness for the identification of relevant parameters in solid mechanics [22]. A method is proposed herein to identify the crack length based on a comparison between experimental displacement fields and computed ones using

<footnote><sup>a</sup> Safran Ceramics, a technology platform of Safran Tech, Le Haillan, France</footnote>

<footnote><sup>b</sup> MAS, ONERA, Université Paris-Saclay, Chatillon, France</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS–Laboratoire de Mécanique Paris-Saclay, Gif sur Yvette, France</footnote>

<footnote>* Correspondence to: ONERA, Chatillon, France.
E-mail address: pierre.bertrand@onera.fr (P. Bertrand).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Studied multilayered system.</figcaption></figure>

<table><caption>**Table 1**

Elastic properties of both constituents.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>CMC</td><td></td></tr><tr><td>E<sub>11</sub> E<sub>22</sub></td><td>275 GPa</td></tr><tr><td>E<sub>33</sub></td><td>130 GPa</td></tr><tr><td>G<sub>12</sub></td><td>96 GPa</td></tr><tr><td>G<sub>13</sub> G<sub>23</sub></td><td>87 GPa</td></tr><tr><td>Poisson's ratio v</td><td>0.2</td></tr><tr><td>YDS</td><td></td></tr><tr><td>E</td><td>112 GPa</td></tr><tr><td>Poisson's ratio v</td><td>0.22</td></tr></tbody></table>

finite element models. Such identification procedure has already been used to follow crack propagation and to estimate fracture mechanics parameters on refractory materials [23], titanium alloy [24], delamination on thermoset composites [25] and silicon carbides [26]. This method is performed on a local model (i.e. with dimensions corresponding to a zone of interest around the crack) to identify the crack lengths during the test and then validated on a global model (i.e. including the full specimen geometry) by comparing reaction forces between the experiment and the simulation. Since closed-form solutions cannot be utilized due to non-symmetric propagation of two cracks [27,28], a numerical technique based on a contour integral is used to calculate the interface fracture energy from the previously identified crack lengths.

In the following, after the description of the studied multilayer system, the instrumented four-point flexural test adapted to the characterization of the CMC/EBC interface fracture energy is presented. Stereocorrelation is used to measure the displacement fields and to follow crack propagation at the interface on a local Finite Element (FE) model. The method to identify the crack length is then validated on a global model. Last, the method to compute the interfacial fracture energy is described and the results are discussed.

## 2. Materials and experimental setup

### 2.1. Materials and specimen geometry

The tested specimen is a parallelepipedic (100 mm × 10 mm × 4.7 mm) multilayer system. The substrate is a 2.7 mm thick SiC/SiC composite. Then, a thin layer of bond coat, made of silicon (Si), and a thick layer (2 mm) of top coat, made of rare earth yttrium disilicate (Y₂Si₂O₇, designated as YDS in the following), are added by thermal spray on only one surface of the CMC substrate. The aim of the bond coat is to improve the chemical and mechanical compatibility between the ceramic coating and the substrate. During fabrication, a very thin layer of silica is formed by oxidation of silicon. The system is then heated in a furnace to stabilize the coating. This multilayer system is illustrated in Fig. 1.

As the thickness of the bondcoat is very small compared to the top coat (40×), it is assumed that its effect on the macroscopic behavior is not significant. Hence, the bond coat is not considered in the following. Table 1 gathers the elastic mechanical properties of the substrate [29] and of the top coat layer.

### 2.2. Experimental setup and performed tests

The four-point flexural tests were performed on a servohydraulic Instrument testing machine with a constant crosshead speed of 0.25 mm/min. The lateral faces of the sample were coarsely polished to eliminate the excess of EBC on the CMC surface. On one lateral side of the sample, a random speckle pattern was deposited using black paint on the YDS top coat and white paint on the CMC substrate. The tests were instrumented with two visible light cameras at an angle of 25° to acquire images at a frequency of 0.5 Hz. Two tridimensional calibration targets were positioned on the specimen top surface, in the camera field of view, to calibrate the stereosystem using Eikotwin DIC software [30]. The testing machine and an example of images acquired during the test are shown in Fig. 2.

The force/displacement curve is shown in Fig. 3. Four different stages occurred. During the first stage, after a nonlinear response due to first contact of the rollers onto the specimen surface, a linear response is observed during loading, no damage has occurred. During the second stage, a force plateau appears, corresponding to crack initiation at the surface of the YDS layer and its propagation through the top coat thickness. This propagation was unstable. Once the interface was reached, the crack bifurcated and two cracks propagated independently. During the third stage, both cracks propagated in a stable manner along the interface until the failure of the CMC substrate occurred (fourth stage). The stage of interest to characterize the interfacial fracture energy is the third one. The fracture energy of the EBC is of secondary interest in this paper. Before the third stage, the goal was to precrack the specimen to initiate cracks at the interface. Usually, this step is carried out in two separate procedures, first the specimen is precracked and then it is reloaded to have re-initiation of the already present cracks and their stable propagation [13,14,18]. In the present test, the procedure was conducted in one single step. The load drop observed between the first and second stages is representative of transverse EBC cracking.

Post-mortem observations using Scanning Electron Microscopy (SEM) on the lateral sides of the polished specimen showed the bifurcation of the crack at the interface (Fig. 4(a)). The gray level residual field obtained via stereocorrelation (Fig. 4(b)) confirms the location of the cracks at the interface during the test [31]. It is also observed that the propagation of the two interfacial cracks was not symmetric (Fig. 4(b)). Thus, analytical formulas to evaluate the fracture energy based on a critical load are not used as the corresponding hypotheses are not fulfilled [27]. It is worth noting that matrix damage was observed at the end of the third stage. This point will be addressed in Section 3.3.

The fractured surface of the specimen and spallation of the EBC are displayed in Fig. 5. It is observed that the crack propagated between the DSY top coat and the Si bond coat. The crack path is rather tortuous along the specimen width. No preferential interface in the multilayered system is identified as the weakest entity where the cracks propagated. This tortuous path as well as microbranching of the interfacial cracks are shown in the SEM picture of Fig. 4(a).
<figure><img src="image_4.png" /><figcaption>Fig. 2. Flexural setup with its instrumentation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Force vs. displacement curve of one four-point flexural test.</figcaption></figure>

## 2.3. Stereocorrelation

In this analysis, global FE-based stereocorrelation [32] was used to compare experimental and simulated kinematic fields. The first step of the procedure to extract kinematic fields is the calibration of the two cameras. It is based on the pinhole camera model [33], which allows the 3D FE nodes to be projected onto 2D image planes of each camera c. This step is achieved by determining the projection matrices $[M^c]$ [21].

First, a pre-calibration step was performed to initialize projection matrices using a Point and Perspective (PnP) method [34]. Ten points were selected for each camera. This step was simplified by the presence of tridimensional calibration targets mentioned earlier, which gave access to characteristic points in different planes. Second, a fine tuning of the projection matrices was performed by identifying the camera parameters that minimized reprojection errors.

Once the calibration of both cameras was performed, global stereocorrelation was run to measure displacement fields over the whole lateral surface of the specimen [35]. This approach facilitates the transfer of measured displacements from experiments to simulations.

# 3. Description of the procedure and results

## 3.1. Methods to follow crack propagation via full-field measurements

Full-field measurements were first used to locate the crack path on the FE model, second, to isolate the left interface crack from the right one with two local models, third, to identify the crack length on the local models. As shown in Fig. 4(b), the presence of cracks induced high gradients in the kinematic fields. Gray level residuals were used to locate the transverse cracks as shown in Fig. 4(b).

It was assumed that the crack, which initiated in the top coat, was straight and perpendicular to the interface. Once the crack in the EBC is located, two local FE models are built, one for each interfacial crack from the transverse crack (i.e. one where the crack propagated to the right and the other to the left as shown in Fig. 6).

Each model represents a part of the sample including its full thickness, starting at the transverse crack location and with a long enough length $L_0$ to contain the left (or right) interfacial crack length $L$ ($L < L_0$) as shown in Fig. 7. The FE analyses were performed using Abaqus Standard [36]. For the first model, an initial crack with length $L$ was introduced in the mesh by node splitting. Boundary conditions corresponding to the experimental displacement fields measured by stereocorrelation were applied on the exterior edges of the local model except for the free edges where no load was applied (Fig. 7). The cost function to minimize

$$ \chi^2(L) = \frac{1}{2N} \sum_{n=1}^{N} \left( \frac{(U_X^c(L) - U_X^e)^2 + (U_Y^c(L) - U_Y^e)^2}{\sigma_u^2} \right) \quad (1) $$

is written in terms of the difference between computed nodal displacements ($U_X^c$ and $U_Y^c$), with X and Y respectively the horizontal and vertical directions, and experimental nodal displacements ($U_X^e$ and $U_Y^e$) measured via stereocorrelation. N is the number of nodes used to evaluate the cost function and $\sigma_u$ the standard uncertainty associated with experimental displacements. Then, the crack length L is increased by an incremental length and the cost function is computed again.

The area close to the crack was not considered in the cost function not to take into account perturbations due to discontinuous fields. The crack length was increased by an increment of 100 µm until the minimum of the cost function was reached.

This method allows both cracks to be separated and to study them one by one. The mesh used in the simulations was refined in comparison to the measurement discretization to achieve a finer resolution for
<figure><img src="image_6.png" /><figcaption>Fig. 4. Illustration of crack bifurcation at the interface before final failure of the specimen via (a) SEM and (b) gray level residual field for a 750 N applied force.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Fractured surface at the interface on the CMC and EBC sides.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Local FE models to identify the crack lengths.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Local model to identify the crack length (example for the left crack).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. (a) Change of the cost function $\chi^2$ as a function of the considered left crack length $L$. (b) Crack lengths identified using full-field measurements.</figcaption></figure>

the detection of the crack tip. The mesh size used for the simulations was 0.1 mm. It is thus the minimal crack length variation over which the cost function was evaluated. To get a better resolution on the identified crack length, the cost function $\chi^2$ was interpolated with a parabolic polynomial (Fig. 8(a)).

The identified crack length corresponds to the minimum of the parabolic interpolation as shown in Fig. 8(a). The evolution of both crack lengths is shown in Fig. 8(b). The procedure is then repeated using nodal displacements measured via stereocorrelation on the next image. The same loop is performed on the other local model to identify the crack length on the other side of the transverse crack. With this procedure, the information on crack propagation is only provided on the surface. It is thus mandatory to make some assumptions when considering the width of the specimen. The shape and the length of the cracks were assumed to be the same through the width of the specimen. This strong assumption is discussed in the next section.

## 3.2. Limit of the method to identify the crack length

The method developed to identify the crack length has been performed on two four-point flexural tests. As the identification procedure is based on experiment/simulation comparisons, the uncertainty level of the crack length is directly linked to those associated with the experimental displacement measurements. The standard uncertainty of the nodal displacement from the stereocorrelation measurement is 0.3 µm. This value was obtained by performing stereocorrelations over 30 reference images (i.e. with no applied load). To assess the uncertainty of the sought crack length $L$, the kinematic sensitivities [37] of the measured data with respect to $L$ are computed

$$S_U(x, L) = \frac{\partial U^c(x, L)}{\partial L} \qquad (2)$$

where $U^c(x, L)$ is the computed nodal displacement. In practice, the sensitivities are computed via forward finite differences between a first simulation where the crack length is the one identified with the method previously described and a second simulation where the crack length is increased by 1% of the identified one. The sensitivities are gathered in a column vector $\{\frac{\partial U^c}{\partial L}\}$. It is now possible to compute the associated scalar Hessian $H_{LL}$

$$H_{LL} = \left\{ \frac{\partial U^c}{\partial L} \right\}^T \left\{ \frac{\partial U^c}{\partial L} \right\} \qquad (3)$$

This quantity is used to determine the standard uncertainty of the crack length denoted $\sigma_L$. In this case, the input data are considered unbiased, Gaussian and white [32]. Therefore, the variance of the identified parameter $\sigma_L^2$ is proportional to the inverse Hessian [38]

$$\sigma_L = \frac{\sigma_U}{\sqrt{H}} \qquad (4)$$

In the present, a value of 1 µm was found.

In Section 3.1, it was mentioned that measured boundary conditions were prescribed on the exterior edges. As the measured displacement has its own uncertainty, its contribution must be taken into account on the final value of the identified crack length $L$. The propagation of uncertainties was conducted by computing sensitivity vectors as previously explained. In practice, for $i$ nodes where experimental boundary conditions are prescribed, $2i$ simulations are run. For each simulation one nodal displacement prescribed as boundary condition was perturbed by 1% of its initial values while the others remained equal to their initial value. The vertical and horizontal components are perturbed independently so two simulations are needed per node. The sensitivity vector (Eq. (2)) becomes a rectangular matrix where each column is the sensitivity to one degree of freedom. The last column is the sensitivity associated with the crack length previously computed. As $\{\eta_U\}$ and $\{\eta_p\}$ respectively the uncertainties on the displacement fields and on the parameters $\{p\}$ (i.e boundary conditions and crack length) are related by

$$\{\eta_p\} = [S]^T [S]^{-1} [S]^T \{\eta_U\} \qquad (5)$$

To evaluate the uncertainty level on the crack length, the full covariance matrix associated with the parameters is needed

$$\langle \{\eta_p\} \{\eta_p\}^T \rangle = \sigma_U^2 [H] \qquad (6)$$

because $\langle \{\eta_U\} \{\eta_U\}^T \rangle = \sigma_U^2 [I]$ and with

$$[H] = [S]^T [S] \qquad (7)$$

where the diagonal terms are the autosensitivities where the correlations are not accounted for. The total variance of the crack length is [37]

$$\sigma_L^2 = \sigma_U^2 \sum_i \frac{V_{kl}^2}{D_{il}} \qquad (8)$$

where the matrix $[V]$ gathers all column eigenvectors, and $[D]$ is the diagonal matrix of eigenvalues of the Hessian matrix, $k$ corresponds to the last column (or line) of the Hessian matrix, $L$ being the last parameter considered to assemble the sensitivity matrix. The standard uncertainty of the crack length was equal to 30 µm with this approach. It is significantly higher than the 1 µm level reported without taking into account the boundary condition contributions on the crack length uncertainty. Yet, this uncertainty level is lower than what has been achieved with crack tip detection algorithms based on direct image processing [39] or convolutional neural networks [40]. Moreover, DIC has the advantage of acquiring kinematic information over a larger region of interest. Some studies proposed to identify the crack length based on the discontinuity in the kinematic fields computed via DIC and
<figure><img src="image_11.png" /><figcaption>Fig. 9. Global model with cracks and contour around the left crack tip used to evaluate J-integrals.</figcaption></figure>

the uncertainty levels reported were significantly higher (270 µm) [41]. Moreover, there is a clear benefit of using this approach because it is mechanically based.

### 3.3. Validation of the method on a global model

The crack lengths previously identified are now used on a global model to compare the measured force by the load cell to the computed reaction forces. The global model includes cracking through the thickness of the EBC and along the interface. Boundary conditions from stereocorrelation are applied on four nodes corresponding to the contact points of the rollers onto the specimen (Fig. 9). They are updated in accordance with the corresponding image (i.e. one simulation is run for each image). The simulation is conducted in 2 dimensions with the assumption of plane strain states. It is worth mentioning that a strong hypothesis was made by considering that the crack lengths and the boundary conditions are identical through the width of the specimen. Such assumption has to be made considering that stereocorrelation gives surface information.

First, only the elastic stage is considered, the experimental force is compared with the sum of the vertical component of reaction forces extracted at the 2 nodes corresponding to the lower contact with the rollers. These quantities are extracted as post-processing for each considered image. This pre-analysis allows for the validation of the model hypothesis (material properties, boundary conditions, 2D assumption).

The next analysis consists in introducing the cracks, identified by the procedure described in Section 3.1 in the global model to simulate the degradation of the global stiffness of the system. As for the local models, transverse and interfacial cracks are introduced by splitting the nodes in the mesh. Thus the simulation is still linear elastic. The comparison between experiment and simulation shows a very good correlation (Fig. 10) until image 80. The loss of global stiffness is well described. After image 80 the simulated response starts to deviate from the experiment.

Let us mention that this curve highlights the importance of using directly displacement fields computed via DIC as boundary conditions instead of “Ideal” boundary conditions (i.e. the displacement prescribed by the machine crosshead on the bottom rollers and no displacement of the top rollers) if one wants to minimize the gap between simulated response and experimental one. “ideal” boundary conditions result in a different force/displacement response as illustrated in Fig. 11(b). That discrepancy is due to the fact that “ideal” boundary conditions do not

<figure><img src="image_12.png" /><figcaption>Fig. 10. Comparison between simulated (simu) and experimental (exp) reaction forces.</figcaption></figure>

take into account the displacement of the top rollers and the left/right side asymmetries as displayed in Fig. 11(a).

This comparison allows us to validate the developed strategy to locate and identify the crack length via its effects on the kinematic fields using local FE models. Moreover the comparison between experiments and simulation for the global model (Fig. 10) shows the limits of validity of this method. One may assume that after image 80 new dissipative phenomena appeared in addition to interface crack propagation, which are not taken into account in the global model (Fig. 9). Among these phenomena, it is very likely that damage arose in the CMC substrate during the third stage in Fig. 3 before its final failure. Some study performed on similar materials showed that as the SiC matrix of the substrate is brittle, microcracking occurs at very low tensile load [10]. In the four-point flexural test presented herein, the bottom part of the CMC is subjected to tensile stress, which is more important in the zone of debonding with the EBC. It is expected to have matrix microcracks in this zone, that is why the overall stiffness of the system is higher in simulations than in the experiment. This subject will be discussed in more details in Section 4.1. For the moment, only the first images of the third stage (Fig. 3), herein from images 50 to 80, were considered in the energy release rate computations presented hereafter.
<figure><img src="image_13.png" /><figcaption>Fig. 11. (a) Experimental displacement applied as boundary conditions on the four nodes corresponding to the rollers. (b) Comparison between reaction forces with “ideal” (ideal BC) and experimental (exp BC) boundary conditions with the experimental (exp) force.</figcaption></figure>

### 3.4. Interfacial energy release rate

With the previously identified crack lengths, it is possible to compute the macroscopic fracture energy of the interface, in the sense that the bond coat is neglected as well as the roughness geometry of the interface. In the present case, the fracture energy corresponds to the critical energy release rate [42]. The critical energy release rate $G_C$ was computed using an interaction integral (J-integral), which has been used for bimaterials [43]. The J-integral allows for the computations of energy release rates for the prescribed boundary conditions. In the present case, as the cracks are propagating in a stable manner for the considered images, the J-integral gives access to the fracture energy $G_C$. The advantage of using such numerical technique is that the fracture energy $G_C$ is computed individually around each interfacial crack, namely, the one propagating on the right and on the left. The evaluation of $G_C$ was conducted after a stable crack propagation regime was established on both sides (i.e after image 50 for both tests). It is worth mentioning that such method does not require the force or any unloading to compute the energy release rate. The force only serves to validate the approach as described in Section 3.3. A preliminary study showed that the J-integral must be evaluated at least around 10 contours of elements around the crack tip to get a converged value of $G_C$. The 10th considered contour is illustrated in Fig. 9. The mesh is refined around the crack tip to capture the stress singularity in this zone (Fig. 9). $G_C$ was computed on the global model previously used to validate the method.

## 4. Discussion

### 4.1. Evolution of the fracture energy

The fracture energy has been evaluated with the J-integral for two four-point flexural tests during the stable crack propagation phase at the interface. The change of $G_C$ during the two tests is reported in Fig. 12.

As previously mentioned, only the part until image 80 is considered for the evaluation of $G_C$. For both tests, the simulation response deviated from the experimental response after image 80. There is a tendency of greater values of $G_C$ when the crack length increases. This may be related to the change in mode mixity during crack propagation [44]. In the present case, the mode mixity $\alpha$ is defined as a ratio of the energy release rate in mode II $G_{II}$ with respect to the total energy release rate

$$ \alpha = \frac{G_{II}}{G_I + G_{II}} \qquad (9) $$

where $G_I$ and $G_{II}$ components are evaluated using the Virtual Crack Closure Technique [45]. As for energy release rate evaluations with the J-integral, the VCCT is applied to the left and right cracks individually

so that the mode mixity is evaluated for both cracks. The ratio of the energy release rate in mode II $G_{II}$ with respect to the total energy release rate is higher for the right crack for both tests. It ranges between 0.25 and 0.28 for the left crack, and between 0.29 and 0.40 for the right crack (Fig. 13).

If each crack is considered separately, $\alpha$ is constant over the domain of propagation with the mean values of $\alpha$ being equal to 0.26 for the left crack for both tests and 0.32 (test1) and 0.34 (test2) for the right crack. Thus, the increase observed in Fig. 12 was due to other phenomena. In particular, as mentioned before, microcracking of the substrate may lead to an overestimation of the energy required to propagate the cracks. To analyze this phenomenon at a more local scale, maximum principal strain levels in the CMC were investigated on the global model (Fig. 9). The evolution of maximum principal strain level in the CMC substrate is shown in Fig. 14(a) between images 50 and 80. The zone where the maximum principal stress is the highest is located above the transverse crack in the EBC (Fig. 14(b)). This region was subjected to the highest tensile stress that gives rise to matrix microcracking of the CMC substrate. In the literature, a study performed on a similar SiC/SiC CMC reported that non linearity of the stress/strain curve appeared at a level of 0.001 [29]. This threshold of principal strain is reached after image 60 for test 1 and after image 55 for test 2. Thus, the elastic assumption made in the FE simulation was valid only for the first images. For these images, before reaching the 0.001 threshold, $G_C$ was constant (Fig. 12(a)). After image 60, it starts to increase for test 1 and the same trend is observed for test 2 after image 55. This increase is due to the presence of microdamage in the CMC substrate. This microdamage does not have an impact on the macroscopic force (Fig. 10) until image 80 for test 1.

The mode mixity change does not explain the increase in fracture energy $G_C$ with respect to the crack length. The mode mixity is slightly different between the cracks propagating on the left and on the right. However, it alone cannot explain such a difference in $G_C$ values between the two sides. It is worth mentioning that, as shown in Fig. 5, the crack at the interface is propagating between the Si bond coat and the DSY layer over the width of the specimen. Such an observation was made for the crack propagating to the right. Surface observations, on the lateral side of the specimen, showed a similar tortuous path for the crack propagating on the left. There is no major difference on the crack path location between left and right sides to explain $G_C$ differences. A better understanding of such a gap may require an analysis at a smaller scale, such as the interface roughness. In the literature, the estimation of $G_C$ at the interface of CMC/EBC systems via flexural tests is not very common. Some values have been reported on the same system but after ageing the specimens, which means that the TGO thickness was higher and the interface was degraded, they range from 12 J m⁻² to 20 J m⁻² [18]. Those values were evaluated on
<figure><img src="image_14.png" /><figcaption>Fig. 12. Evolution of $G_c$ for (a) test 1 (b) test 2.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Evolution of $\alpha$ the ratio of the energy release rate in mode II $G_{II}$ with respect to the total energy release rate for (a) test 1 (b) test 2.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. (a) Change of the maximum principal strain level in the CMC substrate for the two tests. (b) Maximum principal strain field around the transverse crack (image 80).</figcaption></figure>

4-point flexural tests during which unstable propagation of the cracks was reported contrary to the present case for which stable propagation occurred. In the present study, a 4-point flexural configuration was chosen. However, in such a test, one has to deal with two interfacial cracks, which did not propagate symmetrically. To overcome this issue, a modified clamped beam geometry in flexure may be an interesting option since a single interfacial crack then propagates [46].

## 4.2. Impact of uncertainty sources on the fracture energy

To further analyze the uncertainty level associated with $G_c$, a sensitivity analysis is performed. Six parameters {$p$} are identified as sources of uncertainties on $G_c$, namely, the four nodal displacements applied as boundary conditions on the four contact points and two crack lengths at the interface (right and left crack). The sensitivity of $G_c$ with respect to these parameters was evaluated by computing
<table><caption>**Table 2**
Sensitivities of $G_C$ associated with each parameter.</caption><thead><tr><th>Parameter</th><th>Sensitivity <sup>ΔG</sup>/<sub>G<sub>C</sub></sub> (%)</th></tr></thead><tbody><tr><td>Top right displacement</td><td>0.11</td></tr><tr><td>Top left displacement</td><td>0.57</td></tr><tr><td>Bottom right displacement</td><td>0.94</td></tr><tr><td>Bottom left displacement</td><td>1.2</td></tr><tr><td>Right crack length</td><td>0.07</td></tr><tr><td>Left crack length</td><td>0.64</td></tr></tbody></table>

the forward finite differences between the calibrated value of $G_C$ (i.e. computed with initial parameters) and the perturbed value of $G_C$ for which one of the 6 parameters is increased by 1% of its initial value (the other 5 parameters remain unperturbed). This procedure was repeated for all 6 parameters. The results are gathered in Table 2.

This analysis shows that the most sensitive parameters are the boundary conditions applied to the bottom nodes. The right crack has a negligible impact on the $G_C$ value because in this example, $G_C$ is computed for the left crack. The total uncertainty on $G_C$ was computed with the differential of $G_C$

$$dG_C = \sum_p \frac{\partial G_C}{\partial p} \delta p \qquad (10)$$

by taking the variance

$$\sigma_{G_C}^2 = \langle dG_C^2 \rangle = \sum_p \left( \frac{\partial G_C}{\partial p} \right)^2 \sigma_p^2 \qquad (11)$$

gives the standard uncertainty $\sigma_{G_C}$ on $G_C$. With the data reported in Table 2, the standard uncertainty on $G_C$ is $0.5 J m^{-2}$. The uncertainty level of $G_C$ evaluated for the right crack and the left crack are similar and varies from $0.2 J m^{-2}$ to $0.7 J m^{-2}$. The uncertainty associated with the crack length does not have a strong impact on the final value.

## 5. Conclusion

A four-point flexural setup was successfully employed to initiate and propagate stable cracks at the interface between a CMC and its EBC. This test does not need the use of a counter plate or a machined precrack before the test. It is thus easier to prepare the specimen and put it in place in the testing machine.

Crack propagation along the interface was followed via FE-based stereocorrelation, which was found very useful to compare displacement fields from experiments and finite element simulations. This comparison allowed the crack lengths to be identified during the test by their effect on the kinematic fields in the region of interest. This procedure was used to assess the crack length along the interface during the test with two local models in order to study both cracks separately. As the crack length identification was performed on a local model, the relevance and accuracy of this method was validated at a global scale taking into account the full specimen. Simulations were performed with inserted cracks in the mesh and the computed reaction forces were compared to the experimental measurements. This step highlighted the importance of applying boundary conditions from experimental displacements in order to reproduce the test kinematics in FE simulations. With the global simulation, a good correlation was observed with the experimental data, which validated the proposed method to identify the crack length up to a certain level where the two results started to deviate. This difference was attributed to an additional dissipative phenomenon, which corresponds to microdamage in the CMC substrate that was not taken into account in the simulations and gave the domain of validity of this method.

As soon as the crack tip position was determined, the critical energy release rate $G_C$ was computed on the numerical model during crack propagation. The value of $G_C$ was consistent with values found in the literature [18]. The value of $G_C$ was influenced by uncertainties associated with the crack tip position and measured boundary conditions.

The proposed method to identify the crack length and to compute $G_C$ permits to take into account these uncertainties and to quantify their impact. Such method can be easily adapted to various loading conditions as long as a local model can be built around the crack without perturbations of the displacement field caused by the proximity with other cracks for instance. This method has been carried out on two flexural tests where only one transverse crack was present. However, in some cases, multiple transverse cracks initiated in the coating. In these cases, the transverse cracks bifurcated at the interface but the interfacial cracks did not propagate over a distance as long as in the cases where only one transverse crack was present. The local zones around each interfacial crack were very small because the cracks were close to each other. Hence the quantity of kinematic information used during the cost function evaluation was reduced.

## CRediT authorship contribution statement

**Pierre Bertrand:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Conceptualization. **Cédric Huchette:** Writing – review & editing, Methodology, Investigation, Conceptualization. **Thibaut Archer:** Writing – review & editing, Methodology, Investigation, Conceptualization. **Thomas Vandellos:** Writing – review & editing, Methodology, Investigation, Conceptualization. **François Hild:** Writing – review & editing, Methodology, Investigation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported under ATLAAS, French research project cofunded by DGAC, ONERA and SAFRAN Group, involving SAFRAN Group, ONERA and CNRS.

## References

[1] E. Bouillon, Ceramic matrix composite behavior enhancement for gas turbines hot sections, in: Conference of NATO Science & Technology Organization, AVT-356 Research Symposium on Physics of Failure for Military Platform Critical Subsystems, 2021.

[2] J. DiCarlo, H. Yun, G. Morscher, R. Bhatt, SiC/SiC composites for 1200°C and above, in: Handbook of Ceramic Composites, 2004.

[3] K.L. More, P.F. Tortorelli, L.R. Walker, N. Miriyala, J.R. Price, M. van Roode, High-temperature stability of SiC-based composites in high-water-vapor-pressure environments, J. Am. Ceram. Soc. 86 (8) (2003) 1272–1281.

[4] E.J. Opila, J.L. Smialek, R.C. Robinson, D.S. Fox, N.S. Jacobson, SiC recession caused by SiO₂ scale volatility under combustion conditions: II, thermodynamics and gaseous-diffusion model, J. Am. Ceram. Soc. 82 (7) (1999) 1826–1834.

[5] K.N. Lee, Current status of environmental barrier coatings for Si-Based ceramics, Surf. Coat. Technol. 133-134 (2000) 1–7.

[6] M.P.K. Turunen, P. Marjamäki, M. Paajanen, J. Lahtinen, J.K. Kivilahti, Pull-off test in the assessment of adhesion at printed wiring board metallisation/epoxy interface, Microelectron. Reliab. 44 (6) (2004) 993–1007.

[7] T. Vandellos, M. Hautier, N. Carrere, C. Huchette, Development of a new fracture test to identify the critical energy release rate: The tensile flexure test on notched specimen, Eng. Fract. Mech. 96 (2012) 641–655.

[8] W. Zhu, L. Yang, J.W. Guo, Y.C. Zhou, C. Lu, Determination of interfacial adhesion energies of thermal barrier coatings by compression test combined with a cohesive zone finite element model, Int. J. Plast. 64 (2015) 76–87.

[9] S.V. Nair, H.E. Eaton, E.Y. Sun, Measurements of interface strength and toughness in shear of environmental barrier coatings on ceramic substrates at ambient and at elevated temperature, Surf. Coat. Technol. 200 (18) (2006) 5175–5180.

[10] Y. Aoki, J. Inoue, Y. Kagawa, K. Igashira, A simple method for measurement of shear delamination toughness in environmental barrier coatings, Surf. Coat. Technol. 321 (2017) 213–218.

[11] E. Kawai, H. Kakisawa, A. Kubo, N. Yamaguchi, T. Yokoi, T. Akatsu, S. Kitaoka, Y. Umeno, Crack initiation criteria in EBC under thermal stress, Coatings 9 (11) (2019) 697.
[12] J. Yan, T. Leist, M. Bartsch, A.M. Karlsson, On cracks and delaminations of thermal barrier coatings due to indentation testing: Experimental investigations, Acta Mater. 56 (15) (2008) 4080-4090.

[13] P.G. Charalambides, J. Lund, A.G. Evans, R.M. McMeeking, A test specimen for determining the fracture resistance of bimaterial interfaces, J. Appl. Mech. 56 (1) (1989) 77-82.

[14] I. Hofinger, M. Oechsner, H.-A. Bahr, M.V. Swain, Modified four-point bending specimen for determining the interface fracture energy for thin, brittle layers, Int. J. Fract. 92 (1998) 213-220.

[15] Y. Zhao, A. Shinmi, X. Zhao, P.J. Withers, S. Van Boxel, N. Markocsan, P. Nylen, P. Xiao, Investigation of interfacial properties of atmospheric plasma sprayed thermal barrier coatings with four-point bending and computed tomography technique, Surf. Coat. Technol. 206 (23) (2012) 4922-4929.

[16] J.-R. Vaunois, M. Poulain, P. Kanouté, J.-L. Chaboche, Development of bending tests for near shear mode interfacial toughness measurement of EB-PVD thermal barrier coatings, Eng. Fract. Mech. 171 (2017) 110-134.

[17] B.T. Richards, D. Zhu, L.J. Ghosn, H.N. Wadley, Mechanical properties of air plasma sprayed environmental barrier coating (EBC) systems: Preliminary assessments, in: Developments in Strategic Ceramic Materials, John Wiley & Sons, Ltd, 2015, pp. 219-237.

[18] M. Fernandez, G. Couégnat, F. Rebillat, Relation entre oxydation et évolution de l'adhérence d'une barrière environnementale (in French), in: Comptes Rendus des JNC 21 – Bordeaux INP, 2019.

[19] Q. Zhu, W. He, J. Zhu, Y. Zhou, L. Chen, Investigation on interfacial fracture toughness of plasma-sprayed TBCs using a three-point bending method, Surf. Coat. Technol. 353 (2018) 75-83.

[20] G. Irwin, J. Kies, Critical energy rate analysis for fracture strength, Weld. J. Res. Suppl. 19 (1954) 193-195.

[21] M.A. Sutton, J.J. Orteu, H. Schreier, Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications, Springer Science & Business Media, 2009, Google-Books-ID: AlkqMxpQMLsC.

[22] M. Grédiac, F. Hild, Full-Field Measurements and Identification in Solid Mechanics, John Wiley & Sons, 2012, google-Books-ID:JUawSvvgvkcC.

[23] R. Vargas, J. Neggers, R.B. Canto, J.A. Rodrigues, F. Hild, Comparison of two full-field identification methods for the wedge splitting test on a refractory, J. Eur. Ceram. Soc. 38 (16) (2018) 5569-5579.

[24] F. Mathieu, F. Hild, S. Roux, Identification of a crack propagation law by digital image correlation, Int. J. Fatigue 36 (1) (2012) 146-154.

[25] F. Mathieu, P. Aimedieu, J.-M. Guimard, F. Hild, Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations, Composites A 49 (2013) 203-213.

[26] P. Forquin, L. Rota, Y. Charles, F. Hild, A method to determine the macroscopic toughness scatter of brittle materials, Int. J. Fract. 125 (1) (2004) 171-187.

[27] P.G. Charalambides, H.C. Cao, J. Lund, A.G. Evans, Development of a test method for measuring the mixed mode fracture resistance of bimaterial interfaces, Mech. Mater. 8 (4) (1990) 269-283.

[28] J. Hutchinson, Z. Suo, Mixed mode cracking in layered materials, Adv. Appl. Mech. 29 (1992) 193-195.

[29] Y. Gowayed, G. Ojard, R. Miller, U. Santhosh, J. Ahmad, R. John, Correlation of elastic properties of melt infiltrated SiC/SiC composites to in situ properties of constituent phases, Compos. Sci. Technol. 70 (3) (2010) 435-441.

[30] Eikosim, Etalonnage d'une caméra : Principes Et procédures, EikoSim, 2021, https://eikosim.com/articles-techniques/etalonnage-camera-principes-et-procedures.

[31] Y. Chi, W. Zhao, B. Pan, Gray level residual field: an effective metric for pixelwise matching quality evaluation in local digital image correlation, Meas. Sci. Technol. 34 (7) (2023) 075018, publisher: IOP Publishing.

[32] F. Hild, S. Roux, Comparison of local and global approaches to digital image correlation, Exp. Mech. 52 (9) (2012) 1503-1519.

[33] O. Faugeras, Three-Dimensional Computer Vision: A Geometric Viewpoint, MIT Press, Cambridge, Mass, 1993, open Library ID:OL1401278M.

[34] V. Lepetit, F. Moreno-Noguer, P. Fua, EPnP: An accurate O(n) solution to the PnP problem, Int. J. Comput. Vis. 81 (2) (2009) 155-166.

[35] F. Hild, S. Roux, Digital image correlation: from displacement measurement to identification of elastic properties – a review, Strain 42 (2) (2006) 69-80.

[36] D.S. Help 3DS, Simulia user assistance 2022, 2022, https://www.3ds.com/support/documentation/user-guides.

[37] R. Gras, H. Leclerc, F. Hild, S. Roux, J. Schneider, Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization, Int. J. Solids Struct. 55 (2015) 2-16.

[38] S. Roux, F. Hild, 2.02 – Comprehensive Full-Field Measurements via Digital Image Correlation, in: Vadim Silberschmidt (Ed.), Comprehensive Mechanics of Materials, Elsevier, Oxford (UK), 2024, pp. 3-56.

[39] H.L. Gao, S.S. Shen, Y. Yun, Fatigue crack length real time measurement method based on camera automatically tracking and positioning, Appl. Mech. Mater. 130-134 (2012) 3111-3118.

[40] C. Brooks, N. Rajic, Automated visual tracking of crack growth in coupon and component level fatigue testing using thermoelastic stress analysis, Int. J. Fatigue 163 (2022) 107037.

[41] B.V. Farahani, F. Direito, P.J. Sousa, P.J. Tavares, V. Infante, P.P.M.G. Moreira, Crack tip monitoring by multiscale optical experimental techniques, Int. J. Fatigue 155 (2022) 106610.

[42] A. Griffith, The Phenomena of Rupture and Flow in Solids, Royal Society of London, 1921.

[43] C.F. Shih, Cracks on bimaterial interfaces: elasticity and plasticity aspects, Mater. Sci. Eng. A 143 (1) (1991) 77-90.

[44] A.G. Evans, M. Rühle, B.J. Dalgleish, P.G. Charalambides, The fracture energy of bimaterial interfaces, Metall. Mater. Trans. 21 (9) (1990) 2419-2429.

[45] R. Krueger, Virtual crack closure technique: History approach, and applications, Appl. Mech. Rev. 57 (2) (2004) 109-143.

[46] A.K. Mishra, S. Abbas, D. Srinivasan, S. Sampath, B.N. Jaya, Novel, high throughput interface fracture testing method for thermal spray coatings: The modified cantilever bend technique, J. Eur. Ceram. Soc. 44 (12) (2024) 7285-7295.