Image-Based Numerical Modeling of Self-Healing in a Ceramic-Matrix Minicomposite

Abstract

Self-healing, obtained by the oxidation of a glass-forming phase, is a crucial phenomenon to ensure the lifetime of new-generation refractory ceramic-matrix composites. The dynamics of oxygen diffusion, glass formation and flow are the basic ingredients of a self-healing model that has been developed here in 2D in a transverse crack of a mini-composite. The presented model can work on a realistic image of the material section and is able to simulate the healing process and to quantify the exposure of the material to oxygen: a prerequisite for its lifetime prediction. Crack reopening events are handled satisfactorily, and healing under cyclic loading can be simulated. This paper describes and discusses a typical case in order to show the model capabilities.

Highlights

• The numerical model features oxidation kinetics and liquid oxide flow
• The model works on material images or on virtual material meshes (incorporates realistic fiber arrangements and size distributions)
• It describes self-healing restart after crack reopening

1. Introduction

Self-healing Ceramic-Matrix Composites (SH-CMCs) are very promising materials for use as civil aircraft jet engine hot parts [1]. Tests in true conditions of use have been carried out, showing excellent lifetimes [2]. SH-CMCs are constituted of a 3D woven fabric of SiC fibers, infiltrated by an interphase and a multilayer matrix.The Pyrocarbon (PyC) interphase acts as a crack deviator and prevents premature fiber failure while the matrix undergoes multiple cracking [3]. In addition, SiC matrix layers bring stiffness and are chemically inert at moderate temperatures: boron-containing phases produce a liquid oxide, which fills cracks above 450 ${}^{\circ }$C and prevents further oxidation by a diffusion barrier effect [4]. This increases by large amounts the material lifetime under exposition to high-temperature oxidation [5]. Yet, the massive production of SH-CMC parts in engines can only be foreseen if there is a sufficient confidence in the material lifetime duration. The challenge is that the material is expected to last as long as ten years; hence, a material development cycle based only on experimental characterization is not feasible. Modeling is mandatory to incorporate as much understanding of the material physics as possible and to identify and explain the experiments that have been carried out in “accelerated degradation” conditions or on material parts. The present work has been developed in this frame.

The lifetime of SH-CMCs results from a balance between degradation and self-healing [6]. The network of matrix cracks induced by thermomechanical loadings brings the oxygen in contact with all the inner material constituents; they are transformed gradually into oxides, some of them being gaseous (e.g., carbon), liquid (e.g., boron), or solid (e.g., silica at moderate temperatures). The calculation of the chemical reaction rate everywhere in the material is a key point for lifetime prediction. Past reports have described models that give the time evolution of the oxygen concentration in, say, some key part of the composite (i.e., the weakest, most exposed fiber); they deduce from this quantity a global lifetime because the tensile strength of this fiber can be related to oxygen exposure through an oxygen-dependent slow crack growth rate law [7,8]. These models were developed in zero [9,10] or one [11] space dimension.

In reality, the materials have a complex structure, as illustrated in Figure 1. SiC filament bundles are woven in a 3D interlocking pattern. Then, the multi-layered matrix is deposited, usually by Chemical Vapor Infiltration (CVI), around the fiber, with varying thicknesses depending on the location: large close to the part boundaries and small inside the fiber bundles. Pores and cracks are distributed unevenly throughout the material. Accordingly, it appears interesting to design modeling tools that may address the material complexity and multi-scale organization.

Here, in order to bridge the descriptions of the self-healing phenomenon between fiber and tow scales (Figure 1c,d), we propose an image-based approach, in which the resolution domain is a 2D Finite Element (FE) mesh directly obtained from home-made mesh generation software.

The paper is organized as follows. First, the model is set up and translated into equations. Then, the employed numerical schemes are described, and their validation is discussed. Finally, an application example is presented in detail, to illustrate the capabilities of the modeling toolbox.

2. Model Setup

Though the material has a multi-scale character, we will focus here only on minicomposites, that is straight bundles of filaments with multilayer coatings. The concept of minicomposite has been thoroughly used to build up a comprehension of how a full composite works [12], but with a simplified mechanical stress state. Models addressing minicomposites usually take advantage of the separation of the space coordinates into longitudinal (z) and transverse ($x,y$), allowing for 1D models after homogenizing on a transverse section.

