# Damage and lifetime modeling for structure computations

Pierre Ladevèze, Emmanuel Baranger, Martin Genet, Christophe Cluzel

► To cite this version:

Pierre Ladevèze, Emmanuel Baranger, Martin Genet, Christophe Cluzel. Damage and lifetime modeling for structure computations. Narottam P. Bansal; Jacques Lamon. Ceramic Matrix Composites: Materials, Modeling and Technology, Wiley, pp.465-519, 2014, 9781118231166. 10.1002/9781118832998.ch17. hal-01274916

HAL Id: hal-01274916

https://hal.science/hal-01274916

Submitted on 4 Oct 2021

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
# DAMAGE AND LIFETIME MODELING FOR STRUCTURE COMPUTATIONS

PIERRE LADEVEZE, EMMANUEL BARANGER, MARTIN GENET, AND CHRISTOPHE CLUZEL

LMT-Cachan (ENS Cachan / CNRS UMR 8535 / PRES UniverSud Paris), Cachan CEDEX, France

## 17.1 INTRODUCTION

Ceramic matrix composites (SiC fibers, PyC interphase, and mutistacked Si-B-C matrix in this paper) combine good mechanical properties with the capability of sustaining high temperatures in oxidizing environments. The use of a multi-layer self-healing matrix improves the lifetime of ceramic matrix composites (CMCs) considerably [1]. This technology along with the use of SiC fibers paves the way for applications with very long lifetimes (several thousands of hours).

Structural applications of CMCs need to answer to the main and critical question that is at the center of this chapter: how to predict, that is, to compute, the lifetime of engineering structures. The answer is a multiphysic macroscopic model able to reproduce couplings between the mechanical and environmental fields. For such materials, confidence in the simulation is related directly to one's understanding of each degradation mechanism and to the introduction of such mechanisms in the model; lifetimes being too large for a purely experimental validation over long lifetimes (tens of thousands of hours).

The first part of the chapter deals with the modeling of the different microcracking mechanisms. The key phenomenon is the deflection of the matrix cracks in the fiber matrix interphase [2–4]: this mechanical fuse protects the fibers from the cracks and enables the development of multiple cracks in the matrix, giving the composite some level of ductility, even though all of its constituents are brittle. The matrix cracking scenario has already been established for woven CMCs [5]. Since the matrix between the yarns is stiffer, but more brittle, and has larger pores than the yarns themselves, it begins to crack first. For a well-densified composite, the corresponding cracks are oriented by the main loading direction and are orthogonal to the applied stress direction, which is a main modeling difficulty. Then, when this crack network nearly reaches saturation, new cracks mainly appear in the matrix inside the yarns. The corresponding cracks are oriented by the fibers: they are orthogonal to the fibers in the longitudinal yarns, and parallel to the fibers in transverse yarns if the main load direction corresponds to the longitudinal yarns. In summary, there are three types of cracks in the damaged material: (i) inter-yarn cracks, (ii) intra-yarn longitudinal cracks, and (iii) intra-yarn transverse cracks.

Among composites, CMCs are the most difficult situation for damage modeling. Rather complex closure effects occur for the damage mechanism associated with microcracks when direction is orthogonal to the loading direction and, therefore, is not known *a priori*. Here, we are developing our answer with a general anisotropic damage theory including microcracks closure effects [6,7]. In order to differentiate between the crack opening and closing effects, which are very pronounced in CMCs, the strain energy is divided into three parts associated with damage due to pure traction,
damage due to traction and compression, and initial stiffness for compression stresses. First applications to CMCs and C/C were proposed in Reference 8. A major evolution of the model, suitable for composites with self-healing matrix was proposed in References 7,9, and 10. In this approach, the behavior is initially anisotropic and the kinematics of damage is not defined *a priori*; it is the damage evolution laws that define which components of the compliance tensor are affected. Each of the different mechanical degradation mechanisms—inter-yarn matrix cracking, intra-yarn matrix cracking, and wear of the fiber/matrix interface—is represented by a damage variable. Then, an evolution law is associated with each macroscopic damage variable, that is, with each mechanism. Finally, the global damage results from the accumulation of the various degradation mechanisms. A first identification for a CMC family has been done in Reference 11.

The second part of the chapter is devoted to the modeling of couplings with the environment characterized by temperature, oxygen pressure, and humidity. Self-healing mechanism, thanks to multilayer matrix is also taken into account. The idea, here, is to develop micro-multiphysic models that provide information for the computational macromodel. In fact the last version is very suitable for that [11, 12]. The various experimental results obtained by SAFRAN Snecma Propulsion Solide or by Reference 13 were used to identify the mechanical model. In this description of material cracking, the evolution laws and crack opening indicators are derived from assumptions about the physical cracking mechanisms. For any complex static or fatigue loading [14], mechanical macro-modeling enables the calculation of the opening indicator value for the cracks perpendicular to the fiber yarns. This indicator constitutes the link from the mechanical model to the oxidation model. The link from the oxidation model to the mechanical model relies on the fiber's strength evolution as well as a PyC damaging. Indeed, oxygen diffusion toward the fibers, the PyC oxidation, and the creation of an oxide plug that heal the yarn are directly related to the opening of the inter-yarn cracks [15]. SAFRAN Snecma Propulsion Solide proposes a range of self-healing, multilayer CMCs, one of which is CERASEPA410, whose model is presented below. In an oxidizing environment, certain matrix layers become oxidized and form a glass that fills the cracks and plays the role of a plug that slows the diffusion of oxygen [2]. A model of the crack opening, matrix healing, and fiber strength degradation mechanisms as functions of oxygen concentration proposed in Reference 12 are described here. The methodology leading to a macro-multiphysic model giving the lifetime is put forward.

The last part of the chapter deals with computational aspects. Such a damage model allows to predict the macro-crack's initiation and their propagation. Classical numerical difficulties are displayed and answers are given. To end, the lifetime model is applied to several structural tests [16]. The first test concerns a dog-bone specimen subjected to tension. The combined action of the mechanical and thermal fields is analyzed in both hot- and cold-grip experiments (i.e., homogeneous and heterogeneous temperature fields), and elements of validation are presented. The fact that the model is capable of predicting the location of the specimen's fracture zone is found to be very useful for the design of lifetime experiments on structures. The second structural test concerns a plate with a circular opening subjected to tension (open-hole test). This test enables us to assess the effectiveness of the strategy in the case of multiaxial stress fields with high stress gradients. Then, having established a certain degree of confidence in the strategy, we begin to outline a damage tolerance analysis through the prediction of a specimen's residual properties and lifetime after impact. This answers one of the main industrial concerns regarding the development of new composite structures [17, 18]. Of course, this model is also used to *a posteriori* analyze the lifetime for a known state in terms of stresses or strains and temperature. An illustration is also given here.

## 17.2 DAMAGE MODELING BASED ON AN ANISOTROPIC DAMAGE THEORY INCLUDING CLOSURE EFFECTS

The CMCs' damage models given here are based on an anisotropic and unilateral damage theory proposed in References 6,7, and 9.

### 17.2.1 General Notations