As far as CMCs are concerned, one of the most important mechanical test configurations is tensile stress along the fibers, since it leads to matrix crack opening and eventually to material failure. This will be the situation considered here.

2.1. Material Description

The transverse section of a typical SH-CMC minicomposite model is depicted in Figure 1d. Several domains are easily recognized:

• The SiC fibers, which are sensitive to oxidation through a subcritical crack growth behavior. Their diameter is in the range of 8–15 $\mathsf{\mu }$m.
• The thin Pyrocarbon (PyC) interphase coating the fibers. This layer acts a crack deviator, helping the development of multi-cracking, which gives the material its exceptional mechanical behavior [3,13].
• The multi-layer matrix containing two constituents: silicon carbide, which is nearly inert at the temperatures typical of civil aircraft applications (it begins to oxidize appreciably above 1000 ${}^{\circ }$C [14]), and boron carbide, yielding a liquid (glassy) phase upon oxidation [4,5].
• A number of cusp-shaped pores left by the production process. In particular, the PyC, SiC, and boron carbide layers are deposited by Chemical Vapor Infiltration (CVI) [15,16]. This process gives characteristic continuous layers growing in conformal shape with respect to the initial substrate (i.e., the fibers). The layer thicknesses are larger in the outer parts of the fiber bundles, because of the competition between chemical deposition reactions and gas diffusion [17]. As a result, some pores remain between the fibers after infiltration. These imperfections may be connected or not to the ambient atmosphere, depending on whether the crack network intercepts them.

2.2. Phenomena and Scales

Under tensile stress, the minicomposite will undergo matrix transverse multi-cracking [18,19]. In this case, cracks are mostly perpendicular to the fiber direction [20]. Some may exhibit a helicoidal structure [21], but the helix vertical spacing is much larger than the crack thickness. We can thus safely focus on transversal cracks and adjust any micro-mechanical model to the measured behavior through the value of the average spacing between two successive cracks. The thickness of the cracks is of the order of less than 1 $\mathsf{\mu }$m, i.e., 10${}^{-3}$–10${}^{-2}$ times the typical diameter of the minicomposite [20]. Even in the presence of these cracks, fibers continue to provide the material with mechanical continuity and tensile resistance [18]. Nevertheless, the cracks provide a path to oxidizing species toward the fibers [22,23]. This constitutes the root of the weakening of the fibers themselves, and thus of the failure of the material [24,25].

The key phenomenon controlling the material lifetime is oxidation [26]. The most sensitive component is the pyrocarbon coating, which undergoes active oxidation above 500 ${}^{\circ }$C [27]. At temperatures between 450 ${}^{\circ }$C and about 700 ${}^{\circ }$C, typical of civil craft engine operating conditions, the oxidation of the boron carbide matrix component is efficient enough to fill active cracks with a liquid boron oxide [28,29]. This is caused by the higher molar volume of the oxide compared to the carbide. At low enough temperatures, the flow of this liquid is dominated by viscous effects, with respect to which both inertial and capillary effects are negligible [26]. The oxide is thus pushed through the crack by its own growth. Oxygen diffusion through the oxide is considerably slower than in the gas phase. Typically, the time scale of the diffusion through the liquid is 10${}^5$-times larger than that into the gas [30]. Silicon carbide may also produce a solid silica layer upon oxidation [31]; however, the temperatures considered here (<650 ${}^{\circ }$C) are not high enough to allow the development of a significant silica thickness [32].

2.3. Problem Equations

The natural ingredients of the model are the concentration balance for gaseous oxygen and for all oxidized materials, namely, pyrocarbon, undergoing active oxidation and giving gaseous CO${}_2$:

$$\mathrm{C}\left(s\right)+{\mathrm{O}}_2\left(g\right)⟶{\mathrm{CO}}_2\left(g\right),$$

(1)

and boron carbide, yielding a liquid B${}_2$O${}_3$ oxide together with dissolved, then gaseous CO${}_2$ [33]:

$${\mathrm{B}}_4{\mathrm{C}}_{1.6}\left(s\right)+4.6{\mathrm{O}}_2\left(g\right)⟶2{\mathrm{B}}_2{\mathrm{O}}_3\left(l\right)+1.6{\mathrm{CO}}_2\left(g\right).$$

(2)

Note the stoichiometry of boron carbide, experimentally determined as in excess carbon with respect to the known B${}_4$C or B${}_{13}$C${}_2$ formulae. Transport equations for oxygen mass and for the liquid oxide mass were established by resorting to an asymptotic, crack-averaged, two-dimensional description. We exploited the fact that the ratio of the crack half-thickness e over the minicomposite diameter d was small. As recalled in Section 2.2, this ratio assumes in practice values $ϵ=e/d\in [{10}^{-3},{10}^{-2}]$. In our approach, the three-dimensional governing equations, e.g., Fick’s law for oxygen concentration, were integrated over the crack thickness, and higher order variations (in $ϵ$) were neglected. This procedure has two effects:

(i)

the model provides only an averaged two-dimensional description of the transport processes within the crack. The main unknowns are the local crack averaged values of the O${}_2$ concentration and of the oxide flow velocity and the local height of liquid oxide within the crack;

(ii)

mass fluxes associated with chemical reactions, normally appearing in three dimensions as boundary conditions on the upper/lower crack boundaries, are embedded in the equations as source terms, which arise naturally from the integration procedure.

The obtained model is defined on a domain decomposed into several sub-domains and boundaries, as illustrated in Figure 2:

• ${\Omega }_r$: reactive regions of the matrix (indicated by ${\varphi }_r\left(\mathbf{r}\right)=1$ if $\mathbf{r}\in {\Omega }_r$, 0 elsewhere),
• ${\Omega }_i$: inert regions (indicated similarly by ${\varphi }_i\left(\mathbf{r}\right)$),
• ${\Omega }_f$: fibers (indicated similarly by ${\varphi }_f\left(\mathbf{r}\right)$),
• $\partial {\Omega }_{\mathrm{ext}}$: the external boundary,
• $\partial {\Omega }_{\mathrm{in}}$: the internal boundaries delimiting macro-pores, and
• $\partial {\Omega }_{\mathrm{f}}$: the fibers’ boundaries.

The balance equation for the oxygen evolution, defined on the reactive and inert domains, reads as follows:

$${\partial }_t\left(h{C}_{{\mathrm{O}}_2}\right)+\nabla ·{\mathbf{F}}_{{\mathrm{O}}_2}=S_{{\mathrm{O}}_2}$$

(3)

where the effective height h available to oxygen transfer is computed as:

$$h={\varphi }_{l}e(x,y)+{\varphi }_g{h}_g,{\varphi }_l+{\varphi }_g=1$$

(4)

where the phase indicators ${\varphi }_l(x,y)$ and ${\varphi }_g(x,y)$ take the values 1 or 0 depending on whether the considered point is in a region occupied by the liquid plug (${\varphi }_l=1$) or not. The quantity $e(x,y)$ denotes the local thickness of the crack, and the values of $h_l$ and $h_g$ are determined by both chemical reactions and oxide flow, as we will shortly see. The crack height-averaged oxygen flux ${\mathbf{F}}_{{\mathrm{O}}_2}$ is given by:

$${\mathbf{F}}_{{\mathrm{O}}_2}=-({\varphi }_{l}e(x,y){\mathcal{D}}_l+{\varphi }_g{h}_g{\mathcal{D}}_g)\nabla C_{{\mathrm{O}}_2},$$

(5)

with ${\mathcal{D}}_l$ and ${\mathcal{D}}_g$ the oxygen diffusion coefficients in the liquid and gas phases, respectively, and with oxygen consumption driven by the (algebraic) source term:

$$S_{{\mathrm{O}}_2}=-{\varphi }_r{\displaystyle \frac{k_l}{4.6}}{C}_{{\mathrm{O}}_2}$$

(6)

with an effective reaction rate $k_l$ defined using the Deal–Grove model [34] as:

$${\displaystyle \frac{1}{k_l}}={\displaystyle \frac{1}{k_{{\mathrm{B}}_2{\mathrm{O}}_3}}}+{\displaystyle \frac{h_l}{{\mathcal{D}}_l}}.$$