First-order tensors, or vectors (e.g., displacement), will be denoted by $\mathbf{u}$, $\mathbf{v}$, etc. Second order tensors, that is, linear operators over vectors, or matrices (e.g., strains and stresses), will be denoted by $\varepsilon$, $\sigma$, etc. Fourth order tensors, that is, linear operators over matrices (e.g., Hooke's operator, compliance operator) will be denoted by $\mathbf{K}$, $\mathbf{C}$, etc.

The Eulerian transpose will be denoted by $^{t}\mathbf{u}$, $^{t}\mathbf{v}$, etc., such that the scalar product of $\mathbf{u}$ and $\mathbf{v}$ is $^{t}\mathbf{u}\mathbf{v}$. If there is no ambiguity, this scalar product can also be written $\mathbf{u} \cdot \mathbf{v}$.
Work and power bring into play the trace operator:

$$ \mathrm{Tr}(\underline{\sigma} \underline{\epsilon}) = \sum_{i,j=1,2,3} \sigma_{ij} \epsilon_{ij} $$

where $\sigma_{ij}$ and $\epsilon_{ij}$ denote the components of the tensors $\underline{\sigma}$ and $\underline{\epsilon}$ in an orthogonal basis.

For working with anisotropic operators, it is suitable to use the “engineering notation,” which is introduced here for two-dimensional (2D) problems. The stress and strain are represented by

$$ \hat{\underline{\epsilon}} = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \sqrt{2}\epsilon_{12} \end{bmatrix} \qquad \hat{\underline{\sigma}} = \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sqrt{2}\sigma_{12} \end{bmatrix} $$

so that $\mathrm{Tr}(\underline{\sigma} \underline{\epsilon}) = t \hat{\underline{\sigma}} \hat{\underline{\epsilon}}$, and the Hooke operator by

$$ \hat{\underline{K}} = \begin{bmatrix} K_{1111} & K_{1122} & K_{1112} \\ & K_{2222} & K_{2212} \\ \text{sym} & & K_{1212} \end{bmatrix} $$

### 17.2.2 Thermodynamic Framework

#### 17.2.2.1 Energy Potential and State Laws
After choosing a damage kinematics, the construction of a mathematical model does not cause any particular difficulty as long as one uses the classical framework of thermodynamics, which can be found in Reference 19. The scale can be macro or meso or even micro. Small displacement assumption is used here.

Let us consider a damage kinematics defined by the sets of damage variables $d$ and $\delta$. $\underline{\epsilon}$, $\underline{\epsilon}_e$, and $\underline{\epsilon}_i$ designate the total, elastic, and inelastic strains, such that

$$ \underline{\epsilon} = \underline{\epsilon}_e + \underline{\epsilon}_i \qquad (17.1) $$

At time $t$ at point $M$ of the domain $\Omega$, the state of the material is assumed to be correctly defined by the values at time $t$ and point $M$ of $\underline{\epsilon}$, $\underline{\epsilon}_e$, $d$, $\delta$, $X$, where $X$ designates additional internal variables such as hardening variables.

The free energy $\rho\Psi$ has the following functional dependence:

$$ \rho\Psi(\underline{\epsilon}_e, d, \delta, X) \qquad (17.2) $$

To satisfy thermodynamic conditions, it must be continuously differentiable and convex. The damage forces conjugates of $d$ and $\delta$ are, for fixed $\underline{\epsilon}_e$ and $X$:

$$ \left\{ \begin{array}{l} \underline{Y}_d = -\rho \frac{\partial \Psi}{\partial d} \\ \\ \underline{Y}_\delta = -\rho \frac{\partial \Psi}{\partial \delta} \end{array} \right. \qquad (17.3) $$

In many cases, one also has

$$ \left\{ \begin{array}{l} \underline{Y}_d = -\rho \frac{\partial e}{\partial d} \\ \\ \underline{Y}_\delta = -\rho \frac{\partial e}{\partial \delta} \end{array} \right. \qquad (17.4) $$
where e is the strain energy density. One also has

$$ egin{cases} \frac{Y_X}{\underline{X}} = \rho \frac{\partial \Psi}{\partial \underline{X}} \\ \underline{\sigma} = \rho \frac{\partial \Psi}{\partial \underline{\epsilon}_e} \end{cases} \qquad (17.5) $$

**17.2.2.2 Damage Evolution Laws** The damage forces govern the evolution of damage and, subsequently, rupture. The corresponding laws are

$$ egin{cases} \underline{d}(t) = \underline{A}_d(Y_d(\tau), Y_\delta(\tau); \tau \le t) \\ \underline{\delta}(t) = \underline{A}_\delta(Y_d(\tau), Y_\delta(\tau); \tau \le t) \end{cases} \qquad (17.6) $$

that state that the values at time t of the damage variables depend on the past time history of the damage forces. $\underline{A}_d$ and $\underline{A}_\delta$ characterize the material and have to satisfy the following thermodynamic condition:

$$ \underline{Y}_d \cdot \underline{\dot{d}} + \underline{Y}_\delta \cdot \underline{\dot{\delta}} \ge 0 \qquad (17.7) $$

In practice, for composites subject to quasi-static loading, the values of $\underline{d}$ and $\underline{\delta}$ at time t depend on the maximum over $[0, t]$ of some combination of the elementary forces. For fatigue loading, one has

$$ egin{cases} \underline{d} = \underline{d}_s + \underline{d}_f \\ \underline{\delta} = \underline{\delta}_s + \underline{\delta}_f \end{cases} \qquad (17.8) $$

where $(\underline{d}_s, \underline{\delta}_s)$ designate the quasi-static part. The fatigue part $(\underline{d}_f, \underline{\delta}_f)$ follows classical fatigue damage models.

**Remark** It is possible to introduce a variable transformation for the damage variables $\underline{d}$ and $\underline{\delta}$. Thus, the conjugate forces will also be modified. Such a variable transformation does not affect the prediction of the material's behavior.

**Remark** Serious difficulties appear at the structural level concerning localization phenomena and final fracture. Their cause is now rather well understood: the classical damage models present serious shortcomings. For example, they do not contain classical fracture mechanics, which works quite well in many cases. However, let us note that such models provide the correct response before the initiation point for a structure subject to complex loading. One solution is given by the concept of localization limiter [20–23]. These are regularization procedures; the additional terms are built from a nonlocal approach or from a second-gradient approach. Viscoplastic regularization can also be introduced.

Our solution for composites depends on the scale chosen. The best approach, whenever possible, is to use what we call a “damage mesomodel.” It is a semi-discrete modeling approach in which the damage state is locally uniform within meso-constituents [24]. A more general approach is to introduce a damage model with a delay effect combined with a dynamic analysis. Then, a characteristic length—actually a characteristic volume with a strong physical significance—is introduced. Calculations with such mesomodels are given in References 25–27. We will come back to that problem in Section 17.4.

**17.2.2.3 Inelastic Strains: Coupling Damage with Plasticity (or Viscoplasticity)** Microdefects, that is, damage itself, lead to sliding with friction and, thus, to inelastic strains. More generally, it is possible to get coupling phenomena between damage and what we classically call plasticity or viscoplasticity. A method that seems to work quite well for composites is applying plasticity or viscoplasticity type models to quantities that are called effective (effective stress $\underline{\underline{\sigma}}$ and effective inelastic strain $\underline{\underline{\epsilon}}$) [28], which have to satisfy

$$ \mathrm{Tr}(\underline{\underline{\sigma}} \underline{\dot{\epsilon}}_i) = \mathrm{Tr}(\underline{\underline{\sigma}} \underline{\dot{\tilde{\epsilon}}}_i) \qquad (17.9) $$
All other internal variables are, in fact, effective quantities. Then, the dissipation condition becomes

$$ \mathrm{Tr}(\tilde{\sigma} \tilde{\dot{\epsilon}}_i) + \underline{Y_d} \cdot \underline{\ddot{d}} + \underline{Y_{\delta}} \cdot \underline{\dot{\delta}} \geq 0 \qquad (17.10) $$

Here, we do not follow the strain equivalence assumption advocated in Reference 29, which does not work, in general, for composites. A choice that works for many composite materials is

$$ \left\{ \begin{array}{l} \tilde{\sigma} = K_0 K^{-1} \sigma \\ \tilde{\dot{\epsilon}}_i = K^{-1} K_0 \dot{\epsilon}_i \end{array} \right. \qquad (17.11) $$

where $K$ designates the eventually damaged Hooke's tensor, and $K_0$ the initial Hooke's tensor.

In order to predict the state of the material and its time dependence, one uses the following state evolution law:

$$ \frac{d}{dt} \left[ \begin{array}{c} \tilde{\dot{\epsilon}}_i \\ -\underline{X} \end{array} \right] = \underline{B} \left[ \begin{array}{c} \tilde{\sigma} \\ \underline{Y_X} \end{array} \right] \qquad (17.12) $$

where $\underline{B}$ is an operator that depends on the material. It must be positive in order to satisfy the second law of thermodynamics. For noneffective quantities, this type of model is classical for plasticity and viscoplasticity.

### 17.2.3 A First CMC Damage Model

A first characteristic of ceramic composites is that very different responses can be observed for compression and traction loading due to the opening and the closing of microcracks. Another characteristic concerns the microcracking mechanisms: microcracks are not only parallel or orthogonal to the fibers, they can be orthogonal to the loading direction. The family of CMCs studied and identified here are woven SiC/SiC composites with self-healing matrix made by Snecma Propulsion Solide [18].

#### 17.2.3.1 Energy Potential
Modeling closure effects for isotropic materials or for microcracks having a fixed direction is not very difficult [6,30,31]. Similar models have been developed, in particular in References 32–34. The difficulty is to develop models that lead to continuous relations between stresses and strains in the unilateral case and that are also associated with a free energy. Another much more serious difficulty is to model microcracking mechanisms orthogonal to the maximum positive traction stress direction for anisotropic materials. The first step in that direction, introduced in Reference 6 and developed in Reference 25, is to express this unilateral character in terms of energy rather than by distinguishing between “traction” and “compression” behavior using stresses and strains. The latter, at any rate, would not make sense from a general point of view. We split the elastic energy of the undamaged material into three parts: a “traction” energy, a “compression” energy, and an energy identical for traction and compression behavior:

$$ e = e_t + e_c + e_{tc} \qquad (17.13) $$

Moreover, we introduce three sets of damage variables, $(d_t, \delta_t)$, $(d_c, \delta_c)$ and $(d_{tc}, \delta_{tc})$, which now represent the damage state of the material. We obtain

$$ e(\underline{\epsilon}_e, d_t, \delta_t, d_c, \delta_c, d_{tc}, \delta_{tc}) = e_t(\underline{\epsilon}_e, d_t, \delta_t) + e_c(\underline{\epsilon}_e, d_c, \delta_c) + e_{tc}(\underline{\epsilon}_e, d_{tc}, \delta_{tc}) \qquad (17.14) $$

For a large number of materials, and in particular for CMCs, damage does not increase under compression loading, therefore $d_c = \delta_c = 0$. And for many materials, one also has $d_{tc} = \delta_{tc} = 0$.
Following the damage geometrical description introduced in Reference 6, the simplest anisotropic damage model is defined by

$$
\left\{
\begin{aligned}
e_t &= \frac{1}{2E_0} \mathrm{Tr}(\underline{H} \langle \underline{\sigma} \rangle_+^H \underline{H} \langle \underline{\sigma} \rangle_+^H) \\
e_c &= \frac{1}{2E_0} \mathrm{Tr}(\langle \underline{\sigma} \rangle_- \langle \underline{\sigma} \rangle_-) \\
e_{tc} &= -\frac{\nu_0}{2E_0} (\mathrm{Tr}(\underline{\sigma})^2 - \mathrm{Tr}(\underline{\sigma}^2))
\end{aligned}
\right.
\qquad (17.15)
$$

where $E_0$ and $\nu_0$ are the initial Young's modulus and Poisson coefficient, and $\underline{H}$ a positive second symmetric operator directly linked to $\underline{d}_t$ and $\underline{\delta}_t$. Let us introduce the spectral decomposition of $\underline{\sigma}$ relative to $\underline{H}$, where eigenvectors and eigenvalues are denoted $\underline{T}_i$ and $\lambda_i$:

$$
\left\{
\begin{array}{ll}
\underline{\sigma} \underline{T}_i = \lambda_i \underline{H}^{-1} \underline{T}_i & \forall i \\
{}^t \underline{T}_i \underline{H}^{-1} \underline{T}_j = \delta_{ij} & \forall i,j
\end{array}
\right.
\qquad (17.16)
$$

**Definition** Taking the positive part of the eigenvalues, we get the positive part of stress:

$$
\langle \underline{\sigma} \rangle_+^H = \sum_i \langle \lambda_i \rangle_+ (\underline{H}^{-1} \underline{T}_i)^t (\underline{H}^{-1} \underline{T}_i) \qquad (17.17)
$$

This positive part depends on the stress but also on the damage state. (Note that if the stress is positive in the classical sense, one has $\langle\underline{\sigma}\rangle_+^H = \langle\underline{\sigma}\rangle_+ = \underline{\sigma}$. Designates the classical negative part of $\underline{\sigma}$.)

**Remark** An important property of this model is that the energy is positive and can be differentiated [6]:

$$
\delta e_t = \frac{1}{E_0} \mathrm{Tr}(\underline{H} \langle \underline{\sigma} \rangle_+^H \underline{H} \delta \underline{\sigma}) + \frac{1}{E_0} \mathrm{Tr}(\langle \underline{\sigma} \rangle_+^H \underline{H} \langle \underline{\sigma} \rangle_+^H \delta \underline{H}) \quad (17.18)
$$

**Remark** The previous model has been used for CMCs [25,35], as well as for other materials such as concrete [36,37] and metal alloys [38].

**Remark** Other attempts to model closure effects have been made for anisotropic materials [39–41]. It seems, in practice, that the thermodynamics damage force is defined under the assumption that microcracking mechanisms have fixed directions; a critical analysis is given in Reference 42. However, the recent model derived in Reference 43 is consistent; it is based on the homogenization of distributed directional cracks, which has common points with the directional approach introduced in References 7, 24, and 44 and extended in Reference 45.

17.2.4 An Advanced CMC Damage Model

17.2.4.1 Energy Potential The framework of this advanced version has been introduced in Reference 7. Its inter-
est is to facilitate and also to better take into account the damage mechanisms, which at the microscale are rather simple
to describe. The first identification of the mechanical behavior has been done in References 9 and 10. The main idea
is that the damage kinematics is not specified a priori; then the positive part of the stress is defined for any compliance
operator. The damage kinematics is, in fact, specified through the choice of the damage evolution laws.
The starting point is the energy splitting as previously:

$$ \left\{ \begin{aligned} e(\underline{\sigma}, \underline{C}, \underline{Z}) &= \frac{1}{2} \mathrm{Tr}(\underline{C} \langle \underline{\sigma} \rangle_{+}^{C} \langle \underline{\sigma} \rangle_{+}^{C}) + \frac{1}{2} \mathrm{Tr}(\underline{C}_{0} \langle \underline{\sigma} \rangle_{-}^{C_{0}} \langle \underline{\sigma} \rangle_{-}^{C_{0}}) + \frac{1}{2} \mathrm{Tr}(\underline{Z} \underline{\sigma} \underline{\sigma}) \\ \underline{C}(t=0) &= \underline{C}_{0} \\ \underline{Z}(t=0) &= \underline{Z}_{0} = \underline{K}_{0}^{-1} - \underline{C}_{0} \end{aligned} \right. \qquad (17.19) $$

where $\underline{C}$ and $\underline{Z}$ are two compliance operators that can vary. Several options exist for the choice of $\underline{C}_{0}$ and then $\underline{Z}_{0}$. In References 11 and 16 we used $\underline{C}_{0} = \underline{K}_{0}^{-1}$ so that $\underline{Z}_{0} = 0$. However, for an orthotropic material it is also possible to use (example in 2D with $\underline{n}_{1}$ and $\underline{n}_{2}$ the orthotropic axes):

$$ \underline{C}_{0} = \begin{bmatrix} 1/E_{1} & 0 & 0 \\ \text{sym} & 1/E_{2} & 0 \\ 0 & 1/2G_{12} & 0 \end{bmatrix}_{(\underline{n}_{1}, \underline{n}_{2})} \quad \text{and} \quad \underline{Z}_{0} = \begin{bmatrix} 0 & -v_{12}/E_{1} & 0 \\ \text{sym} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}_{(\underline{n}_{1}, \underline{n}_{2})} \qquad (17.20) $$

$\underline{C}$ has the classical restrictions:

$$ \left\{ \begin{aligned} \mathrm{Tr}(\underline{C} \underline{\sigma}_{i} \underline{\sigma}_{j}) &= \mathrm{Tr}(\underline{C} \underline{\sigma}_{j} \underline{\sigma}_{i}) \quad \forall \underline{\sigma}_{i}, \underline{\sigma}_{j} \\ \mathrm{Tr}(\underline{C} \underline{\sigma} \underline{\sigma}) &\ge 0 \quad \forall \underline{\sigma} \end{aligned} \right. \qquad (17.21) $$

but other than that no assumption is made about the compliance; one considers the most general damage kinematics by the internal variable $\underline{C}$. The first point is to consider that $\underline{C}$, rather than a conventional damage variable, is the internal variable defining damage. A second point, again according to the ideas developed previously, is to build the square root of $\underline{C}$.

**Definition** One then introduces the symmetric positive operator $\underline{H}$ such that $\underline{C} = \underline{H}^{2}$. $\underline{H}$ is easy to define by introducing the classical spectral decomposition of $\underline{C}$:

$$ \left\{ \begin{aligned} \underline{C} \underline{\sigma}_{i} &= \lambda_{i} \underline{\sigma}_{i} \quad \forall i \\ \mathrm{Tr}(\underline{\sigma}_{i} \underline{\sigma}_{j}) &= \delta_{ij} \quad \forall i, j \end{aligned} \right. \qquad (17.22) $$

Using the engineering representation, one gets

$$ \left\{ \begin{aligned} \hat{\underline{C}} \hat{\underline{\sigma}}_{i} &= \lambda_{i} \hat{\underline{\sigma}}_{i} \quad \forall i \\ \hat{\underline{\sigma}}_{i}^{\prime} \hat{\underline{\sigma}}_{j} &= \delta_{ij} \quad \forall i, j \end{aligned} \right. \qquad (17.23) $$

Thus, one gets

$$ \hat{\underline{H}} = \sum_{i} \sqrt{\lambda_{i}} \hat{\underline{\sigma}}_{i}^{\prime} \hat{\underline{\sigma}}_{i} \qquad (17.24) $$
**Definition** At this stage, we introduce the positive and negative parts of the stress:

$$ \begin{cases} \langle \underline{\sigma} \rangle_{+}^{C} = \underline{H}^{-1} \langle \underline{H} \underline{\sigma} \rangle_{+} \\ \langle \underline{\sigma} \rangle_{-}^{C_0} = \underline{H}_{0}^{-1} \langle \underline{H}_{0} \underline{\sigma} \rangle_{-} \end{cases} \qquad (17.25) $$

where $\langle \rangle_{+/-}$ designated the classical positive/negative part of a symmetric operator.

**Theorem** The traction energy is positive and can be differentiated [6, 7]:

$$ \delta e_t = \frac{1}{2} \mathrm{Tr}(\underline{C} \langle \underline{\sigma} \rangle_{+}^{C} \delta \underline{\sigma}) + \frac{1}{2} \mathrm{Tr}(\delta \underline{C} \langle \underline{\sigma} \rangle_{+}^{C} \langle \underline{\sigma} \rangle_{+}^{C}) \qquad (17.26) $$

**Remark** With $\underline{H}^2 = \underline{C}$, the traction energy could be written:

$$ e_t = \frac{1}{2} \mathrm{Tr}(\underline{C} \langle \underline{\sigma} \rangle_{+}^{C} \langle \underline{\sigma} \rangle_{+}^{C}) = \frac{1}{2} \mathrm{Tr}(\langle \underline{H} \underline{\sigma} \rangle_{+} \langle \underline{H} \underline{\sigma} \rangle_{+}) \qquad (17.27) $$

**Theorem** Another important result is that if $\underline{C}$ is positive, the traction energy is a convex function. Let us introduce $\underline{A} = \underline{H}\underline{\sigma}$. Thus, the traction energy writes $\mathrm{Tr}(\langle \underline{A} \rangle_{+} \langle \underline{A} \rangle_{+})/2$. The proof comes from the remark that the function $e_t(\underline{A})$ is the Legendre–Fenchel transform of the function

$$ \begin{array}{rcl} l: & \underline{\mathbb{E}} & \rightarrow \underline{\mathbb{R}} \\ & \underline{\epsilon} & \mapsto \frac{1}{2} \mathrm{Tr}(\langle \underline{\epsilon} \rangle_{+} \langle \underline{\epsilon} \rangle_{+}) + \Psi_K(\underline{\epsilon}) \end{array} \qquad (17.28) $$

and then a convex function. $\Psi_K$ is the indicator function of the convex set $K = \{\underline{\epsilon} / \inf_{\parallel \underline{V} \parallel = 1} \underline{V} \underline{\epsilon} \underline{V} \ge 0\}$, that is:

$$ \Psi_K(\underline{\epsilon}) = \begin{cases} 0 & \text{if } \underline{\epsilon} \in K \\ +\infty & \text{otherwise} \end{cases} \qquad (17.29) $$

We have

$$ l^*(\underline{A}) = \sup_{\underline{\epsilon} \in \underline{\mathbb{E}}} (\mathrm{Tr}(\underline{A}\underline{\epsilon}) - \frac{1}{2} \mathrm{Tr}(\langle \underline{\epsilon} \rangle_{+} \langle \underline{\epsilon} \rangle_{+}) - \Psi_K(\underline{\epsilon})) \qquad (17.30) $$

$$ = \sup_{\underline{\epsilon} \in \underline{\mathbb{E}}, \underline{\epsilon} \in K} \left( \mathrm{Tr}(\langle \underline{A} \rangle_{+} \underline{\epsilon}) + \underbrace{\mathrm{Tr}(\langle \underline{A} \rangle_{-} \underline{\epsilon})}_{\le 0} - \frac{1}{2} \mathrm{Tr}(\langle \underline{\epsilon} \rangle_{+} \langle \underline{\epsilon} \rangle_{+}) \right) \qquad (17.31) $$

$$ = \frac{1}{2} \mathrm{Tr}(\langle \underline{A} \rangle_{+} \langle \underline{A} \rangle_{+}) \qquad (17.32) $$

It follows that the admissible domain $\omega$ for the damage variable $\omega = (\underline{C}, \underline{Z})$ is such that the energy $e(\underline{\sigma}, \omega)$ is a convex function of $\underline{\sigma}$. Moreover, for $\omega \in \omega$, the energy is positive or equal to zero. The proof is trivial:

$$ \forall \underline{\sigma} \in \underline{\mathbb{E}} \quad e\left(\frac{\underline{\sigma} + 0}{2}, \omega\right) = \frac{1}{4} e(\underline{\sigma}, \omega) \le \frac{1}{2} e(\underline{\sigma}, \omega) + 0 \qquad (17.33) $$

$$ \Rightarrow \forall \underline{\sigma} \in \underline{\mathbb{E}} \quad 0 \le \frac{1}{4} e(\underline{\sigma}, \omega) \qquad (17.34) $$
It is also easy to see that $\omega$ is a convex set. For $\omega$ at the interior of $\omega$, $e(., \omega)$ is strictly convex. Moreover, from Equations (17.33) and (17.34) we get that $e(., \omega)$ is positive for $\underline{\sigma} \neq 0$.

At the boundary $\partial \omega$, $e(., \omega)$ is no more strictly convex and such $\omega$ states are limit damage states for the material. They are associated to the initiation of a macrocrack. A characteristic of these limit states is

$$ \exists \underline{\sigma} \neq 0 / e(\underline{\sigma}, \omega) = 0 \qquad (17.35) $$

**Remark** Let us define the energy in terms of the elastic strain $\underline{\epsilon}_e$. For $\omega \in \omega$:

$$ g(\underline{\epsilon}_e, \omega) = \sup_{\underline{\sigma} \in \mathcal{E}} (\mathrm{Tr}(\underline{\sigma} \underline{\epsilon}_e) - e(\underline{\sigma}, \omega)) \qquad (17.36) $$

All the properties of the energy $e(., \omega)$ are transfered to the energy $g(., \omega)$ for $\omega \in \omega$.

**Definition** A symmetric and positive compliance operator $\underline{C}$ is totally positive if $\underline{C} \underline{\sigma}$ is a positive operator for any positive $\underline{\sigma}$.

**Definition** A damaged compliance operator $\underline{C}$ is non-transverse relative to the direction $\underline{n}$ if $t^t(\underline{C} - \underline{C}_0)(n^t\underline{n})t = 0 \quad \forall t^t\underline{t}n = 0$.

### 17.2.4.2 Damage Kinematics

Damage evolution laws define the damage kinematics. Here, we consider situations in which it is possible to get information about the damage evolution law (e.g., Mode-I or Mode-II microcracking). In such situations, which are encountered very often, the approach described here leads to rather efficient damage models, of which the previous one related to a H-damage kinematics is a particular case. For the sake of simplicity, 2D problems and quasi-static loadings are considered.

Writing damage evolution laws in terms of stress or strain is mathematically equivalent. However, to the right instability behavior, it is easier to see them in terms of strains. That is why we have introduced some changes regarding our classical CMC model.

#### 17.2.4.2.1 Damage evolution laws describing microcracking mechanisms orthogonal to the loading direction (mechanism M0)

Let us consider that the dissipation due to damage is

$$ \dot{D} = \mathrm{Tr}(\underline{\sigma} \dot{\underline{\epsilon}}_e) - \dot{e}(\underline{\sigma}, \underline{C}, \underline{Z}) = \frac{1}{2} \mathrm{Tr}(\underline{\underline{C}} \langle \underline{\sigma} \rangle_+^C \langle \underline{\sigma} \rangle_+^C) + \frac{1}{2} \mathrm{Tr}(\underline{Z} \underline{\sigma} \underline{\sigma}) \qquad (17.37) $$

The second thermodynamic principle prescribes

$$ \dot{D} \geq 0 \qquad (17.38) $$

Many choices are possible for the damage forces that control the damage evolution. Here the strain is privileged and thus we take

$$ \left\{ \begin{array}{l} \underline{\dot{Y}'} = (\underline{\underline{C}}\langle\underline{\hat{\sigma}}\rangle_+^C)^t (\underline{\underline{C}}\langle\underline{\hat{\sigma}}\rangle_+^C) \\ \underline{\dot{Y}''} = (I_2\underline{\underline{C}}\langle\underline{\hat{\sigma}}\rangle_+^C)^t (I_2\underline{\underline{C}}\langle\underline{\hat{\sigma}}\rangle_+^C) \end{array} \right. \qquad (17.39) $$

where $I_2$ is the $\pi/2$ rotation operator.
Let us assume that the damage evolution depends on only one scalar effective damage force $\bar{z}$. Here, we introduce a mode-I microcracking by taking

$$ \left\{ \begin{aligned} \bar{z}(t) &= \sup_{\tau \le t} z(t) \\ z(t) &= (a \operatorname{Tr}(\hat{Y}')^{n+1} + (1-a) \operatorname{Tr}(\hat{Y}'^{n+1}))^{1/(n+1)} \end{aligned} \right. \qquad (17.40) $$

where a and n are two coefficients that depend on the material. The corresponding damage evolution laws are

$$ \left\{ \begin{aligned} \hat{C}^{-1} \hat{C} \hat{C}^{-1} &= \dot{\alpha} \frac{a \operatorname{Tr}(\hat{Y}')^n \mathbf{1} + (1-a) \hat{Y}'^n}{\bar{z}^n} \\ \hat{C}^{-1} \hat{Z} \hat{C}^{-1} &= \dot{\alpha} \frac{b \hat{Y}''^n}{\bar{z}^n} \end{aligned} \right. \qquad (17.41) $$

where b is a coupling coefficient that depends on the material. The function $\bar{z} \mapsto \alpha(\bar{z})$ is a material function. If a = 0 and n is large, $\hat{Y}'^n$ and $\hat{Y}''^n$ depend on the maximum positive strain and its direction. Then, the damage kinematics described by the evolution laws simulates a microcracking mechanism orthogonal to the maximum positive deformation direction. On the other side, if a = 1 and b = 0, one gets an isotropic damage.

The dissipation is equal to

$$ \dot{D} = \frac{\dot{\alpha}}{z^n} (a \operatorname{Tr}(\hat{Y}')^{n+1} + (1-a) \operatorname{Tr}(\hat{Y}'^{n+1}) + b \operatorname{Tr}((\hat{C}\hat{\sigma})\hat{Y}''^{n'}(\hat{C}\hat{\sigma}))) \qquad (17.42) $$

It is positive as long as $0 \le a \le 1$ and $b \ge 0$. Let us also note that the second is null if $\langle \underline{\sigma} \rangle_+^C = \underline{\sigma}$.

**Remark** Let us study the case where $\underline{n}_1$ and $\underline{n}_2$ denote the orthotropic axes of the material, and where the material has been damaged with a stress $\underline{\sigma} \underline{n}_1 \underline{n}_1$. It follows that the strain energy for any stress $\underline{\sigma}$ has the following expression:

$$ e = \frac{1}{2} \left( \frac{\sigma_{11}^+}{E_1^0} (1+h) + \frac{\sigma_{11}^-}{E_1^0} \right) + \frac{1}{2} \left( \frac{\sigma_{22}^+}{E_2^0} - \frac{2\nu_{12}}{E_1^0} \sigma_{11} \sigma_{22} \right) + \frac{1}{2} \frac{\sigma_{12}^2}{G_{12}^0} (1+\gamma) \qquad (17.43) $$

where h and $\gamma$ depend on the damage state. This expression characterizes a damaged material whose microcracks are parallel to $\underline{n}_2$.

**Remark** The compliance $\underline{C}$ is totally positive if $\underline{C}_0$ is totally positive. Under the hypothesis of the previous Remark, $\underline{C}$ is non-transverse relative to the direction $\underline{n}_2$.

### 17.2.4.2.2 Damage evolution laws describing microcracking mechanisms associated with a fiber direction (mechanism F1)
The damage mechanisms involved here are related to the yarn. Here, we describe a first model for which the damage evolution is controlled by only one effective damage force. More complex behaviors will be introduced in Section 17.3, devoted to the impact of environmental conditions.

Let $\underline{n}_1$ be the fiber direction and $\underline{n}_2$ the orthogonal direction. We define the projectors:

$$ \left\{ \begin{aligned} \underline{P}_1 &= \underline{n}_1' \underline{n}_1, \underline{P}_2 = \underline{n}_2' \underline{n}_2, \underline{P}_{12} = (\underline{n}_1' \underline{n}_2)_{\text{sym}} \\ \underline{Q}_1 &= \underline{P}_1' \underline{P}_1, \underline{Q}_2 = \underline{P}_2' \underline{P}_2, \underline{Q}_{12} = \underline{P}_{12}' \underline{P}_{12} \end{aligned} \right. \qquad (17.44) $$
The damage mechanisms associated with $\underline{n}_1$ are defined by the effective scalar damage force:

$$ \left\{ \begin{aligned} \bar{z}_1(t) &= \sup_{\tau \le t} z_1(t) \\ z_1(t) &= ((\mathrm{Tr}(\underline{Q}_1 \hat{Y}')) + c \mathrm{Tr}(\underline{Q}_2 \hat{Y}'))^2 + e \mathrm{Tr}(\underline{Q}_{12} \hat{Y}')^2)^{1/2} \end{aligned} \right. \qquad (17.45) $$

where c and e are two material constants. The corresponding damage evolution laws are

$$ \left\{ \begin{aligned} \underline{\hat{C}}^{-1} \underline{\hat{C}} \underline{\hat{C}}^{-1} &= \dot{\alpha}_1 \frac{\mathrm{Tr}(\underline{Q}_1 \hat{Y}') + c \mathrm{Tr}(\underline{Q}_2 \hat{Y}')}{\bar{z}_1} (\underline{Q}_1 + c \underline{Q}_2) \\ \underline{\hat{C}}^{-1} \underline{\hat{Z}} \underline{\hat{C}}^{-1} &= \dot{\alpha}_1 \frac{e \mathrm{Tr}(\underline{Q}_{12} \hat{Y}')}{\bar{z}_1} \underline{Q}_{12} \end{aligned} \right. \qquad (17.46) $$

where the function $\bar{z}_1 \mapsto \alpha_1(\bar{z}_1)$ is a material function. The dissipation for these mechanisms is

$$ \dot{D} = \frac{\dot{\alpha}_1}{2\bar{z}_1} ((\mathrm{Tr}(\underline{Q}_1 \hat{Y}')) + c \mathrm{Tr}(\underline{Q}_2 \hat{Y}'))^2 + e (\mathrm{Tr}(\underline{Q}_{12} \hat{Y}'))^2 \qquad (17.47) $$

It is positive as long as $c \ge 0$ and $e \ge 0$.

### 17.2.4.3 Residual Strains

The total strain is divided into elastic and inelastic parts:

$$ \dot{\underline{\epsilon}} = \dot{\underline{\epsilon}}_e + \dot{\underline{\epsilon}}_i \qquad (17.48) $$

where the inelastic part is due essentially to microcracking mechanisms, in fact sliding with friction. The evolution of the inelastic part is modeled thanks to two uncoupled plasticity models, one for each tow direction. For instance, regarding the longitudinal tow, the plasticity model is an associated one with isotropic hardening. For unilateral damage, the definition of inelastic strain rate is not clear. Then, we prefer to operate in a different way. For instance, for the direction $\underline{n}_1$, we introduce a stress criteria that can be differentiated:

$$ g(\underline{\sigma}) = \left( \left( \frac{(\hat{H}_0)_{11}^2}{(\hat{H}_{11}^2)} (\langle \underline{H}\underline{\sigma} \rangle)_+^2_{11} + (\langle \underline{H}\underline{\sigma} \rangle_-)^2_{11} \right) + d \left( \frac{(\hat{H}_0)_{11}^2}{(\hat{H}_{11}^2)} (\langle \underline{H}\underline{\sigma} \rangle)_+^2_{12} + (\langle \underline{H}\underline{\sigma} \rangle_-)^2_{12} \right) \right)^{1/2} \qquad (17.49) $$

Then, one has

$$ \dot{\underline{\epsilon}}_i = \dot{p} \frac{\partial g}{\partial \underline{\sigma}} = \frac{\dot{p}}{R(p)} \left( \left( \frac{(\hat{H}_0)_{11}^2}{(\hat{H}_{11}^2)} (\langle \underline{H}\underline{\sigma} \rangle)_+ \underline{H}P_1 + (\langle \underline{H}\underline{\sigma} \rangle_-) \underline{H}_1 \underline{P}_1 \right) + d \left( \frac{(\hat{H}_0)_{11}^2}{(\hat{H}_{11}^2)} (\langle \underline{H}\underline{\sigma} \rangle)_+ \underline{H}P_{12} + (\langle \underline{H}\underline{\sigma} \rangle_-) \underline{H}_1 \underline{P}_{12} \right) \right) \qquad (17.50) $$
with

$$ \left\{ \begin{array}{l} \dot{p} \ge 0 \\ g(\underline{\sigma}) - R(p) \le 0 \\ \dot{p}(g(\underline{\sigma}) - R(p)) = 0 \end{array} \right. \qquad (17.51) $$

The model depends on the material constant f and the material hardening function R : p ↦ R(p) with

$$ \dot{p} = \frac{\mathrm{Tr}(\underline{\sigma} \dot{\epsilon}_i)}{R(p)} \qquad (17.52) $$

**Remark** Inelastic behaviors such as viscoplasticity have been introduced in Reference 46. However, the previous plastic model can be easily extended to viscoplasticity.

**17.2.4.4 The Final Model** Residual strains are associated to a plasticity model with isotropic hardening. The internal variables are

* inelastic strain rate: $\dot{\epsilon}_i = \dot{\epsilon}_i^{n_1} + \dot{\epsilon}_i^{n_2}$;

* cumulated inelastic strain: $p^{n_1}, p^{n_2}$;

* inelastic threshold: $R^{n_1}, R^{n_2}$.

For an equilibrated material, the complete plasticity model depends on the hardening function R : p ↦ R(p) and the constant f.

The damage behavior is described by the damage variables $C$ and $Z$ such that

$$ \left\{ \begin{array}{l} C = C^{M0} + C^{F1} + C^{F2} \\ Z = Z^{M0} + Z^{F1} + Z^{F2} \end{array} \right. \qquad (17.53) $$

Each damage mechanism has its own variables. For an equilibrated material, the complete damage model for quasi-static loadings depends on two functions h and k and four material constants a, b, c, and e. These material characteristics do not depend upon environmental conditions.

For fatigue behavior, one uses what we have introduced previously:

$$ \left\{ \begin{array}{l} C = C^{\text{quasi-static}} + C^{\text{fatigue}} \\ Z = Z^{\text{quasi-static}} + Z^{\text{fatigue}} \end{array} \right. \qquad (17.54) $$

Damage laws will be seen in the Section 17.2.4.5 for static loadings and for fatigue 17.3.1.5 loadings.

Two rupture criteria ($\sigma_{n1}^f$ and $\sigma_{n2}^f$) are added to describe the brittle failure of the two yarn directions.

$$ \mathrm{Tr}(P_1 \underline{\sigma}) < \sigma_{n1}^f \qquad (17.55) $$

$$ \mathrm{Tr}(P_2 \underline{\sigma}) < \sigma_{n2}^f \qquad (17.56) $$

The two associated strengths depend upon environmental conditions (temperature, oxygen concentration, humidity, etc.), and that will be seen in Section 17.3.3.
**17.2.4.5 Identification at Room Temperature for Quasi-Static Loading** The material studied is woven SiC-SiC; a model using the H-damage kinematics was developed in References 25 and 35. According to the microscopic description derived in Reference 47, one can say that the first damage mechanism to appear is formed by microcracks starting from microvoids. The microcracks are orthogonal to the maximum positive stress direction. Toward the saturation of this mechanism and for a longitudinal stress, a new damage mechanism occurs in the transverse fiber bundle. Near the end, damage appears and increases in the longitudinal fiber bundle. The enhanced model, like the previous one, is able to simulate such phenomena. However, only the enhanced model is able to take into account strong initial anisotropy. Here, it has been identified for a particular composite [9]; one must add the damage contributions of the different damage mechanisms. The identified material functions at room temperature *h* and *k* are given in Figures 17.1 and 17.2. (For details on the calibration of the material parameters *a*, *b*, *c*, etc., see References 14 and 48.) Comparisons between tests and simulations are given in Figure 17.3.

**17.2.5 Some tools for the Implementation of the Anisotropic Framework**

Computation of the positive part of a symmetric operator $\underline{\sigma}$ Let us introduce the associated eigenvalue problem:

$$ \begin{cases} \underline{\sigma} T_i = \lambda_i T_i & \forall i \\ T_i \cdot T_j = \delta_{ij} & \forall i,j \end{cases} \qquad (17.57) $$

The positive part of $\underline{\sigma}$ is

$$ \langle \underline{\sigma} \rangle_+ = \sum_i \langle \lambda_i \rangle_+ T_i^T T_i \qquad (17.58) $$

where

$$ \langle \lambda_i \rangle_+ = \begin{cases} \lambda_i & \text{if } \lambda_i \ge 0 \\ 0 & \text{otherwise} \end{cases} \qquad (17.59) $$

<figure><img src="image_2.png" /><figcaption>**FIGURE 17.1**</figcaption></figure>

Identification of the material function *h* for the matrix damage mechanism.
<figure><img src="image_3.png" /><figcaption>FIGURE 17.2</figcaption></figure>

Identification of the material function k for the fiber bundle damage mechanism.

Computation of the square root of a symmetric operator $\underline{\underline{C}}$ Let us introduce the associated eigenvalue problem:

$$ \begin{cases} \underline{\underline{C}}\underline{\underline{\sigma}}_i = \lambda_i\underline{\underline{\sigma}}_i & \forall i \\ \hat{\underline{\underline{\sigma}}}_i \cdot \hat{\underline{\underline{\sigma}}}_j = \delta_{ij} & \forall i,j \end{cases} \qquad (17.60) $$

The square root of $\underline{\underline{C}}$ is

$$ \hat{\underline{\underline{H}}} = \sqrt{\hat{\underline{\underline{C}}} = \sum_i \sqrt{\lambda_i \hat{\underline{\underline{\sigma}}}_i^I \hat{\underline{\underline{\sigma}}}_i}} \qquad (17.61) $$

<figure><img src="image_4.png" /><figcaption>FIGURE 17.3</figcaption></figure>

Stress-strain response: tests and simulations. Reproduced with permission from Elsevier, Composites Science and Technology, V59, Iss7, 1999.
One also has

$$ \underline{H} \underline{\sigma} = \sum_{i} \sqrt{\lambda_{i}} \mathrm{Tr}(\underline{\sigma}_{i} \underline{\sigma}) \underline{\sigma}_{i} \qquad (17.62) $$

Inversion of the law $\epsilon_e = C_0 \langle \underline{\sigma} \rangle_+^C + C_0 \langle \underline{\sigma} \rangle_-^C + Z \underline{\sigma}$ Let us introduce the energy function:

$$ e(\underline{\sigma}) = \frac{1}{2} \mathrm{Tr}(C_0 \langle \underline{\sigma} \rangle_+^C \langle \underline{\sigma} \rangle_+^C) + \frac{1}{2} \mathrm{Tr}(C_0 \langle \underline{\sigma} \rangle_-^C \langle \underline{\sigma} \rangle_-^C) + \frac{1}{2} \mathrm{Tr}(Z \underline{\sigma} \underline{\sigma}) \qquad (17.63) $$

It is strictly convex as long as the damage state $\omega = (C, Z)$ belongs to the interior of $\omega$. It follows that the minimization problem for $\omega \in \omega^0$:

$$ \text{find } \underline{\sigma} \in \mathbb{E} \text{ minimizing } \left\{ \begin{array}{l} \mathbb{E} \to \mathbb{R} \\ \underline{\sigma}' \mapsto e(\underline{\sigma}') - \mathrm{Tr}(\underline{\sigma}' \underline{\epsilon}_e) \end{array} \right. \qquad (17.64) $$

admits a unique solution for any $\underline{\epsilon}_e \in \mathbb{E}$. Moreover, it is easy to see that the solution $\underline{\sigma}$ satisfies $\underline{\epsilon}_e = C_0 \langle \underline{\sigma} \rangle_+^C + C_0 \langle \underline{\sigma} \rangle_-^C + Z \underline{\sigma}$.

In practice, to get an approximated solution $\underline{\sigma}_h$, we propose an iterative technique to minimize $e$. This is a classical quasi-Newton method where the conditioner is

$$ \underline{\hat{R}} = \frac{1}{2}(\underline{\hat{C}} + \underline{\hat{C}}_0) + \underline{\hat{Z}} \qquad (17.65) $$

The algorithm is

Initialization: $\underline{\sigma}_0$

Iteration m: $\underline{\hat{R}}(\hat{\underline{\sigma}}'_m - \hat{\underline{\sigma}}_{m-1}) + \underline{\hat{C}}\langle \hat{\underline{\sigma}}_{m-1} \rangle_+^C + \underline{\hat{C}}_0\langle \hat{\underline{\sigma}}_{m-1} \rangle_-^C + \underline{\hat{Z}}\hat{\underline{\sigma}}_{m-1} - \underline{\hat{\epsilon}}_e = 0$

$$ \hat{\underline{\sigma}}_m = \lambda \hat{\underline{\sigma}}'_m \text{ with } \lambda = \begin{cases} \mathrm{Tr}(\hat{\underline{\sigma}}'_m \hat{\underline{\epsilon}}_e) / 2e(\underline{\sigma}'_m) & \text{if } \mathrm{Tr}(\hat{\underline{\sigma}}'_m \hat{\underline{\epsilon}}_e) > 0 \\ \mathrm{Tr}(\hat{\underline{\sigma}}'_m \hat{\underline{\epsilon}}_e) / 2e(-\underline{\sigma}'_m) & \text{if } \mathrm{Tr}(\hat{\underline{\sigma}}'_m \hat{\underline{\epsilon}}_e) < 0 \end{cases} $$

$$ \text{Stopping criterion: } \frac{t(\underline{\hat{R}}(\hat{\underline{\sigma}}_m - \hat{\underline{\sigma}}_{m-1})(\hat{\underline{\sigma}}_m - \hat{\underline{\sigma}}_{m-1}))}{t(\underline{\hat{R}}\underline{\hat{\sigma}}_m)(\underline{\hat{R}}\underline{\hat{\sigma}}_m)} \le \epsilon $$

## 17.3 MULTISCALE MODELING OF THE OXIDATION/DAMAGE COUPLING AND THE SELF-HEALING EFFECTS

As introduced before, many degradation and healing mechanisms can occur in the studied material. Handling the competition between these different mechanisms purely experimentally would lead to a very expensive campaign. In fact, the available information is spread over different scales and different physics. In order to have a complete overview of the mechanisms and their competitions, a model able to handle all the competitions is necessary. The framework of the proposed model [11] is depicted in Figure 17.4. It is a macroscopic model enriched with micro-descriptions of the competitions of groups of mechanisms. It relies on the global description of how the material works:

1. different cracks appear in the yarns and between the yarns;

2. oxygen diffuses quickly through the transverse intra-yarn crack network;
<figure><img src="image_5.png" /><figcaption>**FIGURE 17.4**
Modeling framework.</figcaption></figure>

3. the high concentration of oxygen attacks the fibers;

4. matrix layers react with oxygen and form a glass;

5. oxygen diffuses slowly through the glass;

6. the low concentration of oxygen attacks slowly the fiber;

7. the boron carbide layer oxidation decreases the glass plug efficiency;

8. the oxidation under stress of some fibers reduces the threshold of brittle yarns;

9. mechanical loading may lead to yarn breakage and leads to the failure of the composite.

The macroscopic mechanical model presented previously is the first brick of the model. It allows the description of matrix cracks and to separate the influence of each crack network. To this damage modeling is added a brittle failure stress criterion to model the fiber breakage. The other bricks will be described in the present section:

* use of a microcomposite model for the transverse intra-yarn cracks description (summarized as a crack opening indicator);

* use of a diffusion/reaction micromodel for the self-healing description;

* use of a cumulative oxygen concentration indicator for the fiber degradation description.

When available and reliable, the data identified at the microscale are used. When the microscopic data are not available directly due to the complexity of the experimental device or due to a lack of information like the statistical repartition of defects, the complementary microscopic data are identified by an inverse problem from experimental macroscopic tests. In order to increase the robustness of the predictions, the micromodels are used to give the form of the evolution laws and complete the missing macroscopic informations.

It is important to note that the model developed hereafter is used independently in the longitudinal and transverse directions.

### 17.3.1 Modeling of Fatigue and of the Transverse Intra-yarn Crack Opening

During the mechanical loading, different crack networks appear. For a given macroscopic loading, if the porosity is low, the main cracks first appear rather orthogonal to the loading direction. This first network is noted (M0) in Equation (17.53). It corresponds to cracks starting from macropore. In the second network, the crack directions are mainly dependent on the fiber directions. It is noted (Fi) in Equation (17.53). Among this second network, most of the matrix cracks in a given yarn, are orthogonal to the fiber direction. They deviate at the fiber/matrix interphase. This interphase made of Pyrocarbon layers (PyC) brings some ductility to the material, but those transverse cracks allow the oxygen to reach the fibers. The self-healing mechanism is based on the oxidation of boron carbide matrix layers, which leads to the formation of a higher volume of boron oxide used to fill the crack. The efficiency of the healing plug involves the volume of the cracks; therefore, an accurate description of the crack opening is a major key of such a modeling.

Two kinds of information are available on fatigue as well as crack opening of CMC yarns. The first kind is experimental data, almost inexistent for crack openings and for coupled static/fatigue loadings, some data are available on fatigue [13, 49, 50]. The other kind relies on models describing the behavior of a mini- or microcomposite at the fiber scale [3, 4, 51–53]. Our aim being to predict the behavior and lifetime at the macro-scale, the relevant experimental
<figure><img src="image_6.png" /><figcaption>**FIGURE 17.5**
Mini composite modeling for a null consumption of PyC.</figcaption></figure>

information relies at this scale. The models are used here to couple the available parsimonious information to get the form of the different evolution laws and indicators. These mechanical models are of great interest to handle properly the competition between different mechanisms, that is, transverse cracking, crack deviation at the PyC interphase, and wear of the interphase. The aim of such a detailed modeling is to define a relation between macroscopic damage variables and the transverse crack opening values for a given loading.

**17.3.1.1 Micromodel Description and Identification** Let us introduce the modeling of a microcomposite (a single fiber embedded in a matrix layer) as sketched in Figure 17.5. In this first model, the PyC consumption in oxidizing atmosphere is supposed to be neglectible. It is in agreement with models and experimental results in the healing loading range. The matrix and the fiber is denoted by *m* and *f* subcript letters. *L̃* is the crack spacing, which is a parameter in this model. Slipping occurs on a length *L̃*<sub>d</sub>. *σ* represents a homogenized stress, *u* a displacement, and *h̃* is the crack opening. *v*<sub>f</sub> and *v*<sub>m</sub> denote the fiber and matrix volume ratios in a yarn.

The fiber equilibrium local equation reads

$$ \pi \frac{d_f^2}{4} \frac{d\sigma_f}{dx} + \pi d_f \tau = 0 \quad \forall x \in [0, L_d] \qquad (17.66) $$

where $τ(x)$ represents the shear stress in the interface between fiber and matrix at position $x$ in the decohesion zone and $d_f$ the fiber diameter.

With small displacement assumption, the displacement of point O is given by

$$ -u_f(x=0) = \frac{\sigma}{E} \left(\frac{L}{2} - \tilde{L}_d\right) + \int_0^{\tilde{L}_d} \frac{\sigma_f}{E_f} dx \qquad (17.67) $$

$$ = \frac{\sigma}{E} \left(\frac{L}{2} - \tilde{L}_d\right) + \frac{\sigma \tilde{L}_d}{E_f v_f} - \frac{4}{E_f d_f} \int_0^{\tilde{L}_d} \int_0^x \tau(\alpha) d\alpha dx \qquad (17.68) $$

where

$$ E = E_m v_m + E_f v_f \qquad (17.69) $$

The strain is defined from the displacement:

$$ \tilde{\epsilon} = -\frac{u_f(0)}{\tilde{L}/2} \qquad (17.70) $$

$$ = \frac{1}{E}\sigma + \left(\frac{1}{E_f v_f} - \frac{1}{E}\right) \frac{\tilde{L}_d}{\tilde{L}/2}\sigma - \frac{4}{E_f d_f \tilde{L}/2} \int_0^{\tilde{L}_d} \int_0^x \tau(y) dy dx \qquad (17.71) $$

This stress/strain law can be rewritten in the same form as the macroscopic one:

$$ \tilde{\epsilon} = (C^0 \sigma + \Delta \tilde{C}_{stat+fat} \sigma + \tilde{\epsilon}_{stat+fat}^p) \qquad (17.72) $$
The stat and fat are relative to the static and fatigue damages in the yarns. The macroscopic maximal stress of the yarn is defined by

$$ \bar{\sigma}(t) = \sup_{\tau \le t} \sigma(t) \qquad (17.73) $$

The shear stress distribution at the interface when the microcomposite is unloaded is noted $\tau^p$. If the anelastic term $\bar{\epsilon}_{stat+fat}^p$ is defined as the residual strain after unloading, the identification gives

$$ C^0 = \frac{1}{E} \qquad (17.74) $$

$$ \Delta \tilde{C}_{stat+fat} = \frac{E_m v_m C^0 \tilde{L}_d}{E_f v_f \tilde{L}/2} - \frac{8}{E_f d_f \tilde{L} \bar{\sigma}} \int_0^{\tilde{L}_d} \int_0^x (\tau(\alpha) - \tau(\alpha)^p) d\alpha dx \qquad (17.75) $$

$$ \bar{\epsilon}_{stat+fat}^p = -\frac{8}{E_f d_f \tilde{L}} \int_0^{\tilde{L}_d} \int_0^x \tau(\alpha)^p d\alpha dx \qquad (17.76) $$

Like in Reference 51, the fatigue mechanism is attributed to the wear of fiber/matrix interface. To split damage parameters of Equation (17.75), the fatigue part corresponds to the variation of the shear stress, and the static part corresponds to creation of a new crack, or to the evolution of the length of the decohesion zone: $\tilde{L}_d$. To create a new crack, it is necessary to increase the loading. Then, the probability of transverse matrix crack creation is defined using the law:

$$ P = 1 - \exp \left( -2 \frac{S_m}{V_0} \left( \frac{\tilde{L}}{2} - \tilde{L}_d \right) \left( \frac{\sigma_m}{\sigma_0} \right)^n \right) \qquad (17.77) $$

The presented model is identified in order to represent the stress/strain behavior of a minicomposite (Figure 17.6). It has to be noted that for static loading without fatigue, the shear stress is a constant noted $\tau_0$ in the modeling.

### 17.3.1.2 Fatigue Modeling

**Fatigue damage evolution law** For a fatigue loading, the slip conditions at the PyC interface between fiber and matrix along $\tilde{L}_d$ decrease the shear stress. In this modeling, the wear of interface is not uniform, but it is assumed to be driven by the cumulative slip on each point of the interface. This driving parameter $L_{cs}$ is defined by

$$ L_{cs}(x) = \int_{time} |u_f(x) - u_m(x)| dt \qquad (17.78) $$

<figure><img src="image_7.png" /><figcaption>FIGURE 17.6</figcaption></figure>

Behavior of a minicomposite, comparison of the micromodel and experimental data.
<figure><img src="image_8.png" /><figcaption>FIGURE 17.7
Influence of the cyclic loading on the cumulative slip length.</figcaption></figure>

An example of evolution of $L_{cs}$ in case of a fatigue loading is given in Figure 17.7.

The wear effect is a progressive decrease of the interface shear stress $\tau$ for initial value $\tau_0$ to the final value $\tau_{min}$ such as

$$|\tau| = \tau_0 \frac{\exp(-\gamma_1 L_{cs}) + \beta_1}{1 + \beta_1} \quad (17.79)$$

$\tilde{L}_d$ is defined so that the stress in the fibers goes from $\frac{1}{v_f}\sigma$ to $\frac{E_f}{E}\sigma$ by stress transfer ($\tau$) on the matrix through the interface. In order to link the micromodeling to the macroscopic damage modeling, we define a mean shear stress by

$$\tilde{L}_d \tau_{mean} = \int_0^{\tilde{L}_d} \tau(y) dy = \left(\frac{1}{v_f} - \frac{E_f}{E}\right) \frac{d_f}{4} \sigma = \frac{E_m v_m}{E v_f} \frac{d_f}{4} \sigma = c \sigma \quad (17.80)$$

Let us introduce the following functional:

$$T(\tau) = \int_{0}^{1} \int_{0}^{\eta} \tau(\mu \tilde{L}_{d}) d\mu d\eta \quad (17.81)$$

With those functional, the damage of Equation (17.75) reads

$$\Delta \tilde{C}_{stat+fat} = \frac{E_m v_m}{E_f v_f} C^0 \frac{\tilde{L}_d}{\tilde{L}/2} - \frac{8}{E_f d_f \tilde{L} \sigma} \tilde{L}_d^2 T(\tau - \tau^p) \quad (17.82)$$

$$\Delta \tilde{C}_{stat+fat} = \left( \frac{E_m v_m}{E_f v_f} C^0 \frac{c}{\tau_{mean}} - \frac{4}{E_f d_f} \frac{c^2}{\tau_{mean}^2} T(\tau - \tau^p) \right) \frac{\bar{\sigma}}{\tilde{L}/2} \quad (17.83)$$

The function $\alpha$ is defined by

$$\Delta \tilde{C}_{stat+fat} = \alpha \frac{\bar{\sigma}}{\tilde{L}} \quad (17.84)$$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 17.8</figcaption></figure>

Behavior of a minicomposite: normal stress in fiber and matrix and shear stress in PyC interface for initial static loading and unloading.

Taking the time derivative to get the damage evolution law leads to

$$ \dot{\Delta \tilde{C}}_{stat+fat} = \left( \frac{\dot{\bar{\sigma}}}{\bar{L}} - \frac{\bar{\sigma} \dot{\bar{L}}}{\bar{L}^2} \right) \alpha + \frac{\bar{\sigma}}{\bar{L}} \dot{\alpha} \quad (17.85) $$

1. The case of a pure static loading corresponds to the creation of a new crack. The stress distribution in case of loading and unloading is described (Figure 17.8).

During the loading, $\tau_{mean} = \tau_0$ and $T(\tau - \tau^p) = \tau_0/4$ and then $\dot{\alpha} = 0$, which leads to

$$ \dot{\Delta \tilde{C}}_{stat}^{0} = \left( \frac{\dot{\bar{\sigma}}}{\bar{L}} - \frac{\bar{\sigma} \dot{\bar{L}}}{\bar{L}^2} \right) \alpha_0 \quad (17.86) $$

$$ \Delta \tilde{C}_{stat}^{0} = \frac{\bar{\sigma}}{\bar{L}} \alpha_0 \quad (17.87) $$

For this static case, the static damage increment is a function of the stress increment that leads to increase the decohesion zone $\tilde{L}_d$ and a function of the creation of new crack that involves $\tilde{L}$.

2. In case of fatigue after a static loading, $\dot{\bar{\sigma}} = 0$ and $\dot{\bar{L}} = 0$, this leads to

$$ \dot{\Delta \tilde{C}}_{fat}^{0} = \Delta \tilde{C}_{stat}^{0} \frac{\dot{\alpha}}{\alpha_0} \quad (17.88) $$

$$ \Delta \tilde{C}_{fat}^{0} = \Delta \tilde{C}_{stat}^{0} \frac{\alpha - \alpha_0}{\alpha_0} \quad (17.89) $$

In the general case

$$ \dot{\Delta \tilde{C}}_{stat+fat} = \Delta \tilde{C}_{stat}^{0} \left( 1 + \frac{\Delta \tilde{C}_{fat}^{0}}{\Delta \tilde{C}_{stat}^{0}} \right) + \Delta \tilde{C}_{stat}^{0} \frac{\dot{\alpha}}{\alpha_0} \quad (17.90) $$
In the micromodel, the effect of wear is modeled by Equation (17.79). The function $\frac{\alpha_0}{\alpha}$ is a function of the wear of fiber/matrix interface and its evolution follows the same kind of function. Introducing $\phi$ such as, for N fatigue cycles

$$ \phi(N) = \frac{\alpha_0}{\alpha} = \frac{\exp(-\gamma_2 N) + \beta_2}{1 + \beta_2} \qquad (17.91) $$

The macroscopic damage evolution law reads

$$ \frac{d\Delta\tilde{C}_{fat}}{dN} = -\Delta\tilde{C}_{stat}^0 \frac{\phi'(N)}{\phi(N)^2} \qquad (17.92) $$

Introducing the $\phi(N)$ function to describe the evolution of the interface wear leads to represent the evolution of $\alpha$ by a value of number of cycles $N$. The case of a static loading after $N_1$ cycles of fatigue leads to modify the $\tau_{mean}$ value and thus the $\alpha$ value. An equivalent number of cycles has to be recalculated by inversion of Equation (17.89).

**Fatigue inelastic strain evolution law** The inelastic part of the strain reads

$$ \tilde{\varepsilon}_{stat+fat}^p = -\frac{8}{E_f d_f \tilde{L}} \int_0^{\tilde{L}_d} \int_0^x \tau(\alpha)^p d\alpha dx \qquad (17.93) $$

$$ = -\frac{8\tilde{L}_d^2}{E_f d_f \tilde{L}} T(\tau^p) \qquad (17.94) $$

$$ = -\frac{8c^2}{E_f d_f} \frac{\bar{\sigma}^2 T(\tau^p)}{\bar{L} \tau_{mean}^2} \qquad (17.95) $$

$$ = \frac{\bar{\sigma}^2}{\bar{L}} \beta \qquad (17.96) $$

This leads to the following evolution law:

$$ \dot{\tilde{\varepsilon}}_{stat+fat}^p = \left( \frac{2\bar{\sigma}\dot{\bar{\sigma}}}{\bar{L}} - \frac{\bar{\sigma}^2\dot{\bar{L}}}{\bar{L}^2} \right) \beta + \frac{\bar{\sigma}^2}{\bar{L}} \dot{\beta} \qquad (17.97) $$

1. In case of a pure static loading without previous fatigue:

$$ \dot{\tilde{\varepsilon}}_{stat}^{p0} = \left( \frac{2\bar{\sigma}\dot{\bar{\sigma}}}{\bar{L}} - \frac{\bar{\sigma}^2\dot{\bar{L}}}{\bar{L}^2} \right) \beta_0 \qquad (17.98) $$

$$ \tilde{\varepsilon}_{stat}^{p0} = \frac{\bar{\sigma}^2}{\bar{L}} \beta_0 \qquad (17.99) $$

2. In case of a fatigue loading after a static loading:

$$ \dot{\tilde{\varepsilon}}_{fat}^p = \tilde{\varepsilon}_{stat}^{p0} \frac{\dot{\beta}}{\beta_0} \qquad (17.100) $$

$$ \tilde{\varepsilon}_{fat}^p = \tilde{\varepsilon}_{stat}^{p0} \frac{\beta - \beta_0}{\beta_0} \qquad (17.101) $$

In the general case

$$ \dot{\tilde{\varepsilon}}_{stat+fat}^p = \dot{\tilde{\varepsilon}}_{stat}^{p0} \left( 1 + \frac{\tilde{\varepsilon}_{fat}^p}{\tilde{\varepsilon}_{stat}^{p0}} \right) + \tilde{\varepsilon}_{stat}^{p0} \frac{\dot{\beta}}{\beta_0} \qquad (17.102) $$
<figure><img src="image_11.png" /><figcaption>FIGURE 17.9</figcaption></figure>

Crack density of a minicomposite, comparison of the micromodel and experimental data.

Introducing $\Phi$ such as, for N fatigue cycles

$$ \Phi(N) = \frac{\beta_0}{\beta} = \frac{\exp(-\gamma_3 N) + \beta_3}{1 + \beta_3} \qquad (17.103) $$

The $\alpha$ and $\beta$ functions are not the same, but as the source of fatigue is the same wear of interface, we assume that $\gamma_2 = \gamma_3$. The $\beta_2$ and $\beta_3$ parameters affect the level of minimum value of the fatigue functions; thus the two functions can be approximated by a single one and a proportionality factor $\varphi$:

$$ \left(\frac{1}{\Phi(N)}\right)' = \varphi \left(\frac{1}{\phi(N)}\right)' \qquad (17.104) $$

The fatigue anelastic evolution law reads

$$ \frac{d\bar{\varepsilon}_{fat}^p}{dN} = -\bar{\varepsilon}_{stat}^{p0} \frac{\Phi'(N)}{\Phi(N)^2} \qquad (17.105) $$

Introducing the $\Phi(N)$ function to describe the evolution of the interface wear leads to represent the evolution of $\beta$ by a value of number of cycles N. The case of a static loading after $N_1$ cycles of fatigue leads to modify the $\tau_{mean}$ value and thus the $\beta$ value. An equivalent number of cycles has to be recalculated by inversion of Equation (17.101).

### 17.3.1.3 Crack Density Indicator

$$ \tilde{d} = \frac{1}{\tilde{L}} \qquad (17.106) $$

Numerical simulations of the model are in good agreement with the available experimental results (see Figure 17.9).

The crack density does not have a simple form in the general case but in our special case, this law can be approximated as linear (Figure 17.10).

### 17.3.1.4 Crack Opening Indicator

$$ \frac{\tilde{h}}{2} = uf(0) - u_m(0) = \int_{0}^{\tilde{L}_d} \frac{\sigma_f}{E_f} dx - \int_{0}^{\tilde{L}_d} \frac{\sigma_m}{E_f} dx \qquad (17.107) $$

$$ = \int_{0}^{\tilde{L}_d} \frac{\sigma_f}{E_f} dx \left( 1 + \frac{E_f v_f}{E_m v_m} \right) - \tilde{L}_d \frac{\sigma}{E_m v_m} \qquad (17.108) $$
<figure><img src="image_12.png" /><figcaption>FIGURE 17.10</figcaption></figure>

Crack density of a minicomposite versus damage without fatigue.

Note that

$$ \int_{0}^{\tilde{L}_d} \frac{\sigma_f}{E_f} dx = (\Delta \tilde{C}_{stat+fat} \sigma + \tilde{\epsilon}_{stat+fat}^p) \tilde{L}/2 + C^0 \sigma \tilde{L}_d \quad (17.109) $$

Thus:

$$ \frac{\tilde{h}}{2} = ((\Delta \tilde{C}_{stat+fat} \sigma + \tilde{\epsilon}_{stat+fat}^p) \tilde{L}/2 + C^0 \sigma \tilde{L}_d) \left(1 + \frac{E_f v_f}{E_m v_m}\right) - \tilde{L}_d \frac{\sigma}{E_m v_m} \quad (17.110) $$

$$ \frac{\tilde{h}}{2} = \tilde{L}/2(\Delta \tilde{C}_{stat+fat} \sigma + \tilde{\epsilon}_{stat+fat}^p)(1 + \frac{E_f v_f}{E_m v_m}) - \tilde{L}_d \sigma \left(-\frac{1}{E_m v_m} + C^0 \left(1 + \frac{E_f v_f}{E_m v_m}\right)\right) \quad (17.111) $$

$$ \frac{\tilde{h}}{2} = \tilde{L}/2(\Delta \tilde{C}_{stat+fat} \sigma + \tilde{\epsilon}_{stat+fat}^p) \left(1 + \frac{E_f v_f}{E_m v_m}\right) \quad (17.112) $$

### 17.3.1.5 Macroscopic Identification
From that study at the microscale, only the form of the laws is kept as the scale transition and statistics mean definition is not obvious between the studied scale and the macroscale.

Using the Voigt notations, the link between each mechanism and yarn macroscopic damage functions in the fiber longitudinal direction is associated to

1. creation of new crack and increase of decohesion length $\tilde{L}_d$ represented by $\Delta \dot{C}_{stat 11}^0 = \dot{C}$. Its evolution law is defined in Chapter 17.2.4.2.

2. residual crack opening after unloading is represented by $\epsilon_{stat 1}^p$. Its evolution law is defined in Chapter 17.2.4.3.

3. fatigue loading wears the PyC interface and decrease the interface shear stress. This shear is not uniform at the microscopc scale but its effects are taken into account at the macroscopic scale by the following fatigue damage evolution and inelastic strain laws:

$$ \frac{d\Delta C_{fat 11}}{dN} = \Delta C_{stat 11} \left(-\frac{\phi'(N)}{\phi(N)^2}\right) \frac{\Delta\sigma_1}{\sigma_1^{max}} \begin{bmatrix} f & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad (17.113) $$

$$ \frac{d\epsilon_{fat 11}^p}{dN} = \epsilon_{stat 1}^p \left(-\frac{\phi'(N)}{\phi(N)^2}\right) \frac{\Delta\sigma_1}{\sigma_1^{max}} \begin{bmatrix} h \\ 0 \\ 0 \end{bmatrix} \quad (17.114) $$
where $\Delta\sigma_1$ is the longitudinal stress variation in the yarn and $\sigma_1^{max}$ the maximum value of this stress in a loading cycle. And the fatigue function defined by

$$ \phi(N) = \frac{\exp(-\gamma_0 N) + \beta_0}{1 + \beta_0} \qquad (17.115) $$

4. decrease of interface shear stress increases the length of fiber/matrix slipping $\tilde{L}$ by

$$ \Delta\dot{C}_{stat\ 11} = \Delta\dot{C}_{stat\ 11}^0 \left( \frac{\Delta C_{fat\ 11}}{\Delta C_{stat\ 11}^0} + 1 \right) \qquad (17.116) $$

5. when the fatigue loading has different phases with increasing stress levels, all cracks created do not wear out simultaneously. For that case, an additive fatigue damage law is defined based on Miner's law. When there is a variation of fatigue loading amplitude, an equivalent cycle number (in the sense of $\alpha$) is computed from

$$ \frac{\alpha}{\alpha_0} = \frac{\Delta C_{fat\ 11}}{f \Delta C_{stat\ 11}} \frac{\sigma_1^{max}}{\Delta\sigma_1} + 1 \qquad (17.117) $$

$$ \phi(N_{eq\ 1}) = \frac{\alpha_0}{\alpha} \qquad (17.118) $$

The damage macromodel allows to describe the state of degradation of each strand including a fairly accurate description of the longitudinal cracks. It is used to describe the opening of cracks, necessary for modeling the healing process. The macrocrack opening indicator reads

$$ \tilde{h}_1 = \tilde{L}_1 K (\Delta C_{stat+fat\ 11} \sigma_1 + \epsilon_{stat+fat\ 1}^p) \qquad (17.119) $$

where the macrocrack density indicator reads

$$ \tilde{d}_1 = \frac{1}{\tilde{L}_1} = a \Delta C_{stat\ 11} \qquad (17.120) $$

Similar laws are written in the direction 2. The introduced parameters are identified using fatigue tests at the macroscale and a unique value ($K$ over a parameter) of lifetime for the macrocrack opening indicator.

In order to validate the macroscopic representation, different cyclic loading paths (Figure 17.11) have been simulated using the micromodel. The comparison between the results of the micro- and macromodels are shown in Figure 17.12, which shows a good agreement.

### 17.3.2 Modeling of the Healing Process

The main effect of the presence of transverse cracks in the strands is to allow oxygen diffusion to fibers. At high temperature, silicon carbide fibers are sensitive to the presence of oxygen as well as the PyC interphase. The principle of self-healing is thus to limit oxygen concentration around the fibers by the introduction of a plug. The plug is made of glass resulting from the oxidation of matrix layers of boron carbide. When the volume of boron oxide is sufficient to fill the cracks in the yarns, the oxygen flow is limited by diffusion through the plug. The oxygen concentration around the fiber becomes low, which allows a significant increase in service life. Details of this principle and some morphology of self-healing C/SiC materials can be found in References 1,54 (see Figure 17.13). On that understanding of the healing mechanism, the simplified model presented in Figure 17.14 is deduced. It represents a group of several fibers protected by a self-healing matrix barrier on a perimeter noted $P_e$. Assuming a large radius of yarns in comparison to the plug length, the perimeter can be eliminated from the following equations. The material morphology and reaction kinetics are identified from References 55–57.

The diffusion model is solved using a wireframe linear discretization presented in Figure 17.15.
<figure><img src="image_13.png" /><figcaption>**FIGURE 17.11**
Cyclic loading scenarios.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 17.12**
Comparison between micro- and macromodels for various cyclic loadings.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE 17.13**
Morphology of self-healing matrix inside a yarn from Reference 1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 17.14</figcaption></figure>

Sketch of the multilayered matrix and healing scenarios.

**17.3.2.1 Healing Scenario Description** Healing is obtained by the oxidation of boron carbide matrix layers. In that case, the diffusion of oxygen through the oxide plug controls the plug creation and extension as the reaction kinetics are too fast. When the plug is active, the diffusion model allows to estimate the amount of oxygen around the fibers. To compute the diffusion across the plug with a complex geometry that changes during the simulation, a simple solution is to make a wireframe discretization of the plug. The relationship between wireframe and 3D models of the plug is based on the description of the flow of oxide. Thus, several volumes have to be defined.

The total volume of oxide in the plug is

$$V_{B_2O_3} = P_e(e_m e_{B_4C} + (L_{fi} + L_{fe})h) \qquad (17.121)$$

where $e_m \in [0, h/2 + h_m]$ is the thickness of the matrix plug, $e_{B_4C}$ the thickness of the healing matrix layer, $h$ the crack opening, and $e_{SiC}$ the total thickness of the SiC matrix layers, $L_{fi} \in [0, L_{fi0}]$ and $L_{fe} \in [0, L_{fe0}]$ the interior and exterior limit lengths of the oxide plug. $h_m$ is the length of oxidized $B_4C$ matrix. $S = h$. Perimeter is the crack surface perpendicular to the flow of glass (Perimeter being the perimeter of the healing layers around a group of fibers).

<figure><img src="image_17.png" /><figcaption>FIGURE 17.15</figcaption></figure>

Wireframe discretization of the model crack.
Let us introduce some volumes to define the considered healing scenario:

$$V_0 = (h + 2h_m)P_e e_{B_4C} \qquad (17.122)$$

$$V_1 = V_0 + 2L_{fi0}P_e h \qquad (17.123)$$

$$V_2 = V_0 + (L_{fi0} + L_{fe0})P_e h \qquad (17.124)$$

In order to calculate the flow of oxygen through the plug, a wireframe discretization is defined. But two cases are distinguished: the plug creation (case 1 of Figure 17.14) and the extension of the plug in the crack (case 2 of Figure 17.14).

1. $V_{B_2O_3} \le V_0$—closure of the crack by the oxide: G and D' are in B.

$$e_m = \frac{V_{B_2O_3}}{2P_e e_{B_4C}}; \quad L_{fi} = L_{fe} = 0 \qquad (17.125)$$

2. $V_{B_2O_3} > V_0$—extension of the plug within the crack: C is in B.

$$e_m = \frac{h}{2} + h_m \qquad (17.126)$$

a. $V_0 < V_{B_2O_3} \le V_1$—symmetric plug extension.

$$L_{fi} = L_{fe} = \frac{V_{B_2O_3} - V_0}{2P_e h} \qquad (17.127)$$

b. $V_1 < V_{B_2O_3} < V_2$—interior plug truncation.

$$L_{fi} = L_{fi0}; \quad L_{fe} = \frac{V_{B_2O_3} - V_0}{P_e h} - L_{fi0} \qquad (17.128)$$

c. $V_2 = V_{B_2O_3}$—interior and exterior plug truncation.

$$L_{fi} = L_{fi0}; \quad L_{fe} = L_{fe0} \qquad (17.129)$$

The plug truncation corresponds to the impossibility of the plug to grow as it reaches the fibers (interior truncation), or when it flows in the macro-porosity (exterior truncation).

Concerning the oxidative environment, the matrix cracks appear first in the composite. They are noted “M0” in Equation (17.53). They start from macropores and they cross the composite thickness up to the seal coat. These cracks are widely opened and for a temperature less than 700°C, the self-healing mechanism in the seal coat is not efficient. The composite being quite porous, the oxygen concentration in the pores is supposed to be equal to the concentration at the surface of the composite. While the plug in the seal coat is not efficient (cf. Chapter 17.3.4.2), the oxygen concentration at point G is related to the oxygen partial pressure by

$$C_G = \frac{P_{O_2}}{RT} \qquad (17.130)$$

The other boundary conditions are related to the chemical reactions at the PyC interface and $B_4C$ matrix layer. For the moment other healing matrix layers are not considered and the fibers are not supposed to react consequently with oxygen.

$$\Phi_{DE} = C_E P_e e_{pyc} k_{PyC} \left( = \frac{1}{2} \Phi_{BD} \right) \qquad (17.131)$$

$$C_F = 0 \qquad (17.132)$$

where $\Phi_{DE}$ is the flux of oxygen through wire DE, in (mol/s). The flux $\Phi_{BD}$ is used for robustness when $h_p = 0$.
The reaction constant $k_{PyC}$ is dependent on the temperature:

$$k_{PyC} = k_{PyC}^{0} \exp \left( -\frac{E_{PyC}}{RT} \right) RT \quad (17.133)$$

This law is identified from thermo-gravimetric analyses [58] depending on the oxygen partial pressure, this justify the second RT term appearing here.

The evolution of the diffusion coefficient of oxygen in $B_2O_3$ with temperature is given by References 59–63:

$$D_{ox} = k_p^0 \exp \left( -\frac{E_a}{RT} \right) RT \quad (17.134)$$

The parameter identification is obtained by oxidation tests of $B_4C$ pellets described in Reference 56.

The diffusion of oxygen in the air in a Knudsen and Fickian regime is given by Reference 64:

$$D_{gf} = \frac{C_M T^{1.75}}{P_{O_2}} \quad (\text{Fickian regime}) \quad (17.135)$$

$$D_{gk} = \frac{2h}{3} \left( \frac{8RT}{\pi M_o} \right)^{\frac{1}{2}} \quad (\text{Knudsen regime}) \quad (17.136)$$

$$D_g = \frac{1}{\frac{1}{D_{gf}} + \frac{1}{D_{gk}}} \quad (17.137)$$

where $C_M$ is a parameter function of air and oxygen molar mass, $M_o$ the molar mass of oxygen and $h$ the crack opening. The diffusion coefficient in the air must be defined because it allows a robust modeling of oxygen diffusion, particularly to the layer of PyroCarbon, but this is not a sensitive parameter of the model. Indeed, the diffusion rate in the $B_2O_3$ plug is very slow compared to diffusion rate in the air.

The description of the oxidation reaction of boron carbide is used to describe the flow of boron oxide formed as a function of oxygen flow. The main chemical balance equations [56] are

$$B_4C + 4O_2 \rightarrow 2B_2O_3 + CO_2 \quad (17.138)$$

$$B_4C + \frac{7}{2}O_2 \rightarrow 2B_2O_3 + CO \quad (17.139)$$

$$B_4C + 3O_2 \rightarrow 2B_2O_3 + C \quad (17.140)$$

$$C + O_2 \rightarrow CO_2 \quad (17.141)$$

C and CO are supposed to be recombined because O is in excess and so only the first equation is used.

The evolution laws for matrix consumption ($V_{B_4C}$), oxide creation ($V_{B_2O_3}$) and interface consumption ($V_{PyC}$) are function of oxygen flux. It reads

$$\dot{V}_{B_4C} = 2P_e e_{B_4C} \dot{h}_m = \frac{2}{4} V_{m.m} \Phi_{CF} \quad (17.142)$$

$$\dot{V}_{PyC} = \dot{h}_p P_e e_{PyC} = V_{m.pyc} \frac{\Phi_{BD}}{2} \quad (17.143)$$

where: $V_{m.m}$ and $V_{m.pyc}$ are, respectively, the molar volume of $B_4C$ matrix and PyC interphase in ($m^3/mol$). The volume of boron oxide used to create the plug is noted $V_{B_2O_3}$. It is both a function of the oxidation reaction of boron carbide and of oxide removal. The functions describing the evacuation of the oxide depend on the state of healing. The first case corresponds to the closure of the crack with a thickness of plug defined by Equation (17.125). The second case corresponds to the extension of the plug into the crack with sizes defined by equations from (17.126) to (17.129).
### 17.3.2.2 Case 1 Description
The evolution of the plug volume is

$$ \dot{V}_{B_2O_3} = 2V_{m.ox} \frac{\dot{V}_{B_4C}}{V_{m.m}} - Q_v P_e e_{B4C} \quad (17.144) $$

where $Q_v P_e h$ is the $B_2O_3$ oxide volatilization flux. This volatilization is activated by water partial pressure and its evolution on a free surface (in m/s) is given by Reference 65:

$$ Q_v = k_v \exp\left(-\frac{E_v}{RT}\right) P_{H_2O}^n \quad (17.145) $$

The different flux expressions on the branches of the circuit are

$$ \Phi_{ED} = D_g P_e e_{pyc} \frac{C_E - C_D}{h_p} \quad (17.146) $$

$$ \Phi_{BD} = D_g P_e h \frac{C_B - C_D}{L_{fi0}} \quad (17.147) $$

$$ \Phi_{BC} = D_g P_e e_{B4C} \frac{C_B - C_C}{h_m + \frac{h}{2} - e_m} \quad (17.148) $$

$$ \Phi_{CF} = D_{ox} P_e e_{B4C} \frac{C_C - C_F}{e_m} \quad (17.149) $$

The flux balance equations read

$$ 2\Phi_{ED} + \Phi_{BD} = 0 \quad (17.150) $$

$$ \Phi_{BC} - \Phi_{CF} = 0 \quad (17.151) $$

### 17.3.2.3 Case 2 Description
The evolution of the plug volume is

$$ \dot{V}_{B_2O_3} = 2V_{m.ox} \frac{\dot{V}_{B_4C}}{V_{m.m}} - Q_v P_e h - Q_{ej}(e_m - h_m)P_e \quad (17.152) $$

with $Q_{ej} = Q_{ej}^1 \Delta h + Q_{ej}^2$ being the flux of ejected oxide during a cyclic loading and caused by static flow of oxide for high temperature. Both $Q_{ej}^1$ and $Q_{ej}^2$ are parameters related to the viscosity of the oxide and to the wettability properties of the matrix by the oxide. But the modeling of the oxide by a fluid at microscale has so far not succeeded. Identification is performed in Reference 66 for one temperature by differential measurement of cyclic fatigue tests and creep tests. Associated to that law, the mean crack opening indicator over a loading cycle is used, allowing to treat complex fatigue cases. The different flux expressions on the branches of the circuit are

$$ \Phi_{ED} = D_g P_e e_{pyc} \frac{C_E - C_D}{h_p} \quad (17.153) $$

$$ \Phi_{DD'} = D_g P_e h \frac{C_D - C_{D'}}{L_{fi0} - L_{fi}} \quad (17.154) $$

$$ \Phi_{BD'} = D_{ox} P_e h \frac{C_B - C_{D'}}{L_{fi}} \quad (17.155) $$

$$ \Phi_{BF} = D_{ox} P_e e_{B4C} \frac{C_B - C_F}{h_m + \frac{h}{2}} \quad (17.156) $$
The flux balance equations read

$$2\Phi_{ED} + \Phi_{BD} = 0 \qquad (17.157)$$

$$\Phi_{BD'} + 2\Phi_{BF} + \Phi_{BG} = 0 \qquad (17.158)$$

$$\Phi_{BD'} - \Phi_{D'D} = 0 \qquad (17.159)$$

**17.3.2.4 Case Transition, Model Continuity, Integration Scheme, and Identification** Initially, $V_{B_2O_3} = 0$, $V_{B_4C} = 0$, and $V_{PyC} = 0$. The problem to solve is partially nonlinear and time dependent. The balance equations are solved analytically. The evolution laws are solved using a Runge-Kutta 45 algorithm and an adaptive time step. The continuity between the cases 1 and 2 is very important to ensure the convergency of the nonlinear solver. The model is driven by the volume of oxide for that purpose.

**17.3.2.5 Case of Cyclic Loading : Evacuation of Boron Oxide** When the plug heals the crack, if there is a variation of the stress related to mechanical fatigue loading, boron oxide present in the crack can migrate to a macro-pore. It is very complex to model at the microscopic level, the boron oxide flux leaving and returning to a macropore during a variation of crack opening. An average modeling of this mechanism is realized in Reference 66 for very high temperatures. The evacuated glass flow during a cycle is defined proportional to the change of opening by the equation

$$V_{ejection} = C_{ejec.}(h_{max} - h_{min}) \qquad (17.160)$$

where $C_{ejec.}$ is a material coefficient.

### 17.3.3 Modeling of the SubCritical Cracking of Fibers

Once the surrounding matrix has cracked, the fibers are no longer protected from the aggressive environment: air (oxygen, water, etc.) can penetrate the crack network and damage the fibers through oxidation or corrosion [67, 68]. There have been attempts to introduce an environmental protection of the fibers within the composite (using, for example, self-healing multilayered ceramic matrices or fiber treatments, referred to as chemical fuses), but even then the environmental protection is not complete and the fibers remain one of the key factors in the material's lifetime [11, 69]. The study of a fiber's lifetime under given environmental conditions is denoted static fatigue analysis because most works use the Paris law as in classical fatigue analysis, or subcritical cracking. Subcritical propagation, also known as slow crack growth, has been widely studied in the literature, especially for ceramics [67, 69–74]. These works describe the crack's propagation using a Paris-like law relating the crack's velocity to the crack's stress intensity factor calculated under the assumption of elastic behavior. This law, which depends on the temperature, agrees fairly well with experimental results [71, 72], but cannot be applied directly to varying environmental conditions such as oxygen or water concentration around the fibers. Moreover, it was found [67, 70, 72, 75] that as the temperature increases two different propagation stages take place. In the first stage, called the reaction-controlled stage, the crack's tip is damaged directly. The second stage, called the diffusion-controlled stage, is different in the sense that oxygen and water, prior to damaging the crack's tip, must first diffuse through an oxide layer formed on the surface of the fiber. This second stage only appears at high temperatures (above 750°C). Corresponding models can be found in Reference 70 for the fibers and in Reference 75 for the multifilament tow.

The subcritical cracking mechanism depends of the applied stress, the oxygen that has reached the fiber and is activated by temperature. Reference 69 proposes a simple model for constant loadings (applied stress, temperature, oxygen concentration). It links the yarn strength to the lifetime:

$$t\sigma_{yarn}^{n} = A_{0}\exp\left(\frac{E_{a}}{RT}\right) \qquad (17.161)$$

The lifetime is also inversely proportional to the oxygen partial pressure. This law is then simply rewritten introducing the cumulative oxygen concentration:

$$(C t)\sigma_{yarn}^{n} = B_{0}\exp\left(\frac{E_{a}}{RT}\right) \qquad (17.162)$$
From that equation, a failure criteria is deduced:

$$ \sigma < \sigma_{yarn}^{f} \qquad (17.163) $$

$$ \Theta(t, T) \sigma_{yarn}^{f} = B_{0} \qquad (17.164) $$

The strength of the yarn ($\sigma_{yarn}^{f}$) evolution is piloted by the cumulative oxygen concentration seen by the fibers (concentration at point D) and weighted by the temperature T:

$$ \Theta(t, T) = \int_{0}^{t} C_{D}(\tau) \exp \left( -\frac{E_{a}}{RT(\tau)} \right) d\tau \qquad (17.165) $$

This allows to follow the evolution of the strength for variable loadings of the fibers (temperature, healing, crack opening). All the information is condensed in $\Theta$. It can be written as an evolution law and integrated simultaneously as the diffusion/reaction model:

$$ \frac{d\sigma_{yarn}^{f}}{dt} = -\frac{1}{n} B_{0}^{-1} C_{D} \exp \left( -\frac{E_{a}}{RT} \right) (\sigma_{yarn}^{f})^{n+1} \qquad (17.166) $$

$$ \sigma_{yarn}^{f}(t=0) = \sigma_{yarn}^{f,static} \qquad (17.167) $$

$\sigma_{yarn}^{f,static}$ is the strength of the yarn in neutral environment.

To pass to the composite failure criterion, a constant scale factor is introduced, k:

$$ \sigma_{composite}^{f} = k \sigma_{yarn}^{f} \qquad (17.168) $$

$$ \frac{d\sigma_{composite}^{f}}{dt} = -\frac{k^{-n}}{n} B_{0}^{-1} C_{D} \exp \left( -\frac{E_{a}}{RT} \right) (\sigma_{composite}^{f})^{n+1} \qquad (17.169) $$

In practice, transverse cracks completely cross the yarns so that k is taken equal to the fiber volume ratio in the considered direction (either longitudinal or transverse).

The proposed evolution law is in agreement with the available experiments on yarns from Reference 69 at different temperatures and oxygen partial pressures. Nevertheless, one should note that in the available experiments the imposed stress is equal to the strength, that is, $\sigma = \sigma_{composite}^{f}$. As said before, the strength of the fiber must depend of the applied load as in References 76, 77. A variant of the previous evolution law can be postulated to take into account the influence of the applied stress and still respect the available experimental data; the strength is split into two contributions:

$$ \frac{d\sigma_{composite}^{f}}{dt} = -\frac{k^{-n}}{n} B_{0}^{-1} C_{D} \exp \left( -\frac{E_{a}}{RT} \right) (\sigma_{composite}^{f})^{n+1-m} \sigma^{m} \qquad (17.170) $$

In order to identify the exponent m, interrupted static fatigue tests can be used to get the residual strength of the fibers. Without available data, m = 0 is assumed.

### 17.3.4 Simulation of Degradation Mechanisms and Prediction of Lifetime

The first use of this modeling is to describe the competition of different degradation mechanisms leading to rupture of the composite under complex static or fatigue mechanical loadings and under various environmental conditions. To give a good confidence in the predictive simulations, some of the tests of static fatigue on short times are used for additional identification and another part with a different temperature range is used as validation. The associated identification and validation has been developed in Reference 48 for CERASEP A410.

#### 17.3.4.1 Identification
Concerning the fatigue modeling, cyclic fatigue tests have been conducted on coupons in Reference 13 and are directly used for identification of the macromodeling of fatigue.
Concerning the reaction/diffusion model, the experimental data from References 65 and 78 are used for the reaction kinetics, diffusion, and volatilization constants. The morphology is directly taken from micrographs, the sequential layers deposition being controlled and leading to a clear morphology definition. Concerning the crack opening indicator, only the form is kept and one parameter needs to be determined at the macroscale: (K) of Equation (17.119). For that, a static fatigue test is used. This test gives the lifetime of a composite for a constant 200MPa longitudinal loading at 500°C in air [49].

The ejection law due to fatigue described by Equation (17.160) can be identified [79] comparing two experiments: static fatigue lifetimes and cyclic fatigue lifetimes at the same mean stress. Both tests have been performed but not with the same mean crack opening; this does not allow us to identify this law properly. The identification is performed at three scales, micro, meso, and macroscopic.

* Microparameters: multi-layered matrix morphology, chemistry, subcritical crack growth of fibers/yarns;

* Mesoparameters: scale transition (k);

* Macroparameters: mechanical behavior for static and fatigue loadings based on static tests shown in Figure 17.3. The crack opening magnification (K) uses the mean lifetime value for a test under axial stress of 200 MPa, with a temperature of 500°C, 20 kPa of oxygen partial pressure and 1.8 kPa water vapor partial pressure. The corresponding mean lifetime of 250 hours can be seen via the set of tests in the middle of Figure 17.17.

### 17.3.4.2 Validation and Mechanisms Analysis

Further identification of the lifetime modeling for static fatigue loading is realized on coupons at 200 MPa and for a temperature of 500°C. The influence of stress and temperature is present in the whole model via the damage mechanisms (and corresponding crack networks), via the oxidizing of the matrix, via the oxide volatilization, and via the subcritical cracking of the fibers. It is thus interesting to compare the model response with different cases of loading (applied stress and temperature levels) for which experimental results are available (mainly in tension). From an experimental point of view, the test leading to the highest lifetime value (five months and a half or 4000 hours) corresponds to an applied stress of 100 MPa. The experimental lifetimes are quite short, these experiments are complex to control beyond few months while the expected lifetimes for industrial applications is about 5 years (4.4 10⁴ hours). For the sake of simplicity, the first comparison relates to the influence of the applied stress. The experiment/simulation comparison is presented in Figure 17.16 and shows the influence of the tension loading level for a constant temperature of 500°C and oxygen and water vapor partial pressures of 200 kPa and 1.8 kPa.

On the range of the experimental results, the model is in good agreement considering the mean slope as well as the highest variation for low stress levels. On the other hand, for applied stresses close to the static limit, the model predicts a sharp decline of the strength before the initial strength level.

It is interesting to note (in Figure 17.17) that this decline also appears on some experiments. This decline seems to be quite sensitive to the stress variation. The analysis of the mechanisms via the model shows, for high load levels, that the cracks are too widely opened to be healed correctly. To analyze this mechanism for a high load level, the

<figure><img src="image_18.png" /><figcaption>**FIGURE 17.16**

Lifetime versus applied stress for a temperature of 500°C, 20 kPa oxygen partial pressure and 1.8 kPa water vapor partial pressure.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**FIGURE 17.17**
Lifetime versus applied stress for different temperatures.</figcaption></figure>

simulations help to understand the protecting role of the oxide plug on the subcritical cracking of fibers. In Figure 17.18, several simulations have been conducted at 500°C under 20 kPa of oxygen partial pressure and 1.8 kPa water vapor partial pressure. The applied stress is here no more linked to the crack opening, this crack opening is imposed to three different values: low opening (h×1), wide opening (h×2), infinite opening (i.e., without possible healing). In Figure 17.18, the evolution of the failure stress threshold (the strength) of the longitudinal yarns is presented as a function of the exposure time. For example, for a 200 MPa loading, the lifetime corresponds to the intersection between the strength curve and the applied stress horizontal line. For the three prescribed crack openings, it corresponds to 160 hours (h×1), 40 hours (h×2), and 25 hours (without healing). The plateaux on the curves (h×1) and (h×2) correspond to the creation of a plug; as the crack opening increases, the oxide plug creation duration increases. The boron carbide consumption limits its capacity to regenerate, and beyond a certain point, the oxide plug loses its efficiency. The oxygen then reaches quite easily the fibers and the strength declines sharply (in log scale, nevertheless). Widely opened cracks lead to a fiber-driven lifetime while closed cracks lead to healing-driven lifetime. All curves merge for very long lifetimes. This plateau of plug efficiency is almost constant and for high stress levels, a small variation of crack opening (and thus of stress level) can switch a long lifetime to a short one.

<figure><img src="image_20.png" /><figcaption>**FIGURE 17.18**
Strength evolution with times for different crack openings.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIGURE 17.19</figcaption></figure>

Lifetime evolution for different partial pressures of water vapor.

Although macrotests are not available to determine the influence of oxygen and water vapor partial pressures, simulations allow the predictive quantification of their influence. For example, for an applied stress of 200 MPa and 500°C, Figure 17.19 shows the evolution of the lifetime versus the applied water vapor partial pressure. The negative influence of the water vapor is exhibited. In Figure 17.20, the influence of oxygen partial pressure is shown. It presents a non-monotonous evolution due to the competition of the healing efficiency and due to the fiber degradation rate.

In Figure 17.21, the evolution of the lifetime versus temperature can be observed. It is non-monotonic. For an applied stress of 150 MPa, a maximum lifetime can be found around 480°C, which is in agreement with experimental results on coupons [16].

For higher temperatures, the modeling with a plug localized at the heart of the material in the yarns is not sufficient. For high temperatures (around 700°C), a plug can form in the seal coat. Indeed, with thick layers of matrix in the seal coat, the cracks are too widely opened to allow the formation of a plug in the seal coat for temperatures between 450°C and 600°C. At 700°C the Boron oxide is sufficiently abundant to form a protection at the surface of the composite besides the protection at the heart of the composite. To understand this additional mechanism, the crack spacing has been measured directly on morphologies on composite coupons after loading as well as the matrix layers thicknesses in the seal coat. Simulations have been conducted in the presence of these two potential healing zones in a serial framework as shown in Figure 17.22.

Figure 17.23 shows the evolution of the oxygen concentration behind the plug in the seal coat (noted: seal-coat healing), which constitutes the entry for the modeling of the inner plug but also the oxygen concentration near the fibers (noted: yarn healing, Yarn + seal-coat healing). The figure shows that the first active plug occurs at the heart of

<figure><img src="image_22.png" /><figcaption>FIGURE 17.20</figcaption></figure>

Lifetime evolution for different partial pressures of oxygen.
<figure><img src="image_23.png" /><figcaption>**FIGURE 17.21**
Lifetime versus temperature.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**FIGURE 17.22**
Sketch of the double plug healing mechanism.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**FIGURE 17.23**
Influence of the seal coat on oxygen concentration at 700°C.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**FIGURE 17.24**
Influence of the seal coat on oxygen cumulative concentration at 700°C.</figcaption></figure>

the yarns. When the second plug becomes active, it creates a second protecting layer, which limits the yarn degradation in particular by limiting the water vapor content at the heart of the composite.

The simulation of the degradation of the strength of the yarns is obtained from a variable describing the cumulative oxygen concentration. Figure 17.24 shows the evolution of this variable. By fixing a limit of cumulative concentration corresponding to the failure at 150 MPa, this figure shows that the seal-coat plug does not allow an efficient protection by itself. Nevertheless, a combination with the inner plug allows to increase the lifetime of almost two decades.

The competition between the different mechanisms, fiber degradation, plug creation, and evacuation, is obtained from simulation. Figure 17.25 shows the strength evolution curves can cross because the different mechanisms are not affected in the same manner by the temperature. Without modeling these different mechanisms, it is very complex to interpret the experimental values. Indeed, the presented curves present the strength evolution and not only the final failure. The failure corresponds to the intersection of the strength curve and the applied stress history curve. For static fatigue tests, the loading is constant, but the model is able to handle more complex variable loadings.

From simulations, response surfaces can be built as shown in Figure 17.26 that presents the evolution of the lifetime versus applied stress and temperature for constant loadings. It allows to visualize the most pertinent loading ranges. Such a model based on the description of the micromechanisms allows: first to have a better reading and understanding of experimental values, and second, to optimize the material for large lifetimes applications.

<figure><img src="image_27.png" /><figcaption>**FIGURE 17.25**
Strength evolution with times for different temperatures.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>FIGURE 17.26</figcaption></figure>

Influence of the seal coat on lifetime maps.

## 17.4 PREDICTION CAPABILITIES

This section describes the implementation in a finite element code (in our case: ABAQUS/Standard) of the anisotropic and unilateral damage framework introduced in Section 17.2, with several illustrations of its suitability to predict the failure of industrial structures. Two different failure modes are considered here: failure in static fatigue, where the structure's lifetime depends on a coupled action of mechanics and environment; and static failure, where the environment plays a minor role, but where localized mechanical phenomena must be taken into account in order to satisfactorily predict the structure's behavior until final fracture. More details can be found in References 16, 80.

First, the numerical resolution of the model's local equations is described. It allows to compute the stress and damage fields over any structure and over time. And coupled to the environmental modeling introduced in Section 17.3, it allows to predict the lifetime of any structure in static fatigue.

Secondly, extension of the theory to the case of localized failure is described, together with associated numerical techniques. We illustrate the possibility to predict both the initiation and propagation of macroscopic cracks leading to the structure's failure.

### 17.4.1 Calculation of the Local Behavior

Classically, in order to compute mechanical fields over space and time in the framework of continuum mechanics, one usually makes the following assumptions: (i) deformation is small (< 1%), such that one can use a measure of the deformation that is linear with regard to displacement; and (ii) inertia effects are small, so that acceleration terms can be removed from the balance equation, making the problem quasi-static. This assumptions lead to the following problem:

$$ 
\left\{ 
\begin{array}{l} 
\text{find } \underline{u}(x, t), \underline{\epsilon}(x, t), \underline{\sigma}(x, t) \ \& \ \underline{K}(x, t) \ \text{over } \Omega \times [0; T] \ \text{such that} \\ 
\underline{\text{div}}_{\underline{x}}(\underline{\sigma}(x, t)) + \underline{f}^{g}(x, t) = 0 \quad \text{in } \Omega \times [0; T] \\ 
\underline{\sigma}(x, t) = \underline{K}(x, t) \underline{\epsilon}(x, t) \quad \text{in } \Omega \times [0; T] \\ 
\underline{\epsilon}(x, t) = (\underline{\text{grad}}_{\underline{x}} \underline{u}(x, t))_{\text{sym}} \quad \text{in } \Omega \times [0; T] \\ 
\underline{u}(x, t) = \underline{u}^{g} \quad \text{on } \partial\Omega^{u} \times [0; T] \\ 
\underline{\sigma}(x, t) \underline{n} = \underline{F}^{g} \quad \text{on } \partial\Omega^{F} \times [0; T] 
\end{array} 
\right. 
\qquad (17.171) 
$$

where $\Omega$ denotes the material domain and $[0; T]$ the time interval, $\underline{f}^g$ the applied volume load, $\underline{u}^g$ the imposed displacement on the domain boundary $\partial\Omega^u$, $\underline{F}^g$ the applied surface load on the domain boundary $\partial\Omega^F$ (one has
$\partial\Omega^u \cup \partial\Omega^F = \partial\Omega$ and $\partial\Omega^u \cap \partial\Omega^F = 0$), $\underline{n}$ the normal to the outward normal to the domain boundary, and $\underline{K}$ the eventually nonlinear relation between strains and stresses. See for instance, Reference 19 for more details.

Since for most real-world cases problem (17.171) does not have an analytical, or even closed form, solution (because of the complex material behavior as it is the case here, or simply because of the geometry), one has to use numerical methods to compute approximated solutions. If cutting-edge methods seem to tend toward so-called time-space approximations where the whole space and time domain is handled all along the computation in order to optimally capture the key features of the solution [81, 82], more traditional and widespread methods use the following approach. First, the time interval is discretized into several time increments (not necessarily of same duration), and one now looks for the solution at a finite number of time steps, from the initial one to the final one, assuming linear (or quadratic, etc.) evolution of the solution over the time increments. Second, in the same manner, for each time step the spatial domain is discretized into several finite elements (triangles, quadrangles, tetrahedrons, hexahedrons; one usually uses the same discretization for all time steps) over which the solution is assumed to be linear (or quadratic, etc.) and can be expressed from the solution at the element's nodes, so that one now looks for the solution at a finite number of points. This discretization of the equation generates a discrete problem from the initial continuous problem (17.171).

The space resolution for a given time step is usually based on Newton's method, which consists of iterating between two different problems: (i) resolution of global (i.e., structure-wide) equilibrium of the structure with external loading, the variables characterizing the material's behavior being fixed; and (ii) resolution of the local (i.e., point-to-point) material behavior, the displacement field being fixed. Most finite element packages available implement internally the time incrementation as well as the resolution of the structure's global equilibrium, as these are general mechanical equations. Thus, the only part that remains to implement is generally the local resolution of the material's behavior, which is strongly user- and application-dependent. We now describe such implementation for the anisotropic and unilateral damage model introduced in Section 17.2.

**17.4.1.1 Local Loop** The local loop is run for each integration point and at each global iteration of each load increment. The input is $\underline{\underline{\epsilon}}^{l,i}$, the total strain tensor at load increment $l$ and global iteration $i$, and the output consists of the damage tensors $\underline{\underline{C}}^{l,i}$ and $\underline{\underline{Z}}^{l,i}$, the corresponding damage “hardening” variable $\underline{\underline{\alpha}}^{l,i}$, and the stress tensor $\underline{\underline{\sigma}}^{l,i}$. Starting here, in order to simplify the notation, the subscripts $l, i$ will be omitted. The local problem writes

$$ \begin{cases} \text{with } \underline{\underline{\epsilon}}_e \text{ & } \underline{\underline{C}}_0 \text{ known, find } \underline{\underline{\sigma}}, \underline{\underline{C}} \text{ & } \underline{\underline{Z}} \text{ such that} \\ \underline{\underline{\epsilon}}_e = \underline{\underline{C}} \langle \underline{\underline{\sigma}} \rangle_+^C + \underline{\underline{C}}_0 \langle \underline{\underline{\sigma}} \rangle_-^{C_0} + \underline{\underline{Z}} \underline{\underline{\sigma}} \end{cases} \qquad (17.172) $$

**Fixed-point solver** We propose the following fixed-point algorithm.

$$ \left\{ \begin{array}{l} \text{Initialization: } j = 0; \alpha^j = \alpha^{l-1}; \underline{\underline{C}}^j = \underline{\underline{C}}^{l-1}; \underline{\underline{Z}}^j = \underline{\underline{Z}}^{l-1} \\ \text{Loop:} \\ \quad \begin{array}{l} \text{Stress: } \underline{\underline{\sigma}}^j / \underline{\underline{\epsilon}}_e = \underline{\underline{C}}^j \langle \underline{\underline{\sigma}}^j \rangle_+^C + \underline{\underline{C}}_0 \langle \underline{\underline{\sigma}}^j \rangle_-^{C_0} + \underline{\underline{Z}}^j \underline{\underline{\sigma}}^j \\ \text{Residual: } \underline{\underline{R}}^j = \alpha(\underline{\underline{\sigma}}^j) - \alpha^j \end{array} \\ \text{Exit test: } \frac{|\underline{\underline{R}}^j|}{|\alpha^j - \alpha^{l-1}|} < \text{tolerance} \Rightarrow \text{exit} \end{array} \right. \qquad (17.173) \\ \text{Damage: } \begin{cases} \underline{\underline{C}}^{j+1} = \underline{\underline{C}}(\underline{\underline{C}}^{l-1}, \alpha^{j+1} - \alpha^{l-1}) \\ \underline{\underline{Z}}^{j+1} = \underline{\underline{Z}}(\underline{\underline{Z}}^{l-1}, \alpha^{j+1} - \alpha^{l-1}) \end{cases} \\ \text{End loop: } j = j + 1 \end{array} $$

where $\underline{\underline{\alpha}}$ and $\underline{\underline{C}} \& \underline{\underline{Z}}$ are the material function and damage evolution laws defined in Equations (17.41) or (17.46). (In order to simplify the notation we presented the case of a single damage mechanism, but the extension to multiple mechanisms is straightforward.)
**Solver acceleration** The fixed-point method can be viewed as a Newton–Raphson method with a unit search direction and, therefore, can oscillate greatly. This is particularly true in our case because of the presence of unilateral conditions. Therefore, we propose to use Aitken's relaxation methods in order to improve the algorithm's convergence. This can be viewed as a search direction optimization based on previous iterations. The only modification to algorithm (17.173) concerns the damage increase line, which becomes

$$ \alpha^{j+1} = \alpha^j + s^j R^j \qquad (17.174) $$

$$ \text{with } s^j = \begin{cases} 1 & \text{if } j = 0 \\ -s^{j-1} \frac{R^{j-1}}{R^j - R^{j-1}} & \text{if } j > 0 \end{cases} $$

**17.4.1.2 Behavior Loop** In the previous discussion of the local loop, we did not described how we calculate $\underline{\sigma}$ when all the internal variables are fixed (second line of algorithm (17.173)). While this step is straightforward for most existing models, it is not for ours. Indeed, because of the partitioning of $\underline{\sigma}$ into positive and negative parts, the state law is nonlinear even with fixed operators, so the problem can be formulated as follows:

$$ \begin{cases} \text{with } \underline{\epsilon}_e, \underline{C}, \underline{C}_0 \text{ & } \underline{\sigma} \text{ known, find } \underline{\sigma} \text{ such that} \\ \underline{\epsilon}_e = \underline{C} \langle \underline{\sigma} \rangle_+^C + \underline{C}_0 \langle \underline{\sigma} \rangle_-^{\underline{C}_0} + \underline{Z} \underline{\sigma} \end{cases} \qquad (17.175) $$

**Newton–Raphson solver** This nonlinear problem is solved using a Newton–Raphson algorithm.

$$ \left\{ \begin{array}{ll} \text{Initialization: } k = 0; \underline{\sigma}^k = \underline{\sigma}^{k-1} & \\ \text{Loop:} & \\ \quad \text{Residual: } \underline{R}^k = \underline{\epsilon}_e - \underline{C} \langle \underline{\sigma}^k \rangle_+^C - \underline{C}_0 \langle \underline{\sigma}^k \rangle_-^{\underline{C}_0} - \underline{Z} \underline{\sigma}^k & (17.176) \\ \quad \text{Exit test: } \| \underline{R}^k \| < \text{tolerance} \Rightarrow \text{exit} & \\ \quad \text{Stress: } \underline{\sigma}^{k+1} = \underline{\sigma}^k + \underline{D}^k \underline{R}^k & \\ \text{End loop: } k = k + 1 & \end{array} \right. $$

for which several search directions $\underline{D}^k$ can be used. The actual tangent direction is not an option because one cannot derive the state law with respect to $\underline{\sigma}$ in the general case where $\langle \underline{\sigma} \rangle_+^C \neq 0$ and $\langle \underline{\sigma} \rangle_-^{\underline{C}_0} \neq 0$. One can choose, for example, the initial operator $\underline{K}_0 = \underline{C}_0^{-1}$, which is not a good direction, especially when the damage is significant, but which is fast because the calculation of the operator is very inexpensive. Another option is to use the quasi-secant operator, defined as

$$ \underline{D}^k = \begin{cases} (\underline{C}^{k-1} + \underline{Z}^{k-1})^{-1} & \text{if } \text{Tr}(\underline{\sigma}^k) > 0 \\ (\underline{C}_0 + \underline{Z}^{k-1})^{-1} & \text{if } \text{Tr}(\underline{\sigma}^k) < 0 \end{cases} \qquad (17.177) $$

which is a much better search direction, but also a more expensive one because these operators are usually not saved and must be recalculated at each iteration.

**Solver acceleration** In practice, the initial operator converges very slowly, or even does not converge at all if the damage is significant. The secant operator defined in Equation (17.177) generally converges very poorly, too. Therefore,
we propose once again to use Aitken's relaxation. The only modification to algorithm (17.176) is the last line, which becomes

$$ \underline{\sigma}^{k+1} = \underline{\sigma}^{k} + s^{k} \underline{D}^{k} \underline{R}^{k} \qquad (17.178) $$

$$ \text{with } s^k = \begin{cases} 1 & \text{if } k = 0 \\ -s^{k-1} \frac{\text{Tr}(\underline{R}^{k-1}(\underline{R}^k - \underline{R}^{k-1}))}{\text{Tr}((\underline{R}^k - \underline{R}^{k-1})(\underline{R}^k - \underline{R}^{k-1}))} & \text{if } k > 0 \end{cases} $$

### 17.4.2 Lifetime Predictions and Application to Damage Tolerance Analysis

An additional hypothesis is made in order to predict a structure's final failure in static fatigue: we assume that the mechanical fields do not vary significantly after the application of the mechanical load, and therefore consider them constant in time. Consequently, the mechanical damage model presented in Section 17.2 is used to compute the mechanical response of the structure to the mechanical load, and the environmental damage model presented in Section 17.3 is used as a post-treatment of the mechanical data to compute the environmental degradation of the material with time. Final failure, and then lifetime, are predicted thanks to a simple criterion defined now.

#### 17.4.2.1 Failure Criterion
It is well known that local failure criteria cannot predict the final failure of materials properly [83]. Consequently, we use a nonlocal failure criterion, based on Reference 84, and recently applied to polymer matrix composites [85] as well as ceramic matrix composites [86].

The criterion is very simple both in concept and in practice (see Figure 17.27). It consists, for each tow direction, of averaging the stress field over an area $\Omega$ defined by a length $l_0$ considered to be an intrinsic parameter of the material. (In practice, we use the size of a weave pattern, that is, 2 mm.) For instance, for the longitudinal tow, the failure criterion is

$$ \exists \underline{x} / \tilde{\sigma}_1(\underline{x}) \geq \sigma_1^R(\underline{x}) \quad \Rightarrow \quad \text{failure} \qquad (17.179) $$

$$ \text{with } \tilde{\sigma}_1(\underline{x}) = \frac{1}{|\Omega|} \int_{\Omega} \sigma_{11}^{+}(\underline{x}') \, d\underline{x}' \quad \forall \underline{x} $$

The formulation for the transversal tow is similar.

This failure criterion can be rewritten differently based on the normalized residual strength (or strength reserve factor) field [87]. This quantity provides the best illustration of the relation between the stress field (imposed by the

<figure><img src="image_29.png" /><figcaption>**FIGURE 17.27**
Nonlocal failure criterion for a structure with a high stress gradient.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>FIGURE 17.28</figcaption></figure>

The geometry of the dog-bone specimen subjected to tension in the x-direction.

boundary conditions) and the strength field (which evolves over time because of environmental damage) during the failure process, and is defined by

$$ \tilde{\sigma}_1^R(x) = \frac{\sigma_1^R(x)}{\tilde{\sigma}_1(x)} \quad \forall x \qquad (17.180) $$

Then, the failure criterion becomes

$$ \exists x / \tilde{\sigma}_1^R(x) \le 1 \Rightarrow \text{failure} \qquad (17.181) $$

**17.4.2.2 A First Illustration** The first example is a dog-bone specimen made of SAFRAN–Snecma Propulsion Solide’s CERASEP A410 self-healing CMC (Figure 17.28). The FE mesh was generated using GMSH [88]. The specimen was subjected to a static fatigue loading at high temperature: all the loading conditions, that is, the applied stress and temperature and the other environmental conditions (20 kPa partial oxygen pressure; 1.8 kPa partial water pressure) remained constant over time.

We used an identification of the model for this material provided by References 11, 14, and 15. While some parameters (e.g., the reaction/diffusion coefficients) can be found in the literature, most were identified by fitting the mechanical response of macroscopic specimens (e.g., initial behavior, damage evolution laws) or by analyzing the morphology of the material on the fiber’s scale (e.g., the characteristic dimensions of the matrix layers). Still, a few parameters (e.g., the anisotropy of the cracks) remained unidentified because of the lack of experimental data and had to be chosen arbitrarily [14, 15].

We will present first the damage induced by the mechanical loading (i.e., the stiffness reduction), followed by the damage induced by the environment (i.e., the strength reduction).

**Damage fields predictions** As mentioned previously, the key property of the model is that it distinguishes among the contributions associated with each crack network, which, thus, can be analyzed separately. For instance, for a 50 MPa applied load, inter-yarn damage alone is in effect, whereas for a 250 MPa applied load (Figure 17.29a), both inter-yarn matrix damage and intra-yarn matrix damage are in effect (Figures 17.29b and 17.29c): intra-yarn damage takes over after inter-yarn damage has reached saturation.

The mechanical part also provides an intra-yarn transverse crack opening indicator field, which is the key input to the physicochemical model of strength reduction through oxidation (see Section 17.3). Now, let us illustrate the capability of our strategy to simulate the complex competition between mechanical loading and environmental loading through the analysis of two related test cases: a hot-grip experiment (i.e., with a homogeneous temperature field) and a cold-grip experiment (i.e., with a temperature gradient over the specimen).

**Lifetime predictions, case of homogeneous temperature field** The use of hot-grip testing machines greatly simplifies the study of the lifetime of self-healing CMCs because these machines allow the analysis of the influence of the stress and temperature fields on the lifetime to be fully uncoupled. Therefore, let us first consider the use of our model to analyze hot-grip experiments.

Figure 17.30 shows the strength degradation in the specimen after 278 hours under 200 MPa loading at two uniform temperatures (500° and 700°). In this case of a homogeneous temperature field, the result of interaction with the damage field (which is nearly homogeneous too) is that the strength decreases mainly in the middle of the specimen, where the stress is higher. For instance, from the initial value of 340 MPa, the strength reduced to about 120 MPa at
<figure><img src="image_31.png" /><figcaption>FIGURE 17.29</figcaption></figure>

The dog-bone specimen subjected to a 250 MPa traction: (a) longitudinal stress, (b) inter-yarn matrix damage, (c) intra-yarn matrix damage.

500° (Figure 17.30a), and to only about 210 MPa at 700° (Figure 17.30b). Figure 17.31 shows the same type of result, this time for the strength reserve factor field. This shows clearly that at 500° the specimen breaks in the zone of interest whereas at 700° it does not break at all.

This shows that the model is capable of illustrating and quantifying classical types of behavior: (i) an increase in the applied stress decreases the lifetime of the specimen drastically (see Figure 17.32, which indicates a very good match with experimental data); (ii) a modification in the applied temperature enables one to define the domain of applicability of the healing process ([11, 15]). Let us note that the wide scatter in the lifetime data of Figure 17.32 was evidenced only in a recent study [89] and will be analyzed within our modeling framework in a forthcoming paper.

*Lifetime predictions, case of temperature gradient* Hot-grip testing machines are uncommon and most experiments are carried out using cold-grip machines. This makes the analysis complicated because neither the stress field nor the temperature field, both of which have a significant impact on the lifetime, is uniform [89, 93]. This led us to apply our model to the analysis of cold-grip experiments in order to illustrate how it can help improve their design and understanding. The temperature field considered was obtained by References 91 and 92 and is shown in Figure 17.33.

<figure><img src="image_32.png" /><figcaption>FIGURE 17.30</figcaption></figure>

Residual strength after mechanical (200 MPa) and environmental (278 hours, without temperature gradient) loading.
<figure><img src="image_33.png" /><figcaption>FIGURE 17.31</figcaption></figure>

Strength reserve factor after mechanical (200 MPa) and environmental (278 hours, without temperature gradient) loading.

<figure><img src="image_34.png" /><figcaption>FIGURE 17.32</figcaption></figure>

Comparison of the lifetime predictions between model and experiment: stress dependence at 500° (experimental data from References 89–92).

<figure><img src="image_35.png" /><figcaption>FIGURE 17.33</figcaption></figure>

The temperature gradient in the dog-bone specimen in a cold-grip testing machine at two temperatures [89].
<figure><img src="image_36.png" /><figcaption>FIGURE 17.34</figcaption></figure>

With a temperature gradient: the residual strength at two temperatures (200 MPa loading for 278 hours).

Figures 17.34 and 17.35 show the residual strength and the strength reserve factor in the specimen after 278 hours under 200 MPa loading at two nonuniform temperatures (500° max and 700° max). For a 500° prescribed temperature in the center of the specimen, the predicted strengths (Figures 17.34a and 17.35a) are very close to those with a uniform temperature (Figures 17.30a and 17.31a). Indeed, in this temperature range the healing process has almost no influence. Conversely, for a 700° prescribed temperature in the center of the specimen, the predicted strengths (Figures 17.34b and 17.35b) are markedly different from those with a uniform temperature (Figures 17.30b and 17.31b): the strength degradation patterns can be very complex. Indeed, due to the fact that the effectiveness of the oxide tip varies significantly, and especially non-monotonically, as a function of the temperature [11], the strength degradation can be much greater even though the mechanical loading is lesser. It was even possible to predict the fracture zone in the grips (Figure 17.35b), a phenomenon that was observed experimentally and that led to major experimental problems [89, 93]. Thus, our model could be used as an experiment design tool in order to ensure that fracture would take place in the zone of interest of the specimen. Besides, this case illustrates the need to consider the whole structure for the lifetime study, if only because the fracture zone cannot always be determined *a priori*.

#### 17.4.2.3 Application to Damage Tolerance Analysis
In this section, we illustrate the application of our strategy to the damage tolerance analysis of structures, which is one of the main concerns associated with industrial developments [18, 94]. The objective is to characterize the influence of damaging impacts on the residual lifetime of structures.

According to experimental observations, classical quasi-static indentations generate cone-shaped damaged zones, which can be approximated by cylindrical holes [94]. Let us consider the structural example of an open-hole specimen made of the same material and subjected to the same type of loading as the dog-bone specimen of the previous section. Thus, this test case involves multiaxial stress fields and high stress gradients, especially near the hole. Figure 17.36 shows the geometry of the open-hole plate (with the FE mesh generated once again using GMSH [88]) subjected to tension loading in the x-direction.

**Damage fields predictions** The model predicts that with an applied stress equal to 25 MPa, inter-yarn matrix cracks initiate around the hole, whereas with an applied stress equal to 55 MPa, intra-yarn matrix cracks initiate, also around

<figure><img src="image_37.png" /><figcaption>FIGURE 17.35</figcaption></figure>

With a temperature gradient: the strength reserve factor at two temperatures (200 MPa loading for 278 hours).
<figure><img src="image_38.png" /><figcaption>**FIGURE 17.36**</figcaption></figure>

The geometry of the open-hole specimen subjected to tension loading in the x-direction.

the hole. Moreover, with an applied stress equal to 250 MPa (Figure 17.37), the inter-yarn damage becomes saturated within the structure and intra-yarn matrix damage appears near the hole (Figures 17.37b and 17.37c). These predictions, which are consistent with the available experimental data, show that there is no unstable propagation of the damage leading to purely mechanical fracture of the specimen. Therefore, it is now necessary to address the damage tolerance problem associated with delayed fracture, which conditions the structure's lifetime.

<figure><img src="image_39.png" /><figcaption>**FIGURE 17.37**</figcaption></figure>

The open-hole specimen under tension loading: (a) longitudinal stress, (b) inter-yarn matrix damage, (c) intra-yarn matrix damage.
<figure><img src="image_40.png" /><figcaption>**FIGURE 17.38**</figcaption></figure>

The open-hole specimen subjected to tension loading: (a) residual strength, (b) strength reserve factor (200 MPa loading at 500° for $t_0$ = 110 hours).

**Lifetime predictions and damage tolerance analysis** Figure 17.38 shows typical distributions of the indicator. In the uniform temperature case, the model predicts that at any temperature the strength degradation takes place mainly around the hole, where final fracture initiates.

In order to study the effect of the hole on the structure's lifetime, similar simulations were performed on plates with different hole diameters. The influence of the hole was measured by comparing the lifetimes of open-hole and dog-bone specimens under the same effective stress (i.e., the same stress in the midplane of the plates). The results (Figure 17.39) indicate that the influence of the hole is greater at low macroscopic stress levels. Moreover, it appears that even under identical effective stresses large defects have more of an impact on the lifetime than small defects. Indeed, the damage zone is larger, which induces greater chemical degradation. It is important to note that this phenomenon, which is specific to lifetime analysis and does not exist in the case of strength analysis, is inherently part of our strategy, which illustrates its suitability for damage tolerance analysis in an industrial framework.

Also, in Figure 17.39, one should note that the ratio of the open-hole stress (for a given hole) over the dog-bone stress (for a given lifetime) is almost constant. This constant is a stress concentration factor, not in the classical sense, but dedicated to delayed fracture. This stress concentration factor could be used by engineers as a very simple design tool.

<figure><img src="image_41.png" /><figcaption>**FIGURE 17.39**</figcaption></figure>

Damage tolerance analysis: the lifetime as a function of the effective stress at 500° for several hole diameters.
### 17.4.3 Control of the Damage Localization and Static Failure Predictions

We now demonstrate the ability of the mechanical damage model presented in Section 17.2 to predict the initiation and propagation of the macroscopic crack that will induce the structure's final failure. However, it is well known that the calculation of softening damage laws such as the ones of Equations (17.41) and (17.46) over a domain goes through a critical point where a loss of ellipticity, that is, a loss of uniqueness of the solution [20, 21], occurs and leads to pathological mesh dependencies [27, 95]. Several remedies exist for damage models, including nonlocal formulations [20, 95] and delay-effect approaches [27, 95]. The latter was chosen for our work because it is a local method, and then is easy to implement in classical FE codes.

#### 17.4.3.1 The Delay-Effect Method
The delay-effect method for controlling the localization of damage consists of replacing the evolution of the scalar damage variable $\alpha$ (defined in Equation (17.41) for matrix microcracking and in Equation (17.46) for tow microcracking, and calibrated in Section 17.2.4.5) by an evolution with a bounded rate:

$$ \dot{\alpha} = \frac{1}{\tau_c} (1 - e^{-\langle \alpha^{\text{static}} - \alpha \rangle_+}) \quad (17.182) $$

where $\alpha^{\text{static}}$ is the static function calibrated in Section 17.2.4.5 and $\tau_c$ a fictitious time parameter. Equation (17.182) has a strong physical meaning in dynamic [27] (cracks can-not appear and propagate instantaneously), but is used here in quasi-static as a pure numerical method. Thus, the fictitious time parameter is not related to any physical process; it must simply be chosen in order to avoid uncontrolled localization. The numerical resolution of such a simple nonlinear, but scalar, equation is straightforward and will not be discussed here.

The behavior law is modified with the introduction of this viscous-like parameter. Figure 17.40 shows the stress-strain curves obtained with different values of $\tau_c$.

#### 17.4.3.2 Illustrative Examples

**CT specimen** The effectiveness of the delay-effect method in controlling the localization of damage was evaluated using a simple CT-like test, that is, a pre-cracked specimen subjected to tension (see Figure 17.41). The response and the damage fields for different meshes (generated with GMSH [88]) are shown in Figure 17.42: one can clearly observe that the solution is objective, that is, mesh independent, in terms of both the mechanical response and the damage fields.

**Plate with an open hole** Now let us consider a first structural example that consists of a plate with an open hole (see Figure 17.43). The model was the same as for the CT specimen, except for $\tau_c = 110^{-2} \text{s}$. Figure 17.44 shows the response of the model using different meshes (again generated with GMSH [88]). Once again, one can clearly observe that the method enables an objective prediction of the final failure of the structure.

<figure><img src="image_42.png" /><figcaption>**FIGURE 17.40**</figcaption></figure>

The response of the model under unidirectional tension loading using different characteristic times for the delay-effect law.
<figure><img src="image_43.png" /><figcaption>**FIGURE 17.41**</figcaption></figure>

The problem used for the evaluation of localization control based on delay-effect damage ($\epsilon_1^g$ denotes the applied strain).

<figure><img src="image_44.png" /><figcaption>**FIGURE 17.42**</figcaption></figure>

The response and the damage field for the CT specimen at $t = 0.75$ s (dark grey: no damage; light grey: full damage): objectivity of the mesh and prediction of the crack's propagation ($\epsilon_1^g = 0.2\%/s$), DOFs: Degrees of freedom.
<figure><img src="image_45.png" /><figcaption>FIGURE 17.43</figcaption></figure>

Structural case: plate with an open hole ($\epsilon_1^g$ denotes the applied strain).

<figure><img src="image_46.png" /><figcaption>FIGURE 17.44</figcaption></figure>

The response and the damage field for the open hole specimen at $t = 1$ s (dark grey: undamaged; light grey: full damage): objectivity of the mesh and prediction of the final failure ($\epsilon_1^g = 0.1\%$), DOFs: Degrees of freedom.

REFERENCES

1. F. Lamouroux, S. Bertrand, R. Pailler, R. Naslain, and M. Cataldi. Oxidation-resistant carbon-fiber-reinforced ceramic-matrix composites. *Comp. Sci. Technol.*, 59(7), 1073–1085 (1999).

2. P. Forio and J. Lamon. Fatigue behavior at high temperatures in air of a 2 D SiC/Si-B-C composite with a self-healing multilayered matrix. *Ceram. Trans. (USA)*, 128, 127–141 (2001).

3. J. Aveston and A. Kelly. Theory of multiple fracture of fibrous composites. *J. Mater. Sci.*, 8(3), 352–362 (1973).

4. A. Evans, J. Domergue, and E. Vagaggini. Methodology for relating the tensile constitutive behavior of ceramic-matrix composites to constituent properties. *J. Am. Ceram. Soc.*, 77(6), 1425–1435 (1994).
5. D. Rouby and N. Louet. The frictional interface: a tribological approach of thermal misfit, surface roughness and sliding velocity effects. *Compos. Part A: Appl. Sci. Manufact.*, 33(10), 1453–1459 (2002).

6. P. Ladevèze. On an anisotropic damage Theory [in French]. Internal Report 34 (LMT-Cachan, March 1983).

7. P. Ladevèze. An anisotropic damage theory with unilateral effects: applications to laminate and three- and four-dimensional composites, in *Continuum Damage Mechanics of Materials and Structures*, edited by O. Allix and F. Hild (Elsevier, 2002), pp. 205–233.

8. A. Gasser, P. Ladeveze, and M. Poss. Damage mechanisms of a woven SiC/SiC composite: modelling and identification. *Compos. Sci. Technol.*, 56(7), 779–784 (1996).

9. P. Ladevèze and S. Letombe. Recent advances on an anisotropic damage theory including unilateral effects. Proceedings of the International Bimestre on Damage Mechanics. Symposium on Continuous Damage and Fracture, Cachan, France, October 23–27, 2000.

10. P. Ladevèze, S. Letombe, and C. Cluzel. A CMC damage model based on micro and macromechanics for high-temperature and complex loading. Proceedings of the 4th International Conference on High Temperature Ceramic Matrix Composites [HTCMC4], Munich, Germany, October 1–3, 2001.

11. C. Cluzel, E. Baranger, P. Ladevèze, and A. Mouret. Mechanical behaviour and lifetime modelling of self-healing ceramic-matrix composites subjected to thermomechanical loading in air. *Compos. Part A: Appl. Sci. Manufact.*, 40, 976–984, (2009).

12. S. Letombe, C. Cluzel, and P. Ladevèze. A macroscopic model coupling oxidation and damage for CMCs. Proceedings of JNC13 Conference, Strasbourg, France, March 12–14, 2003, pp. 713–722.

13. O. Penas, P. Reynaud, D. Rouby, and G. Fantozzi. Self-healing SiC/SiC composite behaviour under high-temperature cyclic fatigue in air, in *High Temperature Ceramic Matrix Composites*, edited by W. Krenkel, R. Naslain, and H. Schneider. (Wiley-VCH, 2001), pp. 480–485.

14. E. Baranger, C. Cluzel, P. Ladevèze, and A. Mouret. Prediction of the lifetime of self-healing ceramic matrix composites: I—macroscopic modeling of cracking. Proceedings of the 15th National Conference on Composite Materials (JNC15) [in French], Marseille, France, June 6–8, 2007.

15. C. Cluzel, E. Baranger, P. Ladevèze, and A. Mouret. Prediction of the lifetime of self-healing ceramic matrix composites: II—analysis of cracking and oxidation mechanisms. Proceedings of the 15th National Conference on Composite Materials (JNC15) [in French], pp. 885–894, 2007.

16. M. Genet, L. Marcin, E. Baranger, C. Cluzel, P. Ladevèze, and A. Mouret. Computational prediction of the lifetime of self-healing CMC structures. *Compos. Part A: Appl. Sci. Manufact.*, 43, 294–303 (2012).

17. B. Cox and Q. Yang. In quest of virtual tests for structural composites. *Science*, 314(5802), 1102–1107 (2006).

18. M. Bourgeon. Thermostructural materials in aerospace industry: applications and standardization. Proceedings of the 7th International Conference on High Temperature Ceramic Matrix Composites [HTCMC7], Bayreuth, Germany, September 20–22, 2010.

19. J. Lemaître, J.-L. Chaboche, R. Desmorat, and A. Benallal. *Solid Materials Mechanics*, 3rd edition [in French]. (Dunod, 2009).

20. Z. Bazant and G. Pijaudier-Cabot. Nonlocal continuum damage, localization instability and convergence. *J. Appl. Mech.* [Transactions of the ASME], 55(2), 287–93 (1988).

21. A. Needleman. Material rate dependance and mesh sensitivity in localization problems. *Comput. Meth. Appl. Mech. Eng.*, 67, 69–85 (1988).

22. D. Lasry and T. Belytschko. Localization limiters in transient problems. *Int. J. Solids Struct.*, 24, 581–597 (1988).

23. L. Sluys and R. De Borst. Wave propagation and localization in a rate-dependent cracked medium—model formulation and one-dimensional examples. *Int. J. Solids Struct.*, 29, 2945–2958 (1992).

24. P. Ladevèze. On a damage mechanics approach, in *Mechanics and Mechanisms of Damage in Composites and Multi-Materials*, edited by D. Baptiste (MEP, SaintÉtienne, France, 1989), pp. 119–141.

25. P. Ladevèze. Modeling and Simulation of the Mechanical Behavior of CMCs. *High Temp. Ceram. Matrix Compos.*, 47, 53–63 (1995).

26. O. Allix and J.-F. Deü. Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading. *Eng. Trans.*, 45(1), 29–46 (1997).

27. P. Ladevèze, O. Allix, J.-F. Deü, and D. Lévêque. A mesomodel for localisation and damage computation in laminates. *Comput. Meth. Appl. Mech. Eng.*, 183, 105–122 (2000).

28. P. Ladevèze. Inelastic Strains and Damage, in *Damage Mechanics of Composite Materials*, edited by R. Talreja (Elsevier, Amsterdam, The Netherlands, 1994), pp. 117–138.

29. J. Lemaître. How to use damage mechanics. *Nucl Eng. Des.*, 80, 233–245, (1984).

30. P. Ladevèze and J. Lemaître. Damage effective stress in quasi unilateral conditions. IUTAM Congress, Lyngby, Denmark, August 19–25, 1984.
31. P. Ladevèze. On the damage mechanics of composites [in French]. Comptes rendus des Cinquièmes Journées nationales sur les composites, Paris, France, September 9-11, 1986, edited by C. Bathias and D. Menkès. pp. 667-683.

32. M. Ortiz. A constitutive theory for the inelastic behavior of concrete. *Mech. Mater.*, 4, 67-93, (1985).

33. J. Simo and J. Ju. Strain-based and stress-based continuum damage models. 1. Formulation. *Int. J. Solids Struct.*, 23(7), 821–840 (1987).

34. J. Simo and J. Ju. Strain-based and stress-based continuum damage models. 2. Computational Aspects. *Int. J. Solids Struct.*, 23(7), 841–869 (1987).

35. P. Ladevèze, A. Gasser, and O. Allix. Damage mechanisms modelling for ceramic composites. *J. Eng. Mater. Technol.*, 116, 331–336 (1994).

36. J. Dubé, G. Pijaudier-Cabot, and C. Laborderie. Rate dependent damage model for concrete in dynamics. *J. Eng. Mech.*, 122(10), 939–947 (1996).

37. R. Desmorat, F. Gatuingt, and F. Ragueneau. Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials. *Eng. Fract. Mech.*, 74, 1539–1560 (2007).

38. J. Lemaître, R. Desmorat, and M. Sauzay. Anisotropic damage law of evolution. *Eur. J. Mech.*, 19(2), 187–208 (2000).

39. J.-L. Chaboche. Damage induced anisotropy: on the difficulties associated with the active/passive unilateral condition. *Int. J. Damage Mech.*, 1, 148–171 (1992).

40. J.-L. Chaboche. Development of continuum damage mechanics for elastic solids sustaining anisotropic and unilateral damage. *Int. J. Damage Mech.*, 2, 311–329 (1993).

41. D. Halm and A. Dragon. A model of anisotropic damage by mesocrack growth; unilateral effect. *Int. J. Damage Mech.*, 5, 384–402 (1996).

42. F. Cormery and H. Welemane. A critical review of some damage models with unilateral effect. *Mech. Res. Commun.*, 29, 391–395 (2002).

43. F. Cormery and H. Welemane. A stress-based macroscopic approach for microcracks unilateral effect. *Comput. Mater. Sci.*, 47, 727–738, (2010).

44. P. Ladevèze. On an anisotropic damage theory. in *CNRS International Colloquium on Failure Criteria of Structured Media*, Villard de Lans, France, June 21–24, 1983, edited by J. Boehler (Balkema, Rotterdam, 1993), pp. 335–365.

45. Q.-C. C. He and A. Curnier. A more fundamental approach to damaged elastic stress-strain relations. *Int. J. Solids Struct.*, 32, 1433–1457 (1995).

46. C. Rospars, J.-L. Chermant, and P. Ladevèze. On a first creep model for a 2D SiC-SiC composite. *Mater. Sci. Eng.: A*, 250, 264–269 (1998).

47. L. Guillaumat and J. Lamon. Multi-fissuration de composites SiC/SiC. *Revue des composites et des matériaux avancés*, 3, 159–171 (1993).

48. E. Baranger, C. Cluzel, P. Ladevèze, and A. Mouret. Identification and validation of a multi-physic macro model for the lifetime prediction of self-healing ceramic matrix composites. Proceedings of the 13th European Conference on Composite Materials [ECCM 13], Stockholm, Sweden, June 2–5, 2008.

49. M. Moevus, M. Mécanismes d'endommagement, émission acoustique et durées de vie en fatigue statique du composite SiC/[Si-BC] aux températures intermédiaires (< 800 C). PhD thesis. INSA-Lyon (2007).

50. H. Mei and L. Cheng. Stress-dependence and time-dependence of the post-fatigue tensile behavior of carbon fiber reinforced SiC matrix composites. *Compos. Sci. Technol.*, 71, 1404–1409 (2011).

51. D. Rouby and P. Reynaud. Fatigue behaviour related to interface modification during load cycling in ceramic-matrix fibre composites. *Compos. Sci. Technol.*, 48(1-4), 109–118 (1993).

52. P. Reynaud, D. Rouby, and G. Fantozzi. Effects of temperature and of oxidation on the interfacial shear stress between fibres and matrix in ceramic-matrix composites. *Acta Mater.*, 46(7), 2461–2469 (1998).

53. N. Lissart and J. Lamon. Damage and failure in ceramic matrix minicomposites: experimental study and model. *Acta mater.*, 45(3), 1025–1044 (1997).

54. F. Rebillat, X. Martin, E. Garitte, and A. Guette. Overview on the self-sealing process in the SiCf/[Si, C, B]m composites under wet atmosphere at high temperature, in *Design, Development, and Applications of Engineering Ceramics and Composites*, edited by D. Singh, D. Zhu, and Y. Zhou. Vol. 215: Ceramic Transactions, edited by M. Singh (John Wiley & Sons, 2010), pp. 151–166.

55. X. Martin, F. Rebillat, and A. Guette. Oxidation behavior of a multilayered (Si-BC) ceramic in a complex atmosphere N2/O2/H2O. *Proc. 4th Int. Symp. High Temp. Corros. Mater. Chem.*, 339 (2003).

56. X. Martin, Oxidation/corrosion de matériaux composites (SiCf/SiBCm) à matrice auto-cicatrisante. PhD thesis. Université Bordeaux I (2003).

57. E. Garitte, Etude de l'oxydation/corrosion des composites céramiques. PhD thesis. Université Bordeaux I (2007).
58. L. Filipuzzi and R. Naslain. Oxidation mechanisms and kinetics of 1D-SiC/C/SiC composite materials: II. Modeling. *J. Am. Ceram. Soc.*, 77(2), 467-480 (1994).

59. F. Norton. Permeation of gaseous oxygen through vitreous silica. *Lett. Nat.*, 191, 701 (1961). Doi:10.1038/191701a0.

60. B. Deal and A. Grove. General relationship for the thermal oxidation of silicon. *J. Appl. Phys.*, 36(12), 3770-3778 (1965).

61. J. Schlichting and S. Neumann. GeO₂/SiO₂-glasses from gels to increase the oxidation resistance of porous silicon containing ceramics. *J. Non-Cryst. Solids*, 48(1), 185-194 (1982).

62. J. Schlichting. Oxygen transport through glass layers formed by a gel process. *J. Non-Cryst. Solids*, 63(1-2), 173-181 (1984).

63. J. Viricelle, P. Goursat, and D. Bahloul-Hourlier. Oxidation behaviour of a multi-layered ceramic-matrix composite (SiC)f/C/(SiBC)m. *Compos. Sci. Technol.*, 61(4), 607-614 (2001).

64. R. Reid, J. Prausnitz, and B. Poling. *The Properties of Gases and Liquids* (McGraw Hill Book Co., New York, 1987).

65. F. Rebillat, X. Martin, and A. Guette. Kinetic oxidation laws of boron carbide in dry and wet environments. *Proceedings of HTCMC*, 5, 321-326 (2004).

66. S. Letombe, Modélisation du couplage oxydation/endommagement des Composites à Matrice Céramique Autocicatrisante. PhD thesis. ENS-Cachan (2005).

67. A. G. Evans, F. W. Zok, R. M. McMeeking, and Z. Z. Du. Models of high-temperature, environmentally assisted embrittlement in ceramic-matrix composites. *J. Am. Ceram. Soc.*, 79, 2345-2352 (1996).

68. P. Forio, F. Lavaire, and J. Lamon. Delayed failure at intermediate temperatures (600 degrees-700 degrees C) in air in silicon carbide multifilament tows. *J. Am. Ceram. Soc.*, 87, 888-893 (2004).

69. W. Gauthier and J. Lamon. Delayed failure of Hi-Nicalon and Hi-Nicalon S multifilament tows and single filaments at intermediate temperatures (500 degrees-800 degrees C). *J. Am. Ceram. Soc.*, 92, 702-709. (2009).

70. B. R. Lawn. Diffusion-controlled subcritical crack growth in the presence of a dilute gas environment. *Mater. Sci. Eng.*, 13, 277-283 (1974).

71. S. M. Wiederhorn, E. R. Fuller, and T. Thomson. Micromechanics of crack growth in ceramics and glasses in corrosive environments. *Met. Sci.*, 14, 450-458 (1980).

72. J. B. Wachtman. *Mechanical Properties of Ceramics*. (John Wiley & Sons, Inc., 1996).

73. F. Hild, D. Marquis, O. Kadouch, and J. P. Lambelin. Analysis of the failure of ceramics due to subcritical crack growth. *J. Eng. Mater. Technol.*, 118, 343-348 (1996).

74. P. Ladevèze and M. Genet. A new approach to the subcritical cracking of ceramic fibers. *Compos. Sci. Technol.*, 70, 1575-1583, (2010).

75. A. Laforêt and J. Lamon. Static fatigue of multifilament tows at high temperatures above 900°C. *Proceedings of the 13th European Conference on Composite Materials* [ECCM 13], Stockholm, Sweden, June 2-5, 2008.

76. P. Ladevèze and M. Genet. A new approach to the subcritical cracking of ceramic fibers. *Compos. Sci. Technol.*, 70(11), 1575-1583 (2010).

77. W. Gauthier, F. Pailler, J. Lamon, and R. Pailler. Oxidation of silicon carbide fibers during static fatigue in air at intermediate temperatures. *J. Am. Ceram. Soc.*, 92(9), 2067-2073 (2009).

78. F. Rebillat, E. Garitte, and A. Guette. Quantification of higher SiC fiber oxidation rates in presence of B₂O₃ under air, in *Design, Development, and Applications of Engineering Ceramics and Composites*, edited by D. Singh, D. Zhu, and Y. Zhou. Vol. 215: Ceramic Transactions, edited by M. Singh (John Wiley & Sons, 2010), pp. 135-149.

79. S. Letombe. Modelling of the damage/oxydation coupling in ceramic matrix composites [in French]. PhD thesis. ENS-Cachan, France (2005).

80. M. Genet, L. Marcin, and P. Ladevèze. On structural computations until fracture based on an anisotropic and unilateral damage theory. *Int. J. Damage Mech.*, 23, 483-506 (2014). Doi:10.1177/1056789513500295.

81. P. Ladevèze. *Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation*. Mechanical Engineering Series (Springer, 1999).

82. P. Ladevèze, J.-C. Passieux, and D. Néron. The LATIN multiscale computational method and the proper generalized decomposition. *Comput. Meth. Appl. Mech. Eng.*, 199(21-22), 1287-1296 (2010).

83. R. J. Whitney and J. M. Nuismer. Stress fracture criteria for laminated composite laminates containing stress concentrations. *J. Compos. Mater.*, 8, 253-265 (1974).

84. J. M. Nuismer and R. J. Whitney. Uniaxial failure of composite failure laminates containing stress concentrations, in *Fracture Mechanics of Composites* (ASTM STP 593, 1975), pp. 117-142.

85. S. Miot, C. Hochard, and N. Lahellec. A non-local criterion for modelling unbalanced woven ply laminates with stress concentrations. *Compos. Struct.*, 92, 1574-1580 (2010).
86. S. Flores, A. G. Evans, F. W. Zok, M. Genet, B. N. Cox, D. B. Marshall, O. Sudre, and Q. D. Yang. Treating matrix nonlinearity in the binary model formulation for 3D ceramic composite structures. *Compos. Part A: Appl. Sci. Manufact.*, 41, 222–229 (2010).

87. Department of Defence. *Composite Materials Handbook*, vol. 5. (Ceramic Matrix Composites, Department of Defence, USA, 2002).

88. C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Numer. Meth. Eng.*, 79, 1309–1331 (2009).

89. O. Loseille and J. Lamon. Prediction of lifetime in static fatigue at high temperatures for ceramic matrix composites. *Adv. Mater. Res.*, 112, 129–140, (2010).

90. O. Penas, Study of SiC/SiBC multi-layered matrix composites in cyclic fatigue at high temperature under air [in French]. Phd thesis. INSA-Lyon, France (2002).

91. M. Moevus, D. Rouby, N. Godin, M. R'Mili, P. Reynaud, G. Fantozzi, and G. Farizy. Analysis of damage mechanisms and associated acoustic emission in two SiC/[Si-B-C] composites exhibiting different tensile behaviours. Part I: damage patterns and acoustic emission activity. *Compos. Sci. Technol.*, 68, 1250–1257 (2008).

92. M. Moevus, N. Godin, M. R'Mili, D. Rouby, P. Reynaud, G. Fantozzi, and G. Farizy. Analysis of damage mechanisms and associated acoustic emission in two SiCf/[Si-B-C] composites exhibiting different tensile behaviours. Part II: unsupervised acoustic emission data clustering. *Compos. Sci. Technol.*, 68, 1258–1265 (2008).

93. S. Momon, M. Moevus, N. Godin, M. R'Mili, P. Reynaud, G. Fantozzi, and G. Fayolle. Acoustic emission and lifetime prediction during static fatigue tests on ceramic-matrix-composite at high temperature under air. *Compos. Part A: Appl. Sci. Manufact.*, 41, 913–918 (2010).

94. V. Herb, G. Couégnat, and E. Martin. Damage assessment of thin SiC/SiC composite plates subjected to quasi-static indentation loading. *Compos. Part A: Appl. Sci. Manufact.*, 41, 1677–1685 (2010).

95. R. Desmorat, M. Chambart, F. Gatuingt, and D. Guilbaud. Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading. *Eng. Fract. Mech.*, 77, 2294–2315 (2010).