(7)

Here, the liquid-phase diffusion coefficient ${\mathcal{D}}_l$ is related to the parabolic kinetic constant of oxidation. The evolution equation for the liquid height $h_l$, without flow, is an ordinary differential equation:

$$\frac{d{h}_l}{dt}=-2V_{m,{\mathrm{B}}_2{\mathrm{O}}_3}{S}_{{\mathrm{O}}_2}$$

(8)

Simultaneously, the height of available boron carbide decreases according to the following equation:

$$\frac{d{h}_{{\mathrm{B}}_x\mathrm{C}}}{dt}=V_{m,{\mathrm{B}}_x\mathrm{C}}{S}_{{\mathrm{O}}_2}$$

(9)

As a result, the available height for gas circulation above the liquid in the reactive zone obeys:

$$\frac{d{h}_g}{dt}=\left(2V_{m,{\mathrm{B}}_2{\mathrm{O}}_3}-V_{m,{\mathrm{B}}_x\mathrm{C}}\right)S_{{\mathrm{O}}_2}$$

(10)

The pyrocarbon consumed height obeys the following equation:

$$\frac{d{h}_p}{dt}=V_{m,\mathrm{C}}{k}_p{C}_{{\mathrm{O}}_2}$$

(11)

with an effective reaction rate $k_p$ defined as for the liquid:

$${\displaystyle \frac{1}{k_p}}={\displaystyle \frac{1}{k_{\mathrm{PyC}}}}+{\displaystyle \frac{h_p-e}{{\mathcal{D}}_g}}.$$

(12)

In order to close the model, we must specify:

• an initial condition, here uniform in space:

  $$C_{{\mathrm{O}}_2}(t=0)=C_{{\mathrm{O}}_2}^e,$$

  (13)
• a boundary condition at $\partial {\Omega }_{\mathrm{ext}}$, translating exchange with the outer atmosphere through a boundary layer with thickness ${\delta }_g$:

  $${\mathbf{F}}_{{\mathrm{O}}_2}·\mathbf{n}=\frac{e\mathcal{D}}{{\delta }_g}\left(C_{{\mathrm{O}}_2}-C_{{\mathrm{O}}_2}^e\right),$$

  (14)
• a boundary condition on $\partial {\Omega }_{\mathrm{in}}$, here chosen equal to Equation (14), i.e., assuming that the outer atmosphere also enters the intra-bundle pores,
• and a boundary condition on $\partial {\Omega }_{\mathrm{f}}$, representing the oxygen consumption by pyrocarbon oxidation:

  $${\mathbf{F}}_{{\mathrm{O}}_2}·\mathbf{n}=e_p{k}_p{C}_{{\mathrm{O}}_2}$$

  (15)

When flow occurs, the spreading of the oxide is accounted for by the following propagation equation:

$${\displaystyle \frac{\partial {\varphi }_l}{\partial t}}+\mathbf{u}·\nabla {\varphi }_l=0$$

(16)

in which the fluid velocity $\mathbf{u}$ appears. To be able to solve this equation, one needs to define or compute somehow the liquid velocity $\mathbf{u}$. In our modeling approach, we have considered a potential approximation for the flow, i.e., that there exists a function $\Phi (x,y)$ such that $\mathbf{u}=\nabla \Phi$. Indeed, the flow is very slow and fully laminar, so we can safely assume that it is irrotational. Starting from a classical 3D potential formulation in the (moving boundary) liquid domain, with mass fluxes given on the reactive boundaries, we arrive at a 2D asymptotic crack-averaged potential equation reading:

$$\nabla ·\left(h\nabla \Phi \right)=Q,$$

(17)

where all the mass fluxes on the 3D boundaries lead to the appearance of the source term Q condensing reading:

$$Q=-\left(2V_{m,{\mathrm{B}}_2{\mathrm{O}}_3}-V_{m,{\mathrm{B}}_x\mathrm{C}}\right)S_{{\mathrm{O}}_2}$$

(18)

Note that the main unknowns of the model now represent average values across the crack opening of the material concentrations (in particular oxygen) and of the flow speed computed as the gradient of the potential obtained by solving (18).

3. Numerical Solution of the Coupled System

3.1. Overall Coupled Strategy

The global algorithm is schematized in Figure 3. The resolution involves three blocks, one for initialisation and two solvers for elasticity and oxidation, respectively.

The initialization block generates two Finite Element (FE) meshes from the initial data containing the position and diameter of all fibers, plus the different matrix (inert and reactive) layers.

The first FE solution block addresses only mechanical equilibrium by solving the elasticity equations on a 3D mesh describing the fibers, the crack in the matrix, and introducing decohesion between fibers and matrix. The crack width at every point of the matrix domain (${\Omega }_i\cup {\Omega }_r$) is utilized in the next block.

The second block addresses oxidation and contains three subsets devoted respectively to the oxygen mass balance (diffusion and reaction, Equations (3)–(7), solid and liquid mass balances (8)–(12)), and liquid flow (Equations (16)–(18). The solution was obtained by chaining to each other the three subsets sequentially. The same 2D mesh was used for the solution of the three subsets. The crack width field was necessary for this block and was recovered from the first block. The output of the second block contained, among others, the total exposure of the fibers to oxygen, which was given as a criterion for life duration of the material.

Additionally, it is possible to restart the computation of the second block after a fiber breakage event. In this case, the mechanical equilibrium was updated and a new crack height field obtained, which may lead to a partial or complete opening of the crack. Then, the oxygen diffusion and oxide production resumed. The area corresponding to the broken fiber was integrated in the domain ${\Omega }_i$ (inert), but the former fiber boundary $\partial {\Omega }_f$ had still a reactive character, so that there was an “internal boundary” condition there:

$$\left[{\mathbf{F}}_{{\mathrm{O}}_2}\right]·\mathbf{n}=e_p{k}_p{C}_{{\mathrm{O}}_2}$$

(19)

where the brackets denote the flux discontinuity between exterior and interior sides of the boundary.

3.2. Mesh Generation

The starting point for the problem resolution was a 2D image, containing the fibers described by discs; their size distribution and space correlation function were retrieved from the analysis of transverse optical micrographs of a minicomposite or of a bundle inside a composite. Once the fibers were located, the multi-layer matrix was generated by geometric dilation, either homogeneous or biased with a lesser dilation factor in the bundle center. The output of matrix image analysis or of CVI modeling is also a possible alternative as a starting point.

Two meshes were generated: the first was a fine 2D mesh for the oxidation and liquid spreading part; the second one was a coarser 2D mesh of the same image, but extruded in the third direction, in order to perform elasticity computations. An initial fiber/matrix decohesion can be introduced over a given length $l_d$.

3.3. Space and Time Discretization

The species mass balances were discretized following a classical finite element scheme: Equations (3) and (7) were written in weak (integral) formulation, then the Green–Ostrogradsky formula was applied and the boundary conditions (14)–(15) inserted; the result was discretized over a basis function set and transformed into a linear algebraic system, the unknowns of which were the basis function parameters, chosen for convenience as the variable values at the integration points of the elements. Triangular P1 elements were chosen in this study.

The solution of the fluid spreading equations (Equations (16) and (17) or (17)) was achieved with an inexpensive approximate numerical scheme, consisting of distributing uniformly the amount of fluid created in any connected component of the liquid phase (a “droplet”) during a time increment $\Delta t$ on all neighboring nodes surrounding the existing droplet. The continuity Equation (17) in its integral form can be written as a mass balance over $\Delta t$:

$$\Delta t{\int }_{\partial {\Omega }_d/{\Omega }_f}h\mathbf{u}·\mathbf{n}dS=\Delta t{\int }_{{\Omega }_d\cap {\Omega }_r}QdV$$

(20)

The droplet ${\Omega }_d$ was discretized over triangles T with faces f. Accordingly, Equation (20) can be rewritten as:

$$\Delta t\sum _{f\in \partial {\Omega }_d\backslash {\Omega }_f}h\mathbf{u}·\mathbf{n}\left|f\right|=\Delta t\sum _{T\in {\Omega }_d\cap {\Omega }_r}Q\left(T\right)\frac{\left|T\right|}{3}$$

(21)

where $\left|f\right|$ and $\left|T\right|$ are the measures of the face f (length) and triangle T (area).

By analogy with this last equation, one can write:

$$\Delta t\sum _{f\in \partial {\Omega }_d\backslash {\Omega }_f}h\mathbf{u}·\mathbf{n}\left|f\right|=\sum _{i\in {\Omega }_{nd}\backslash {\Omega }_f}{\left(\Delta h_l\right)}_i\frac{\left|T_{fi}\right|}{3}$$

(22)

where ${\Omega }_{nd}$ designs the set of nodes i lying outside the droplet, but forming a triangle $T_{fi}$ with a face f of the droplet boundary $\partial {\Omega }_d$ (see Figure 4).

The height increase ${\left(\Delta h_l\right)}_i$ on the $i^{\mathrm{th}}$ mesh node neighboring the face f and lying outside of the droplet will be:

$${\left(\Delta h_l\right)}_i=\Delta t\times \frac{h\mathbf{u}·\mathbf{n}\left|f\right|}{{\sum }_{j\in \left[1;N_{nt,i}\right]}\left|T_j\right|}$$

(23)

where $N_{nt,i}$ is the number of triangles having node i as one summit, as illustrated in Figure 4.

The time discretization of Equation (3) followed a first-order implicit Euler scheme, which is simple but robust. The time integration of the ODEs (9)–(12) describing the evolution of the height variables has been achieved with a Newton–Raphson scheme. For the sake of convenience, the inverses of the heights has been chosen as the variables; an arbitrarily large value was chosen for the initialization.

The choice of the time step was rather difficult to address because of the existence of two largely different time scales in the liquid and the gas phases ($\Delta t_l=\frac{L^2}{D_l}$ and $\Delta t_g=\frac{L^2}{D_g}$, respectively, where L is the characteristic width of the crack) arising from the ratio $\frac{D_g}{D_l}\approx {10}^5$ between the diffusion coefficients in both phases. Choosing a time step linked to gas diffusion $t_g$ would lead to extremely slow computations. In a first attempt, it had been chosen to consider the diffusion Equation (3) as steady, while the time integration would only rely on Equations (9)–(12) [33], but after the creation of a liquid domain ($h_g=0$), this was no longer possible. Therefore, here, the initial time step was set to the geometrical average of both time scales:

$$\Delta t=\sqrt{\Delta t_l\Delta t_g}$$

(24)

and then, after the healing event occurred, the time step was set to $\Delta t_l$. Figure 5 shows that the choice of time step (Equation (24) allows having a solution within 5% accuracy with respect to the finest time step, with a $\approx 40\times$ solution time acceleration; doubling this time step still preserved a good quality of the results, even though the precision of the implicit Newton scheme did not guarantee it in principle.

The mechanical equilibrium equations were solved by a classical elastic finite element scheme over P1 tetrahedral elements; the matrix solution made use of the MUMPS (MUltifrontal Massively Parallel Sparse direct solver) linear algebra library [35].

4. Numerical Example

We will now present and discuss in detail an example, in order to illustrate the capabilities of the developed modeling framework. First, the construction of the resolution domain is addressed; then, all physico-chemical parameters will be given. Finally, the results will be displayed and discussed.

4.1. Problem Statement

The 2D resolution domain, illustrated in Figure 6a, is representative of a small minicomposite, with 60 fibers, the sizes and positions of which have been extracted from a transverse micrograph of an actual composite. The multilayer matrix was modeled as concentric layers around each fiber; the thicknesses of each layer is reported in

Table 1. Input data for the example case: geometrical parameters.

Parameter	Symbol	Value	Unit
Outer inert layer thickness		3.0	$\mathsf{\mu }$m
Outer healing layer thickness		2.0	$\mathsf{\mu }$m
Second inert layer thickness		3.0	$\mathsf{\mu }$m
Second healing layer thickness		6.0	$\mathsf{\mu }$m
Third inert layer thickness		3.0	$\mathsf{\mu }$m
Third healing layer thickness		6.0	$\mathsf{\mu }$m
Inner inert layer thickness		3.0	$\mathsf{\mu }$m
Pyrocarbon thickness	$e_p$	0.1	$\mathsf{\mu }$m
Average fiber diameter	$〈d_f〉$	13.3	$\mathsf{\mu }$m
Fiber diameter std. deviation	$\sigma \left(d_f\right)$	2.2	$\mathsf{\mu }$m
Fiber volume fraction	${\Phi }_f$	48.6	%
Minicomposite section	$S_0$	0.0167	mm${}^2$

. As previously mentioned, this 2D image has been converted to an FE mesh using a dedicated Python script and the Triangle software [36]. This mesh will serve as a numerical support for the physico-chemistry part of the simulation.

Then, by a simple mesh extrusion technique, a 3D mesh was produced to fully represent a segment of the minicomposite, and a matrix crack was introduced in the middle of the segment, as illustrated in Figure 6b). Note that, for the elastic computations, the matrix layers were not meshed individually, We rather considered a homogenized matrix with averaged properties (rendered as the gray domain in Figure 6b)). This allowed reducing the size of the computational mesh, without any significant influence on the results.

4.2. Geometrical and Physico-Chemical Parameters

All the parameters of the virtual material have been chosen as representative of a realistic situation, such as the ones analyzed in previous reports [26].

We report in particular the geometrical parameters characterizing the crack simulated in Table 1, while the parameters related to physics and chemistry are listed in

Table 2. Input data for the example case: physico-chemical parameters.

Parameter	Symbol	Value	Unit
Temperature	T	700	${}^{\circ }$C
Pressure	p	1.015	hPa
External oxygen mole fraction	$x_{O_2}$	0.21	-
Oxygen/air gas diffusion coefficient	${\mathcal{D}}_g$	$1.67\times {10}^{-4}$	m${}^2$ s${}^{-1}$
Oxygen/boria liquid diffusion coefficient	${\mathcal{D}}_l$	$3.75\times {10}^{-10}$	m${}^2$ s${}^{-1}$
External boundary layer size	$\delta$	10${}^{-1}$	m
Pyrocarbon oxidation rate constant	$k_{\mathrm{PyC}}$	$1.7\times {10}^{-3}$	m s${}^{-1}$
Boron carbide oxidation rate constant	$k_{{\mathrm{B}}_2{\mathrm{O}}_3}$	$1.0\times {10}^{-3}$	m s${}^{-1}$
Boria molar volume	$V_{m,{\mathrm{B}}_2{\mathrm{O}}_3}$	$4.49\times {10}^{-5}$	m${}^3$ mol${}^{-1}$
Boron carbide molar volume	$V_{m,{\mathrm{B}}_x\mathrm{C}}$	$2.33\times {10}^{-5}$	m${}^3$ mol${}^{-1}$
Pyrocarbon molar volume	$V_{m,\mathrm{C}}$	$7.06\times {10}^{-6}$	m${}^3$ mol${}^{-1}$

.

The oxygen/air gas diffusion coefficient has been evaluated thanks to the classical Chapman–Enskog relationship [37]. The oxygen/boria liquid diffusion coefficient and the boron carbide linear oxidation rate constant have been evaluated from experimental data on $B_{x}C$ oxidation [38].

4.3. Results and Discussion

4.3.1. Oxygen Diffusion and Plug Formation

The initial mechanical equilibrium computation yielded a 3D displacement field, from which can be extracted a 2D crack height field projected on the crack surface. This is illustrated in Figure 6c. Once this field was obtained, the oxidation computation took place. Figure 7 summarizes the evolution of the system during the first 5 h of exposition. The upper graph gives the oxygen concentration at three points pt1, pt2, and pt3, one on the external boundary, one in the middle layer, and the other close to a fiber (see also Figure 6a for a localization on the resolution domain). The inset of Figure 7b is a rendering of the concentration field at 0.05 h (i.e., 3 min), illustrating the important oxygen concentration gradient all around the same boundaries. The evolution of the concentration at point $pt3$ illustrates well the whole physico-chemistry at work during this stage: first, a very rapid decrease was linked to the primary oxygen consumption by the pyrocarbon interphase and to a lesser extent by the boron carbide layers; after 3 min, the decrease eventually stopped, because the glass formation had already managed to hinder partially the diffusion of gaseous oxygen through the formation of a continuous external barrier (see the bottom images of Figure 7); at $t=$ 18 min, the glass reached the fiber, and the oxygen concentration was consequently at its lowest. In the later period, up to 5 h, this concentration gradually increased again, because of diffusion through the liquid barrier.

Figure 8 is a map of the consumed pyrocarbon height, showing a marked gradient, not only from the outside to the interior of the minicomposite, but also on the boundary of the most external fibers, with a 1:4 ratio between the less and most exposed parts of the interphase. Even at equal distances from the outer boundary, significant differences appear between points of the interphase, depending on their position with respect to neighboring fibers.

4.3.2. After a Fiber Breakage Event: Crack Re-Opening and Re-Healing

To illustrate the coupling capabilities of the software, a fiber breakage event has been manually introduced at 22 min. The computation has been resumed following the procedure described above, i.e., considering the fiber interior as an inert, wettable domain and letting the pyrocarbon interphase proceed reacting. Figure 9 shows that after less than 7 min, the broken fiber had been fully invaded by the liquid. Figure 10 is a map of the liquid height at 29 min (end of the computation). This quantity was equal to the crack width in all inert domains and was superior to it in all reactive domains, with a marked gradient from the exterior (larger height) to the interior (smaller height); the zones in which the liquid barrier did not yet fill the crack completely are visible in the interior of the minicomposite section.

4.3.3. Towards Lifetime Prediction

So far, the presented model is able to simulate correctly the formation and spreading of the healing oxide in a given preexisting crack, with possible time variations of the crack width and morphology. However, the complete lifetime prediction still requires additional modeling elements, suitably coupled to the current oxide behavior model.

Indeed, the final lifetime duration is controlled by the diffusive processes of oxygen through the healing fluid, with the presence of strong gradients throughout the minicomposite, combined to some law describing the progressive strength reduction of the fibers under the effect of oxygen exposure. This can be addressed, e.g., using an oxygen-sensitive Subcritical Crack Growth (SCG) model, in which the classical SCG model parameters are made dependent on the cumulated oxygen exposure for every fiber [8,9]. Other models for SiC fiber behavior under oxygen exposure are now also available, accounting for phenomena like silica scale formation, crystallization, and residual stresses buildup or relaxation [39]. The presence of water vapor can also be handled in these models, with its main effect of volatilizing the oxides instead of having a liquid phase growth.

One of the most interesting expected outcomes of such a coupling is the possibility of accounting for the interplay between the initial strength distribution and the oxygen exposure gradients. It appears clearly from the simulations shown here that there are severe radial gradient of oxygen exposure inside the bundles (see for instance Figure 7); accordingly, it can be expected that the strength decrease will be larger on the bundle periphery. However, it is still unclear whether this has a strong impact on the life duration, as compared to the natural strength dispersion present in the fibers, which is classically represented by a Weibull distribution of strengths. A comparison of classical load-sharing models containing such a Weibull distribution and of models altered by the presence of an oxygen exposure gradient should be carried out in order to quantify this interaction.

5. Conclusions

A detailed modeling of the behavior of self-healing ceramic-matrix composite has been developed at the scale of a unidirectional bundle. This 2D image-based model takes precisely into account the inner structure of the composite. A series of original numerical tools has been chained in a general algorithm able to describe the progressive creation of the healing fluid, even after a fiber breakage event. A realistic application example has been shown and discussed. It clearly shows significant gradients in liquid height and pyrocarbon consumption throughout the minicomposite section. New insights can be obtained on the interplay of the mechanical and physico-chemical mechanisms at work in the material during its life in use. Many further work directions can be defined from this point. First, varying conditions would allow having a more precise view of the influence of every phenomenon. Second, coupling the current model to a fiber strength reduction model [8,9,39] would allow simulating the material degradation and lifetime and determining the most critical conditions for the composite life duration. This involves coupling the degradation models to mechanical equilibrium computations, allowing for a prediction of stress transfer upon each fiber breakage event until final failure of the bundle. Another improvement to the model is the incorporation of the influence of water vapor, which has the effect of partly volatilizing the boron oxide and reducing the efficiency of the self-healing mechanism [28,38]. Finally, the approach may be extended to other types of CMCs, with different phases, like silicon-infiltrated SiC${}_f$/SiC composites [40].

Such extensions to the presented approach are currently under investigation.