# Finite element modelling of the oxidation gradients of epoxy-diamine matrices below and above their glass transition temperature

Juan Pablo Márquez Costa, Xavier Colin

► To cite this version:

Juan Pablo Márquez Costa, Xavier Colin. Finite element modelling of the oxidation gradients of epoxy-diamine matrices below and above their glass transition temperature. Polymer Degradation and Stability, 2025, 234, pp.111194. 10.1016/j.polymdegradstab.2025.111194. hal-04973228

HAL Id: hal-04973228

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

Submitted on 2 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Finite element modelling of the oxidation gradients of epoxy-diamine matrices below and above their glass transition temperature

Juan Pablo Márquez Costa and Xavier Colin *

PIMM, Arts et Métiers Institute of Technology, CNRS, CNAM, HESAM, 151 Boulevard de l'Hôpital, 75013 Paris, France

* Corresponding authors: xavier.colin@ensam.eu; Tel.: +33 1 44 24 61 47

## Abstract

This article outlines a numerical calculation code developed under finite elements to simulate as accurately as possible the thermal aging of epoxy-diamine (EPO-DA) matrices below and above their glass transition temperature ($T_g$). The code is based on a kinetic model that couples, in a single balance equation, the diffusion of oxygen ($O_2$) and its consumption by the chemical reaction with the epoxy network. The reliability of the oxidation rate expression was carefully checked for six EPO-DA matrices over wide ranges of temperature and $O_2$ partial pressure in previous publications. The code gives access to the oxidation gradients within the sample thickness whose validity is successfully checked through a literature compilation (completed by our own recent results) of more than 40 thicknesses of oxidized layer (TOL), measured with different physico-chemical techniques on various EPO-DA matrices thermally aged in air below and above their $T_g$. In the glassy state, where the coefficient of $O_2$ diffusion ($D_{O_2}$) obeys a well-established Arrhenius' law depending on the chemical structure of the diamine hardener (either aliphatic or aromatic), the estimates of TOL are in nice agreement with the experimental data. In the rubbery state, in contrast, the calculations largely underestimate the experimental data for the stiffest EPO-DA matrices if a sharp increase in $D_{O_2}$ is not also considered when passing above $T_g$. This positive jump in $D_{O_2}$ places all EPO-DA matrices at the same level of
molecular mobility in the rubbery state, thus completely erasing the stiffening effect of aromatic
rings observed in the glassy state.

Keywords: Epoxy-diamine matrices, Thermal oxidation, Oxidation gradients, Thickness of oxidized layer, Kinetic modelling, Diffusion/reaction coupling, Molecular mobility.
# 1. Introduction

Carbon/epoxy composites are often considered for applications as structural parts in the aeronautic field due to their high specific mechanical properties. However, there are potential risks to consider during the design and certification stages, particularly when these materials are exposed to high temperatures close to their glass transition temperature ($T_g$). Among these risks, those related to the thermal oxidative degradation of the epoxy matrix are particularly critical [1-7]. Although restricted to the superficial layers of the part due to its kinetic control by oxygen ($O_2$) diffusion [8, 9], thermal oxidation leads simultaneously to the development of tensile stress gradients and a catastrophic drop in fracture properties. These two combined effects are responsible for spontaneous cracking [10]. The new matrix surfaces freshly exposed to the environment can then oxidize, thus allowing this scenario to repeat and the oxidation-induced cracks to reach the core of the part [10]. Thenceforth, the failure of the part can occur prematurely.

In order to better address this complex issue of durability, successively involving oxidation and initiation then propagation of cracks in a material with gradients of physico-chemical and mechanical properties, but also involving a coupling between oxidation and cracking, it is first necessary to develop new scientific concepts and calculation tools. The first tool to develop is a diffusion/reaction model giving accurately access to the oxidation gradients within the part before cracking. As already shown by numerous research teams for many unfilled hydrocarbon polymers since the 1950s [8, 9, 11-20], such a model can be formulated as follows in the case of the unidirectional diffusion of $O_2$ along the z axis (e.g. through the thickness of the part):

$$ \frac{\partial C(z,t)}{\partial t} = D_{O2} \frac{\partial^2 C(z,t)}{\partial z^2} - r(C(z,t)) \quad (\text{Eq. 1}) $$

where $C(z,t)$ is the local $O_2$ concentration, $r(C(z,t))$ is the local rate of the $O_2$ chemical consumption (commonly called the “oxidation rate”), $D_{O2}$ is the coefficient of $O_2$ diffusion into
the thickness of the part, $t$ is the time of exposure, and $z$ is the spatial co-ordinate (i.e. the depth beneath the surface of the part).

Therefore, the challenge consists in determining, on the one hand, the mathematical expression of the local oxidation rate $r(C)$ and, on the other hand, the temperature dependences of its different parameters, as well as the temperature dependence of the coefficient of $O_2$ diffusion $D_{O_2}$. Such a challenge is not trivial, particularly at high temperatures, due to the possible triggering of additional thermal decomposition reactions involving chemical bonds generally considered stable at lower temperatures [21], but also due to the well-known impact of the molecular mobility on the oxidation rate [22] and $O_2$ diffusivity [23, 24] when approaching physical transitions such as $T_g$.

Fortunately, a sound expression of $r(C)$ was already proposed for epoxy-diamine (EPO-DA) matrices by analytically solving the system of differential equations derived from the so-called “closed-loop” mechanistic scheme [25]. The reliability of this expression was carefully checked for six different EPO-DA matrices over wide ranges of temperature (typically between 25 and 200 °C) and $O_2$ partial pressure (between 0.16 and 20 bars) [22, 25, 26]. It should be pointed out that this checking was not only performed in the glassy domain, i.e. from temperatures close to $T_g$ up to temperatures very far below $T_g$ (typically $T_g - 140$ °C), but also within the glass transition zone and at the beginning of the rubbery domain ($T_g + 40$ °C). Thanks to a careful kinetic analysis of all experimental data, it was shown that the expression of $r(C)$ can take into account the effect of the molecular mobility through the temperature dependences of its different parameters. This expression can be written as follows [25]:

$$r(C) = 2r_0 \frac{\beta C}{1+\beta C} \left(1-\frac{\beta C}{2(1+\beta C)}\right) \left(\frac{1}{1+b \text{Exp}(-Kt)}\right)^2 \quad (\text{Eq. 2})$$

In Equation 2, $\beta^{-1}$ roughly corresponds to the critical $O_2$ concentration $C_c$ above which $O_2$ excess is reached, and $r_0$ is the oxidation rate in $O_2$ excess (i.e. when $C \gg C_c$) in the steady
state regime ($t \to \infty$), as schematically shown in Figure 1. These two parameters depend on the concentration of oxidizable sites ([PH]) and the propagation ($k_2$ and $k_3$) and termination rate constants ($k_5$ and $k_6$) of the “closed-loop” mechanistic scheme as follows [25]:

$$ \beta = \frac{k_6 k_2}{k_5 k_3 [PH]} \quad (\text{Eq. 3}) $$

$$ r_0 = \frac{k_3^2 [PH]^2}{4k_6} \quad (\text{Eq. 4}) $$

In addition, $K$ is an apparent first-order rate constant depending not only related to these first two parameters (i.e. $\beta$ and $r_0$), but also to the initiation rate constant $k_{1b}$ [25]:

$$ K = 2(r_0 k_{1b})^{1/2} \left( \frac{\beta C}{1 + \beta C} \right)^{1/2} \quad (\text{Eq. 5}) $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Shape of the curve: oxidation rate versus $O_2$ concentration.</figcaption></figure>

As $k_{1b}$ is an elementary rate constant, its temperature dependence obeys a general Arrhenius' law for all EPO-DA matrices regardless of their physical state (i.e. glassy or rubbery) [22]:

$$ k_{1b} = 1.9 \times 10^9 \text{Exp} \left[ \frac{-90\,000}{RT} \right] \quad (\text{in L.mol}^{-1}.\text{s}^{-1}) \quad (\text{Eq. 6}) $$

where $T$ is the temperature expressed in $K$, and $R$ is the universal gas constant approximated to 8.314 J.mol$^{-1}$.K$^{-1}$.
For this reason, it was concluded that $k_{1b}$ is not (or only very slightly) impacted by the molecular mobility.

Although $\beta$ and $r_0$ depend on several elementary rate constants (each of them characterized by its own activation energy), their temperature dependence also seems to obey a general Arrhenius' law for all EPO-DA matrices in the glassy state [26]. According to recent numerical calculations [22], these Arrhenius' laws should be written as follows:

$$ \beta = \frac{2.5 \times 10^3}{[PH] \times F(T_g)} \times \text{Exp} \left[ \frac{48\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in L.mol}^{-1}) \quad (\text{Eq. 7}) $$

$$ r_0 = 6.4 \times 10^{-10} \times [PH]^2 \times F^2(T_g) \times \text{Exp} \left[ \frac{-112\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in mol.L}^{-1}.\text{s}^{-1}) \quad (\text{Eq. 8}) $$

where $T_g$ is the glass transition temperature expressed in K.

where $F(T_g)$ is an exponential function of $T_g$ previously defined in reference [22] as:

$$ F(T_g) = \text{Exp} \left( \frac{T_g - 323}{30} \right) \quad (\text{Eq. 9}) $$

When approaching $T_g$, however, the Arrhenius' graphs of $\beta$ and $r_0$ of some EPO-DA matrices clearly display a beginning of discontinuity [26]. In fact, such a behavior is expected to be general for all EPO-DA matrices because a jump of few decades was clearly evidenced when passing $T_g$ for two elementary rate constants: $k_3$ and $k_6$ [22]. According to recent numerical calculations [22], in the rubbery state, the Arrhenius' laws could rather be written as follows:

$$ \beta = \frac{4.7 \times 10^7}{[PH] \times F(T_g)} \times \text{Exp} \left[ \frac{64\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in L.mol}^{-1}) \quad (\text{Eq. 10}) $$

$$ r_0 = 1.8 \times 10^{-13} \times [PH]^2 \times F^2(T_g) \times \text{Exp} \left[ \frac{-128\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in mol.L}^{-1}.\text{s}^{-1}) \quad (\text{Eq. 11}) $$

It should be noted, however, that this is only a very rough estimate as the glass transition zone is certainly very wide (at least between $-2 \times 10^{-3} K^{-1}$ and $+2 \times 10^{-3} K^{-1}$) so that the
rubbery domain is not very well defined. In addition, there are still too few experimental data on the thermal oxidation of EPO-DA matrices in the rubbery state to have an absolute confidence in these Equations 10 and 11.

One of the main objectives of this article is to now clearly highlight the effect of molecular mobility on $\beta$ and $r_0$ for all EPO-DA matrices. The question that naturally arises is whether this discontinuity alone can explain all the differences in behavior observed on either side of $T_g$ during the thermal ageing of EPO-DA matrices, such as differences in the penetration depth of oxidation, but also possible differences in the shape of the oxidation gradients. Indeed, it should be recalled that the molecular mobility is also expected to significantly impact the $O_2$ diffusivity, since $O_2$ permeation measurements tabulated in polymer handbooks suggest an increase in $D_{O_2}$ between one and two decades when moving from glassy to rubbery state for completely amorphous polymers [23, 24]. If this trend is confirmed, then an additional discontinuity should be considered to accurately simulate the oxidation gradients in EPO-DA matrices.

Unfortunately, there are still too few values of $D_{O_2}$ available in the literature for EPO-DA matrices to answer this issue. In addition, $O_2$ permeation measurements are usually performed at low temperature close to ambient, i.e. mainly in the glassy state, presumably to avoid undesired chemical consumption of $O_2$ during the test. The few $D_{O_2}$ values reported for EPO-DA matrices in the glassy state have been compiled in a recent publication [22], which led to distinguish two main families depending on the chemical structure of their diamine hardener (either aliphatic or aromatic):

a) Epoxies crosslinked with an aromatic diamine hardener for which:

$$D_{O_2} = 4.7 \times 10^{-10} \text{Exp} \left( -\frac{20\,000}{RT} \right) \text{m}^2 \cdot \text{s}^{-1} \quad (\text{in m}^2 \cdot \text{s}^{-1}) \quad (\text{Eq. 12})$$

b) Epoxies crosslinked with an aliphatic diamine hardener for which:
$$D_{O2} = 3.0 \times 10^{-5} \text{Exp} \left( -\frac{44,600}{RT} \right) \text{m}^2 \cdot \text{s}^{-1} \quad (\text{in m}^2 \cdot \text{s}^{-1}) \quad (\text{Eq. 13})$$

Until now, the impact of very stiff chemical structures, namely aromatic rings, on the $O_2$ diffusivity in the glassy state is not very well understood. Does it also occur in the rubbery state? This question remains completely open.

Finally, it should be emphasized that, to our knowledge, there is only one publication in which $O_2$ permeation has also been investigated in the rubbery state, but only on EPO-DA matrices crosslinked by an aliphatic diamine hardener [27]. According to the results found in this study, it would seem that the Arrhenius' law determined in the glassy state (i.e. Equation 13) could remain valid above $T_g$. Many more values of $D_{O2}$ in the rubbery state are needed to verify this trend.

The objective of this article is twofold:

i) First of all, outline a numerical calculation code based on Equations 1 and 2 and developed on the background of the finite element method (FEM) in order to simulate as accurately as possible the oxidation gradients in EPO-DA matrices.

ii) Then, compare the predictions of this code to experimental data compiled from the literature and completed by our own recent research works, in particular to the thicknesses of oxidized layer (TOL) measured with different physico-chemical techniques, such as: optical microscopy, Fourier transform Infra-Red (FTIR) mapping, and (nano or micro)-indentation.

It should be pointed out that, in the glassy state, the calculations will be done with the values of parameters $k_{1b}$, $\beta$, $r_0$ and $D_{O2}$ determined respectively with Equations 6, 7, 8 and 12 (or 13), i.e. established in previous publications [22, 25, 26]. In contrast, in the rubbery state, the values of parameters $\beta$, $r_0$ and $D_{O2}$ will gradually be changed if necessary, in order to adjust as faithfully as possible the oxidation gradients and thus, reveal the expected discontinuity for these three parameters in their Arrhenius' graph on either side of $T_g$.
Once the reliability of this code fully checked, it will then be possible to complete it with a
series of structure/property relationships (see for instance in reference [28]) in order to then
compute the gradients of several properties of practical interest for the aeronautical field, such
as T<sub>g</sub> and elastic properties (e.g. Young's modulus). Thus, the future publications will be
devoted to the careful analysis of these property gradients with the final objective of
determining the critical conditions for the initiation and/or propagation of cracks in the oxidized
layer.

## 2. Description of the calculation code

As explained in the introduction, the numerical calculation code is based on Equations 1 and 2
in order to simulate as accurately as possible the oxidation gradients in the EPO-DA matrices
regardless of their physical state (i.e. glassy or rubbery). In addition, this code is developed on
the background of FEM in order to write a strong coupling between the chemical consumption
of O<sub>2</sub> and its diffusion through the thickness of the part. This will allow us to consider the
possible effects of the mechanical loading state on the oxidation kinetics in future publications.
It should be pointed out that, thanks to this development choice, the code can then be easily
introduced into common commercial software for mechanical calculation and sizing, such as
Abaqus®.

### 2.1. Formulation and solving by a finite element method (FEM)

The Figure 2 shows the classical discretization scheme under finite elements for the calculation
of the gradient of O₂ concentration C(z, t) in the case of the unidirectional diffusion of O₂ along
the z axis (e.g. in the thickness of the part). In the z direction, the part is split into a high
number of elementary slices (or elements denoted e) with a thickness Δz and whose external
surfaces are positioned on two successive nodes, typically at z = zᵢ and z = zⱼ with j = i + 1.
According to this formalism, zᵢ ∈ [z₀ = 0, zₙ₋₁ = L - Δz] with i ∈ [0, n - 1], and zⱼ ∈
$[z_1 = \Delta z, z_n = L]$ with $j \in [1, n]$, $L$ and $n$ being the total thickness and the total number of elements of the part, respectively.

A temporal discretization is also performed over the entire period of thermal aging with a time increment $\Delta t$. Thus, thereafter, the time $t$ will appear with an index $k$, such as: $t_k \in [t_0 = 0, t_m = t_f]$ with $k \in [0, m]$, $t_f$ and $m$ being the final time and the total number of time increments, respectively.

<figure><img src="image_3.png" /><figcaption>Figure 2: Scheme for discretizing the spatio-temporal problem of thermal oxidation.</figcaption></figure>

Equation 1 can be rearranged into the following strong derivative form:

$$ \frac{\partial C}{\partial t} - D_{02} \frac{\partial^2 C}{\partial z^2} + f(C) = 0 \quad \text{in } \Omega = (0, L) \quad (\text{Eq. 14}) $$

where $\Omega$ refers to the spatial domain under study, and $f(C)$ represents the external forces of the problem to be solved. In the present case, this function is analogous to the oxidation rate $r(C)$:

$$ f(C) = 2r_0 \frac{\beta C}{1+\beta C} \left(1 - \frac{\beta C}{2(1+\beta C)}\right) \left(\frac{1}{1+b \text{Exp}(-Kt)}\right)^2 \quad (\text{Eq. 15}) $$

Equation 14 can be written as a variational minimization problem with boundary conditions whose weak integral formulation can be written as follows:
$$ \int_0^L \psi \operatorname{Res}(C) dz = \int_0^L \psi \left( \frac{\partial C}{\partial t} - D_{02} \frac{\partial^2 C}{\partial z^2} + f(C) \right) dz = 0 \quad (\text{Eq. 16}) $$

where $\psi$ is a weighting function at the nodes suitably chosen to minimize the residual $\operatorname{Res}(C)$ according to the Galerkin's approach.

The term in the middle of Equation 16 can be integrated by parts as follows:

$$ \int_0^L \psi D_{02} \frac{\partial^2 C}{\partial z^2} dz = \int_0^L D_{02} \psi \frac{\partial}{\partial z} \left( \frac{\partial C}{\partial z} \right) dz = \int_0^L D_{02} \left( \frac{\partial}{\partial z} \left( \psi \frac{\partial C}{\partial z} \right) - \frac{\partial \psi}{\partial z} \frac{\partial C}{\partial z} \right) dz \quad (\text{Eq. 17}) $$

$$ \text{i.e.} \quad \int_0^L \psi D_{02} \frac{\partial^2 C}{\partial z^2} dz = D_{02} \left[ \psi \frac{\partial C}{\partial z} \right]_0^L - \int_0^L D_{02} \frac{\partial \psi}{\partial z} \frac{\partial C}{\partial z} dz \quad (\text{Eq. 18}) $$

The introduction of Equation 18 into Equation 16 gives:

$$ \int_0^L \left( \psi \frac{\partial C}{\partial t} + D_{02} \frac{\partial \psi}{\partial z} \frac{\partial C}{\partial z} + \psi f(C) \right) dz - D_{02} \left[ \psi \frac{\partial C}{\partial z} \right]_0^L = 0 \quad (\text{Eq. 19}) $$

The last term of Equation 19 corresponds to the natural boundary conditions of the problem to be solved.

The local $O_2$ concentration $C(z_i, t) = C_i(t)$ at each node $i$ is given by the shape function $N_i$.

Consequently, the gradient of the $O_2$ concentration $C(z, t)$ across the thickness of the part can be approximated by a linear combination of all nodal values as follows:

$$ C(z, t) \approx \sum_{i=0}^{n} N_i(z) C_i(t) \quad (\text{Eq. 20}) $$

The same formalism can be used for the gradient of the weighting function $\psi(z, t)$:

$$ \psi(z, t) \approx \sum_{i=0}^{n} N_i(z) \psi_i(t) \quad (\text{Eq. 21}) $$

Let us now use the vector notation, which has the advantage of being much more compact, to represent the nodal values:

$$ N = \begin{bmatrix} N_0(z) \\ \vdots \\ N_n(z) \end{bmatrix} \quad (\text{Eq. 23}) $$

$$ C = \begin{bmatrix} C_0(t) \\ \vdots \\ C_n(t) \end{bmatrix} \quad (\text{Eq. 24}) $$
$$ \dot{\mathbf{C}} = \begin{bmatrix} \frac{\partial c_0(t)}{\partial t} \\ \vdots \\ \frac{\partial c_n(t)}{\partial t} \end{bmatrix} \qquad (\text{Eq. 25}) $$

$$ \mathbf{\psi} = \begin{bmatrix} \psi_0(t) \\ \vdots \\ \psi_n(t) \end{bmatrix} \qquad (\text{Eq. 26}) $$

$$ \mathbf{g} = D_{02} \begin{bmatrix} \left( \frac{\partial C(z,t)}{\partial z} \right)_0 \\ 0 \\ \vdots \\ 0 \\ \left( \frac{\partial C(z,t)}{\partial z} \right)_n \end{bmatrix} \qquad (\text{Eq. 27}) $$

Thus, Equations 20 and 21 become:

$$ C(z,t) \approx \mathbf{N}^T \mathbf{C} \qquad (\text{Eq. 28}) $$

$$ \psi(z,t) \approx \mathbf{N}^T \psi \qquad (\text{Eq. 29}) $$

Using these approximations, Equation 19 becomes:

$$ \begin{gathered} \psi^T \left( \int_0^L \mathbf{N} \mathbf{N}^T dz \right) \dot{\mathbf{C}} + D_{02} \psi^T \left( \int_0^L \frac{\partial \mathbf{N}}{\partial z} \frac{\partial \mathbf{N}^T}{\partial z} dz \right) \mathbf{C} + \psi^T \int_0^L \mathbf{N} f(C) dz \\ - \psi^T (\mathbf{N}(z=L) \mathbf{g}_n - \mathbf{N}(z=0) \mathbf{g}_0) = \mathbf{0} \end{gathered} \qquad (\text{Eq. 30}) $$

$$ \text{i.e.} \quad \left( \int_0^L \mathbf{N} \mathbf{N}^T dz \right) \dot{\mathbf{C}} + D_{02} \left( \int_0^L \frac{\partial \mathbf{N}}{\partial z} \frac{\partial \mathbf{N}^T}{\partial z} dz \right) \mathbf{C} = - \int_0^L \mathbf{N} f(C) dz + \mathbf{N}(z=L) \mathbf{g}_n - \mathbf{N}(z=0) \mathbf{g}_0 \qquad (\text{Eq. 31}) $$

where $\mathbf{N}^T$ is the transpose of the vector $\mathbf{N}$.

Finally, Equation 31 can be rewritten in the following matrix form:

$$ M \dot{\mathbf{C}} + K \mathbf{C} = \mathbf{F} + \mathbf{G} \qquad (\text{Eq. 32}) $$

where $M$ and $K$ are respectively the global mass and stiffness matrices of the problem to be solved, which must be computed from the nodal shape functions as follows:

$$ M = \int_0^L \mathbf{N} \mathbf{N}^T dz \qquad (\text{Eq. 33}) $$

$$ K = D_{02} \int_0^L \frac{\partial \mathbf{N}}{\partial z} \frac{\partial \mathbf{N}^T}{\partial z} dz \qquad (\text{Eq. 34}) $$
The same computation must be performed for the last two vectors **F** and **G**, describing respectively the chemical consumption of O₂ and the natural boundary conditions of the problem to be solved:

$$F = - \int_{0}^{L} N f(C) dz \quad (\text{Eq. 35})$$

$$G = N(z = L)g_n - N(z = 0)g_0 \quad (\text{Eq. 36})$$

For the sake of simplicity, linear functions of z were chosen for the nodal shape functions $N_i$ (and $N_j$). For practical convenience, an elementary basis is preferred. Then, it is assumed that each element e contains two half nodal shape functions ($N_i^e$ and $N_j^e$). Likewise, each nodal shape function can be calculated by a simple contribution of elementary shape functions of the two connected elements. Thus, discretizing the spatial domain in n elements, it is possible to apply Equations 33, 34 and 35 in order to calculate the elementary matrices and column vectors for each element e as follows:

$$m_{ij}^{e} = \begin{bmatrix} \int_{z_i}^{z_j} N_i^e(z) N_j^{eT}(z) dz & \int_{z_i}^{z_j} N_j^e(z) N_i^{eT}(z) dz \\ \int_{z_i}^{z_j} N_i^e(z) N_j^{eT}(z) dz & \int_{z_i}^{z_j} N_j^e(z) N_i^{eT}(z) dz \end{bmatrix} \quad (\text{Eq. 37})$$

$$k_{ij}^{e} = D_{O2} \begin{bmatrix} \int_{z_i}^{z_j} \frac{\partial N_i^e(z)}{\partial z} \frac{\partial N_j^{eT}(z)}{\partial z} dz & \int_{z_i}^{z_j} \frac{\partial N_j^e(z)}{\partial z} \frac{\partial N_i^{eT}(z)}{\partial z} dz \\ \int_{z_i}^{z_j} \frac{\partial N_i^e(z)}{\partial z} \frac{\partial N_j^{eT}(z)}{\partial z} dz & \int_{z_i}^{z_j} \frac{\partial N_j^e(z)}{\partial z} \frac{\partial N_i^{eT}(z)}{\partial z} dz \end{bmatrix} \quad (\text{Eq. 38})$$

$$f_{ij}^{e} = \begin{cases} - \int_{z_i}^{z_j} N_i^e(z) f(C) dz \\ - \int_{z_i}^{z_j} N_j^e(z) f(C) dz \end{cases} \quad (\text{Eq. 39})$$

The corresponding global matrices and column vectors over the entire spatial domain were then obtained by assembling these elementary bricks under study as follows:

$$M = \sum_{e=1}^{n} m_{ij}^{e} = \frac{\Delta z}{6} \times \begin{pmatrix} 2 & 1 & 0 & \cdots & 0 \\ 1 & 4 & \cdots & \cdots & \vdots \\ 0 & \cdots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & 4 & 1 \\ 0 & \cdots & 0 & 1 & 2 \end{pmatrix} \quad (\text{Eq. 40})$$
$$ K = \sum_{e=1}^{n} k_{ij}^{e} = \frac{D_{O2}}{\Delta z} \times \begin{pmatrix} 1 & -1 & 0 & \cdots & 0 \\ -1 & 2 & \cdots & \cdots & \vdots \\ 0 & \cdots & \cdots & \cdots & 0 \\ \vdots & \ddots & \ddots & 2 & -1 \\ 0 & \cdots & 0 & -1 & 1 \end{pmatrix} \quad (\text{Eq. 41}) $$

$$ F = \sum_{e=1}^{n} f_{ij}^{e} \quad (\text{Eq. 42}) $$

It should be noted that $M$ and $K$ are both tridiagonal matrices with a size $(n+1) \times (n+1)$, whereas $F$ is a column vector with a size $(n+1)$.

The Dirichlet's initial conditions (at $t=0$ at any node $i, j$) and the Neumann's boundary conditions (at any time $t$ on both surfaces of the part, i.e. at nodes $i=0$ and $j=n$) were chosen to solve Equation 32. In particular, the values of $C_i(0)$, $C_0(t)$ and $C_n(t)$ were directly be determined with the common Henry's law, and the value of $F_i(0)$ was easily deduced from $C_i(0)$ by applying Equation 15:

$$ C_i(0) = S_{O2} \times p_{O2} \quad (\text{Eq. 43}) $$

$$ C_0(t) = C_n(t) = S_{O2} \times p_{O2} \quad (\text{Eq. 44}) $$

$$ F_i(0) = -2r_0 \frac{\beta C_i(0)}{1+\beta C_i(0)} \left(1 - \frac{\beta C_i(0)}{2(1+\beta C_i(0))}\right) \left(\frac{1}{1+b}\right)^2 \quad (\text{Eq. 45}) $$

where $S_{O2}$ is the coefficient of $O_2$ solubility into the material under study, and $p_{O2}$ is the partial pressure of $O_2$ in the exposure environment.

It should also be noted that solving Equation 32 at $t=0$ gives direct access to the two terms

$$ g_0 = D_{O2} \left( \frac{\partial C(z,t)}{\partial z} \right)_0 \quad \text{and} \quad g_n = D_{O2} \left( \frac{\partial C(z,t)}{\partial z} \right)_n $$

Equation 32 was solved iteratively using the Crank-Nicolson's method. The gradients of $O_2$ concentration $C(z, t)$ in the thickness of the part, and their corresponding gradients of spatial $\partial C(z, t)/\partial z$ and temporal derivatives $\partial C(z, t)/\partial t$, were all determined in double precision according to the iterative scheme shown in Figure 3.
<figure><img src="image_4.png" /><figcaption>Figure 3. Iterative scheme for solving the spatio-temporal problem of thermal oxidation.</figcaption></figure>

## 2.2. Calculation of oxidation gradients and estimation of TOL

Then, the corresponding gradients of oxidation rate $r(C(z, t))$ were determined by applying Equation 2. It should be recalled that this latter quantity corresponds to the number of moles of $O_2$ consumed per volume unit of material and per time unit. This is the reason why the corresponding oxidation gradients $Q(C(z, t))$ were finally deduced by a simple integration of $r(C(z, t))$ with respect to time:

$$Q(C) = \int_{0}^{t} r(C(z, t))dt \quad (\text{Eq. 46})$$

In order to clearly explain in the next paragraphs how the thickness of the oxidized layer TOL was determined from the oxidation gradients, it seemed relevant to us to illustrate the demonstration with the simulations obtained for an EPO-DA matrix whose values of parameters $k_{1b}$, $\beta$, $r_0$ and $D_{O_2}$ have already been determined in previous publications [22, 26]. As an example, Figure 4 reports the $C(z, t)$ gradients computed in air at 200°C for the commercial EPO-DA matrix known under the reference: Tactix 123/Tactix 742-DDS.

Due to the choice of initial conditions (Eq. 43), $O_2$ is first homogeneously distributed throughout the thickness of the part (no gradient). Exposing the sample to high temperature causes the chemical consumption of $O_2$ throughout the thickness of the part. In the core of the
part, O₂ quickly reaches zero concentration. In the superficial layers of the part, in contrast, the supply of O₂ by diffusion from the external surfaces ultimately counteracts this chemical consumption and thus, generates the *C(z, t)* gradients. The shape of the final gradient (at *t* = *t*<sub>f</sub>) depends both on the boundary conditions (Eq. 44) and on the thermal oxidation kinetics, i.e. the reaction/diffusion coupling described by Equation 32.

It can be observed that, at any time *t*, the *C(z, t)* gradient is symmetrical due to the application of the same boundary conditions on the two external surfaces of the part (Eq. 44). As already observed by previous authors [8, 9, 11-20], *C(z, t)* is clearly chain-shaped. In addition, as shown in the insert at the right of Figure 4, *C(z, t)* reaches rapidly a limiting curve (after about 200 hours of exposure to air at 200°C), indicating that a steady-state regime is established. Of course, the duration to reach this steady-state (often called “oxidation induction time” in the literature) depends on the values of the input parameters, in particular on the values of *k*₁b *β*, *r*₀. In a previous publication, it was shown that this duration is a decreasing function of both temperature and O₂ partial pressure [25].

<figure><img src="image_5.png" /></figure>
Figure 4. Computation of the gradients of O₂ concentration in air at 200°C for different times of exposure in a 2 mm thick plate of Tactix 123/Tactix 742-DDS matrix.

Figure 5 reports the corresponding oxidation gradients $Q(z, t)$ obtained in air at 200°C for Tactix 123/Tactix 742-DDS matrix. As expected, at any time $t$, the $Q(z, t)$ gradient is also symmetrical in the thickness of the part. Contrarily to $C(z, t)$, it can be observed that $Q(z, t)$ is an increasing function of time because its derivative function $r(z, t)$ is positive. In the core of the part, as the initial concentration of O₂ is very low and quickly vanishes, $Q(z, t)$ quickly takes a very low and almost constant value such as:

$$Q_{min} = C(z, 0) = S_{O2} \times p_{O2} \quad (Eq. 47)$$

In the superficial layers of the part, in contrast, a fast oxidation can clearly be observed as a direct consequence of the diffusion control of the oxidation kinetics. Thus, $Q(z, t)$ is maximum on the two external surfaces of the part and minimum in its centre.

<figure><img src="image_6.png" /><figcaption>Figure 5. Computation of the oxidation gradients in air at 200°C for different times of exposure in a 2 mm thick plate of Tactix 123/Tactix 742-DDS matrix.</figcaption></figure>
At this point in the investigations, it should be recalled that the oxidation rate expression $r(C)$ was established for the low conversion ratios of the thermal oxidation reaction ($[PH] = [PH]_0 = constant$ in Eq. 2), i.e. without considering the progressive disappearance of the oxidizable sites (PH) during the exposure. As a direct consequence, once $r(z, t)$ has increased up to its steady-state value, it remains constant throughout the rest of the exposure. That is the reason why the increase in $Q(z, t)$ never slows down and thus, reaches much higher values than those usually measured experimentally in the superficial layers of the part.

Of course, the most relevant way to solve this problem would be to couple Equation 2 with a second equation evaluating the consumption of PH at any time $t$ at each position $z$ in the thickness of the part. In a first approach, it was decided to explore a simpler approach that consists in limiting $Q(z, t)$ to its maximum theoretical value, corresponding to the total concentration of available oxidizable sites in the EPO-DA matrix under study:

$$Q_{max} = [PH] \quad (Eq. 48)$$

In order to clearly show that the diffusion of $O_2$ from the external surfaces of the part is the main contribution to the oxidation gradient, well before the initial $O_2$ concentration within the part, the difference $Q(z,t) - Q(L/2, t)$ was computed in air at 200°C and plotted in Figure 6 for different depths in the thickness of the Tactix 123/Tactix 742-DDS matrix. All these curves clearly highlight the first three stages of the thermal oxidation kinetics, namely: the induction period, the sudden self-acceleration of oxidation, and the steady-state regime. The insert at the right of Figure 6 clearly shows that oxidation penetrates to a depth of at least 200 µm into the part.
<figure><img src="image_7.png" /><figcaption>Figure 6. Computation of the oxidation kinetics in air at 200°C at different depths in a 2 mm thick plate of Tactix 123/Tactix 742-DDS matrix.</figcaption></figure>

In addition, in order to better monitor the progress of the oxidation front in the thickness of the part and determine the corresponding values of TOL, it seemed relevant to us to normalize the oxidation gradients as follows:

$$q(z, t) = \frac{Q(z,t)-Q(L/2,t)}{Q_{max}-Q_{min}} \quad (\text{Eq. 49})$$

where $Q(z,t) - Q(L/2,t)$ is the contribution of $O_2$ diffusion to the oxidation gradient, and $Q_{max} - Q_{min}$ is its maximum amplitude.

Indeed, such a normalization has several advantages:

i) The maximum amplitude of the oxidation gradient is a value quite easy to determine for any EPO-DA matrix, as shown in Equations 47 and 48.

ii) The normalized quantity $q(z, t)$ can be easily related to the conversion ratio $\phi(z, t)$ of the oxidation reaction (see Equations 50 to 52).
iii) A critical value of $q(z, t)$, denoted for instance $q_C$, can then be selected to assess the penetration depth of oxidation (i.e. the TOL value) in the thickness of the part. To be relevant and robust for all the EPO-DA matrices, this oxidation criterion should be chosen within the range of the oxidation detection thresholds reported for the most common laboratory techniques.

Although the sensitivity threshold may vary from one experimental technique to another (typically between 0.1 and 1% for FTIR spectrometry and DSC), for the sake of simplicity, it was decided to choose a single value of $q_C$ equal to 0.2% for all the EPO-DA matrices in this study. Of course, this criterion value may be corrected later if it does not correctly account for all the experimental measurements obtained for this family of materials.

Figure 7 reports the corresponding normalized gradients $q(z, t)$ computed in air at 200°C for the Tactix 123/Tactix 742-DDS matrix. The insert at the right of Figure 7 shows clearly the progress of the oxidation front in the thickness of the part, and the horizontal dashed line indicated the position of the chosen value of oxidation criterion, i.e. $q_C = 0.2\%$, from which can be deduced the corresponding value of TOL. It is found that in air at 200°C the steady-state value of TOL is within the range 200-220 µm.
<figure><img src="image_8.png" /><figcaption>Figure 7. Computation of the normalized oxidation gradients in air at 200°C for different times of exposure in a 2 mm thick plate of Tactix 123/Tactix 742-DDS matrix. Application of the oxidation criterion of $q_c = 0.2\%$ for the estimate of TOL.</figcaption></figure>

However, as explained previously, values of $q(z,t)$ higher than unity have no physical meaning. That is the reason why a new quantity $\phi(z,t)$ was define to better describe the conversion ratio of the oxidation reaction. In fact, $\phi(z,t)$ correspond to $q(z,t)$ but without its non-physically realistic values, as shown in Figure 8. Schematically, the part can be split in the through-thickness direction into three successive zones with distinct oxidation conversion ratios (Figure 9):

- The “fully oxidized” zone where:

$$ \phi(z, t) = 1 \quad \text{if } 1 \le q(z, t) \quad (\text{Eq. 50}) $$

$O_2$ from the external environment passes through this zone without reacting with the polymer because there is no more oxidizable sites. This is an extremely brittle zone that
will easily crack under low mechanical loading, but may also crumble to give rise to the well-known phenomenon of chalking.

- The “still active oxidized” zone where:

$$ \phi(z,t) = q(z,t) \quad \text{if } q_c \le q(z,t) < 1 \quad (\text{Eq. 51}) $$

$O_2$ reacts with the polymer when it reaches this zone. The oxidation gradients and the corresponding gradients of the physical and mechanical properties continue to develop until the conversion ratio of the oxidation reaction reaches its maximum possible value given by Equation 50. Cracks will initiate and easily propagate in this zone when critical loading conditions (to be determined) are reached. $q_c$ represent the lowest threshold criterion corresponding to a virgin state. Then, it is assumed that material does not undergo aging before reaching $q_c$.

- The “healthy” zone where:

$$ \phi(z,t) = 0 \quad \text{if } 0 \le q(z,t) < q_c \quad (\text{Eq. 52}) $$

$O_2$ is just beginning to penetrate this zone. The conversion ratio is too low to have an impact on the physical and mechanical properties. It can be considered that this zone has not yet evolved during aging.
<figure><img src="image_9.png" /><figcaption>Figure 8. Normalized oxidation gradients of Figure 7 plotted with a limitation to unity for the values computed in the superficial layer of the plate. The horizontal dashed line in the insert at the right corresponds to the oxidation criterion $\phi(z, t) = q_c = 0.2\%$ for the estimate of TOL.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Schematization of the three successive zones with distinct oxidation conversion ratios in the superficial layer of the part.</figcaption></figure>

It should be noted that this concept of three successive zones was already proposed in the literature for polyimide matrices, which are generally used under harsher aging conditions than
epoxy matrices, with the objective of establishing relationships between the conversion ratio of the oxidation reaction and several properties of practical interest, such as the mass loss [29, 30], but also the oxidation induced shrinkage and the elastic properties [31, 32].

With this latest representation of the oxidation gradients, the steady state regime already observed for the *C*(z, t) gradient in Figure 4, is now also clearly visible for the *φ*(z, t) gradient in Figure 8. The corresponding changes in *TOL*, computed in air at 200 °C for the Tactix 123/Tactic 742 matrix, have been plotted in Figure 10 where they are compared to the experimental values determined by optical microscopy in a previous publication [28]. A nice agreement can be observed between the theory and experiment, thus highlighting the robustness of the proposed oxidation criterion: *φ*(z,t) = *q*<sub>c</sub> = 0.2%. It is now clarified that the steady-state value of *TOL* is about 215 µm.

<figure><img src="image_11.png" /><figcaption>Figure 10. Changes in the *TOL* of the Tactix 123/Tactix 742-DDS matrix in air at 200 °C.
Comparison between the experimental (full squares) and computed values (full line).</figcaption></figure>
### 3. Checking the validity of the calculation code

#### 3.1. Experimental values of TOL

The values of TOL measured on a wide variety of EPO-DA matrices after their thermal ageing in air between 110 and 200°C have been compiled from the literature and completed with our own recent research works (in particular on the DGEBA/DDS matrix), then plotted in the Arrhenius' diagram of Figure 11. It was decided to differentiate each epoxy matrix with a code referring to the chemical structure of the starting monomers, written as follows: "epoxy-hardener monomers". When the monomers are unknown, the EPO-DA matrix is designated with its trade name. It should be noted that these values of TOL were determined with three main physico-chemical techniques, commonly used to reveal the significant changes in the chemical composition or in the physical or mechanical behavior of the material after oxidation:

- Optical microscopy in polarized light allows a better observation of an oxidized layer, because it appears much brighter than the still intact core of the part,

- IR mapping allows the detection of the superficial layers where the oxidation products are concentrated,

- And (nano or micro)-indentation allows the detection of the changes of the elastic properties (in particular, the increase in the Young's modulus) in the oxidized layer.

Finally, it should be pointed out that, in order to clearly highlight the sole effect of molecular mobility, the values of TOL have not been plotted versus the reciprocal temperature, but rather versus the reciprocal gap to $T_g$. In each physical domain (i.e. both in glassy and rubbery states), it is found that TOL is a decreasing function of temperature because rising the temperature accelerates the chemical consumption of O₂ faster than its diffusion. In other words, if assuming that these two phenomena obey really an Arrhenius' law in these two domains, the activation energy is higher for the chemical consumption of O₂ than its diffusion. Indeed, it should be
recalled that, in each physical domain, the temperature dependence of TOL can be satisfyingly
described with a simple scaling law [9]:

$$
TOL = TOL_0 \exp \left( \frac{E_{TOL}}{RT} \right)^{1/2} \quad (\text{Eq. 53})
$$

$$
\text{with} \quad TOL_0 \approx \left( \frac{D_0 S p_{O2}}{r_0} \right)^{1/2} \quad \text{and} \quad E_{TOL} = \frac{1}{2}(E_r - E_D) \quad (\text{Eq. 54})
$$

$E_r$, $r_0$, $E_D$ and $D_0$ being the activation energy and the pre-exponential factor of the chemical
consumption of O₂ on the surface of the part and the diffusivity of O₂, respectively. As an erratic
variation of the solubility of O₂ with the temperature is often reported in the literature, it will
be assumed, as a first approximation, that this transport property remains constant regardless of
the temperature.

In Figure 11, a value of about 33 kJ.mol⁻¹ was determined for *E*<sub>TOL</sub> from the slope of the straight-line in the glassy domain. Knowing that *E*<sub>D</sub> is ranged between 20 and 44 kJ.mol<sup>-1</sup> for EPO-DA matrices (see Equations 12 and 13), it was thus deduced that *E*<sub>r</sub> is ranged between 86 and 110 kJ.mol<sup>-1</sup> below *T*<sub>g</sub>. This result is in nice agreement with previous literature findings for EPO-DA matrices [25].
<figure><img src="image_12.png" /><figcaption>Figure 11. Arrhenius' graph of the values of *TOL* measured on a wide variety of EPO-DA matrices after their thermal ageing in air. Experimental data recovered from 23 references [18, 28, 33-53] and completed with our own recent research works.</figcaption></figure>

For all the EPO-DA matrices, when the temperature exceeds the *T*<sub>g</sub> value, a positive jump of almost 2 decades in the pre-exponential factor of *TOL* is observed. As a result, *TOL* suddenly increases from about 20 µm to values approaching 1 mm when transitioning from the glassy to rubbery state. No doubt, such a jump is the result of molecular mobility. What is its exact origin? Can it only be explained by the discontinuities previously observed in the two Arrhenius' graphs of the parameters *β* and *r*<sub>0</sub>? Or should an additional discontinuity in the Arrhenius' graph of the O<sub>2</sub> diffusivity *D*<sub>O2</sub> need to be considered? The question is completely open.

In the rubbery state, roughly the same variation in temperature (i.e. with roughly the same *E*<sub>TOL</sub> value) seems to be observed for *TOL*. But, there are still too few experimental data above *T*<sub>g</sub> to be able to conclude with certainty.
## 3.2. Determination of the model parameters

Using the calculation code as an inverse solving method allowed the identification of the parameters $k_{1b}$, $\beta$, $r_0$ and $D_{O2}$ from the compilation of TOL values reported in Figure 11. Firstly, this simulation campaign allowed to check and confirm the orders of magnitude previously found for these parameters, mainly in the glassy state, by simulating as faithfully as possible the changes over time of exposure of two chemical quantities: the chemical consumption of $O_2$ and the build-up of carbonyl products [22, 25, 26]. Then, this simulation campaign allowed to determine new parameter values in temperature ranges that had little been studied until now, in particular within the glass transition zone and in the rubbery state. The Arrhenius' diagrams of Figure 12 to Figure 15 show all these values reported in both previous and current studies for parameters $k_{1b}$, $\beta$, $r_0$ and $D_{O2}$, respectively.

<figure><img src="image_13.png" /><figcaption>Figure 12. Arrhenius' graph of the values of $k_{1b}$ determined for a wide variety of EPO-DA matrices.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13. Arrhenius' graph of the values of β determined for a wide variety of EPO-DA matrices.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14. Arrhenius' graph of the values of r₀ determined for a wide variety of EPO-DA matrices.</figcaption></figure>

<figure><img src="image_16.png" /></figure>
Figure 15. Arrhenius' graph of the values of $D_{O2}$ determined for a wide variety of EPO-DA matrices.

The values of these parameters call for the following comments:

i) As already observed in a previous article [22], there is no impact (or only a very small impact, not detected in Figure 12) of the molecular mobility on the initiation of the oxidation reaction for EPO-DA matrices. The new results obtained in this study confirm the validity of the Arrhenius' law reported in introduction (Equation 6) for the corresponding rate constant $k_{1b}$.

ii) In contrast, there is a very significant impact of the molecular mobility on the termination and much less significant on the propagation of the oxidation reaction, thus leading to jumps of a few decades but of opposite directions when passing above $T_g$ for parameters $\beta$ and $r_0$. The directions of these jumps and their amplitudes can be predicted with Equations 3 and 4, respectively. For EPO-DA matrices, it can clearly be observed a positive jump of about 3 decades for $\beta$ in Figure 13, and a negative jump of about 2 decades for $r_0$ in Figure 14.

iii) Regarding the coefficient $D_{O2}$, the situation is more complex as $D_{O2}$ also depends on the chemical structure of the diamine hardener (either aliphatic or aromatic) of the EPO-DA matrices. In Figure 15, a positive jump between 1 and 2 decades can be observed when passing above $T_g$ for epoxies crosslinked with an aromatic hardener, whereas no jump is detected for epoxies crosslinked with an aliphatic hardener. It can therefore be concluded that, in addition to molecular mobility, aromatic rings have also an influence on the transport properties of $O_2$ in EPO-DA matrices in the glassy state, and that their influence disappears completely in the rubbery state.

The literature report that aromatic rings strongly stiffen the macromolecular chains, in particular the elastically active chains located between two crosslinking nodes, thus leading to higher
values of $T_g$. In networks, the chain stiffness is usually assessed using the flex parameter $F$, which is clearly an increasing function of the concentration of aromatic rings [54]. $F$ typically ranges from 14 g.mol⁻¹ for EPO-DA matrices without aromatic ring such as DGEBD-DA2 [54], to up to 33 g.mol⁻¹ for highly aromatic EPO-DA matrices such as Tactix 123/Tactix 742-DDS [28] or DGEBA-DDM [54]. Unfortunately, at this stage of the investigations, the effect of $F$ on $D_{O2}$ is not very well understood. Elucidating this effect and witting structure/O₂ diffusivity relationships would require performing O₂ permeation tests on a series of EPO-DA matrices with distinct concentrations of aromatic rings.

The following Arrhenius' laws can be determined from the slopes of the straight-lines in the glassy and rubbery states for parameters $\beta$ and $r_0$ in Figures 13 and 14, respectively:

- In glassy state:

$$ \beta = \frac{3.0 \times 10^3}{[PH] \times F(T_g)} \times \text{Exp} \left[ \frac{40\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in L.mol}^{-1}) \quad (\text{Eq. 55}) $$

$$ r_0 = 6.7 \times 10^{-10} \times [PH]^2 \times F^2(T_g) \times \text{Exp} \left[ \frac{-105\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in mol.L}^{-1}.\text{s}^{-1}) \quad (\text{Eq. 56}) $$

- In rubbery state:

$$ \beta = \frac{3.2 \times 10^6}{[PH] \times F(T_g)} \quad (\text{in L.mol}^{-1}) \quad (\text{Eq. 57}) $$

$$ r_0 = 2.7 \times 10^{-12} \times [PH]^2 \times F^2(T_g) \times \text{Exp} \left[ \frac{-91\,000}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in mol.L}^{-1}.\text{s}^{-1}) \quad (\text{Eq. 58}) $$

As expected, Equations 55 and 56 are not very far from Equations 7 and 8 reported in the introduction. Indeed, Equations 7 and 8 had been established from a lot of experimental data on the thermal oxidation of EPO-DA matrices in the glassy state [22]. The new experimental data reported in this study only serve to define these Arrhenius laws with even more precision.
In contrast, large discrepancies can be observed between Equations 57 and 10, as well as between Equations 58 and 11. As already mentioned in the introduction, the confidence in Equations 10 and 11 was very relative because they had provisionally been established while the rubbery domain was not yet well defined and contained too few experimental data. Thanks to the new experimental data provided in the rubbery state in this study, much more precise Arrhenius laws are now determined, which can of course be improved with the accumulation of experimental data in the future.

Regarding the coefficient $D_{O2}$, it is confirmed that, in the glassy state, $D_{O2}$ obeys the two different Arrhenius's law reported in the introduction (i.e. Equations 12 and 13) depending on the chemical structure of the diamine hardener (i.e. aliphatic or aromatic). However, in order to take into account the impact of the molecular mobility, these equations need to be reformulated from the slopes of the straight-lines shown in Figure 15. Their new expressions are:

- For epoxies crosslinked with an aliphatic diamine hardener:

$$D_{O2} = 7.6 \times 10^{-11} \text{Exp} \left[ \frac{-41800}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in } m^2 \cdot s^{-1}) \quad (\text{Eq. 59})$$

- For epoxies crosslinked with an aromatic diamine hardener:

$$D_{O2} = 2.7 \times 10^{-12} \text{Exp} \left[ \frac{-12900}{R} \left( \frac{1}{T} - \frac{1}{T_g} \right) \right] \quad (\text{in } m^2 \cdot s^{-1}) \quad (\text{Eq. 60})$$

However, this behavioral difference completely disappears when passing above $T_g$ in Figure 15. In a first approach, Equation 59 can be used to estimate the value of $D_{O2}$ for both EPO-DA families in the rubbery domain.

### 3.3. Comparison between numerical and experimental results

The experimental values of TOL reported in Figure 11 for a wide variety of EPO-DA matrices have been estimated with the calculation code detailed in Sections 2.1 and 2.2 using the values of parameters $\beta$, $r_0$ and $D_{O2}$ reported in Section 3.2, and an oxidation criterion $q_c = 0.2\%$ for
delimiting the boundary up to the oxidized layer in the sample. Figure 16 gives a comparison between these numerical calculations and the experiment.

This comparison highlights the reliability of the calculation code over the entire temperature range investigated. Considering the estimation uncertainties on the different parameters (namely on $\beta$, $r_0$ and $D_{O2}$) allows the satisfying prediction of the measurement uncertainties on TOL. Moreover, the positive jump of almost 2 decades observed for TOL when passing above $T_g$ is also correctly predicted by the calculation code. For epoxies crosslinked with an aliphatic diamine hardener, this jump results from both a positive jump of about 3 decades for $\beta$ (Figure 13) and a negative jump of about 2 decades for $r_0$ (Figure 14). For epoxies crosslinked with an aromatic diamine hardener, this jump is also the consequence of a positive jump between 1 and 2 decades for $D_{O2}$ (Figure 15).

<figure><img src="image_17.png" /></figure>
Figure 16. Arrhenius's graph of TOL values for a wide variety of EPO-DA matrices after their thermal ageing in air. Comparison between the experimental measurements (full circles) and the computed values (empty triangles).

**4. Conclusions and prospects**

A numerical calculation code has been developed by FEM to solve a diffusion/reaction problem and simulate the thermal oxidation gradients of EPO-DA matrices below and above their $T_g$. An experimentally-based oxidation criterion $q_c = 0.2\%$ was chosen to calculate the penetration depth of the oxidation in the sample thickness. The robustness of this criterion was proven by comparing the predicted and experimental values of TOL for a wide variety of EPO-DA matrices over a wide temperature range (Figure 4). The calculation code has been used as an inverse method in order to estimate the values of $D_{O2}$ from the different experimental values of TOL. A pronounced increase in $D_{O2}$ between one and two decades was evidenced for epoxy matrices crosslinked with an aromatic diamine hardener when the temperature exceeds their $T_g$, whereas no change in the linear Arrhenius' behavior was observed between the glassy and rubbery states for epoxy matrices crosslinked with an aliphatic diamine hardener. It was thus concluded that, in addition to molecular mobility, aromatic rings have also an influence on the transport properties of O₂ in EPO-DA matrices in the glassy state, and this influence disappears completely in the rubbery state. However, the lack of TOL measurements around $T_g$ requires more experiences in order to confirm this jump in $D_{O2}$ and enrich the Arrhenius graph of TOL. A deeper parametric study should be performed to take into account the uncertainties on the model parameters (namely on $\beta$, $r_0$ and $D_{O2}$) and examine their impact on the confidence interval of TOL estimates.

A series of perspectives can be considered for this study. First of all, the consequences of thermal oxidation on the thermomechanical properties of practical interest for the aeronautical
field, such as $T_g$ and elastic properties (in particular, Young's modulus $E$) could be calculated using well-known structure/property relationships [28]. Schematically, it is expected that predominant chain scissions (over post-crosslinking and oxidative cross-linking) cause a catastrophic reduction in $T_g$, while the antiplasticization effect leads to a significant increase in $E$. The development of $T_g$ and $E$ gradients in the oxidized layer should lead to its embrittlement and spontaneous cracking.

Then, these gradients could carefully be analysed in order to determine the critical conditions for crack initiation and/or propagation in the oxidized layer in the absence of external mechanical loading. Will these conditions correspond to critical concentrations of structural events directly affecting the mechanical behaviour of the EPO-DA matrices, *i.e.* chain scissions? Or rather to a combination of critical concentrations of several structural events, *e.g.* chain scissions and oxidative cross-linking? Will these conditions also depend on the local tensile stress state generated by oxidation? These questions are completely open.

In order to extend the calculation code to composite structures, the influence of carbon fibers and the residual thermal stresses (due to the mismatch of the thermal expansion coefficients of the matrix and the fibers) on the oxidation and damage kinetics could finally be considered.

All these steps are of paramount importance in the general approach for lifetime prediction of composite structures in the aeronautical field.

## References

[1] J.R. Kerr, J.F. Haskins, Effects of 50,000 h of thermal aging on graphite/epoxy and graphite/polyimide composites, AIAA J. 22(1) (1984) 96-102.

[2] K.N. Street, A.J. Russel, F. Bonsang, Thermal damage effects on delamination toughness of a graphite/epoxy composite, Compos. Sci. Technol. 32 (1988) 1-14.
[3] T.K. Tsotsis, S.M. Lee, Long-term thermo-oxidative aging in composite materials: Experimental methods, Compos. Sci. Technol. 32(11) (1998) 1115-1135.

[4] T.K. Tsotsis, S.M. Lee, Long-term thermo-oxidative aging in composite materials: Failure mechanisms, Compos. Sci. Technol. 58(3-4) (1998) 355-368.

[5] T.K. Tsotsis, S. Keller, K. Lee, J. Bardis, J. Bish, Aging of polymeric composite specimens for 5000 hours at elevated pressure and temperature, Compos. Sci. Technol. 61 (2001) 75-86.

[6] M.C. Lafarie-Frenot, S. Rouquie, Influence of the oxidative environments on damage in C/epoxy laminates subjected to thermal cycling, Compos. Sci. Technol. 64(10-11) (2004) 1725-1735.

[7] X. Colin, J. Verdu, Aging of organic matrix composite materials, In: Wiley encyclopedia of composites, 2nd edition, vol. 1/5, L. Nicolais, A. Borzacchielo and S.M. Lee eds, John Wiley & Sons Ltd, New York, Chap. 4, pp. 35-49, 2012.

[8] K.T. Gillen, R.L. Clough, Techniques for monitoring heterogeneous oxidation of polymers, In: Handbook of polymer science and technology, Vol. 2, M.P. Cheremisinoff ed., Marcel Dekker Inc., New York, pp. 167-202, 1989.

[9] L. Audouin, V. Langlois, J. Verdu, J.C.M. De Bruijn, Review: Role of oxygen diffusion in polymer ageing: Kinetic and mechanical aspects, J. Mater. Sci. 29 (1994) 569-583.

[10] X. Colin, A. Mavel, C. Marais, J. Verdu, Interaction between cracking and oxidation in organic matrix composites J. Compos. Mater. 39(15) (2005) 1371-1389.

[11] A.V. Tobolsky, D.J. Metz, R.B. Mesrobian, Low temperature oxidation of hydrocarbons. The phenomenon of maximum rates, J. Amer. Chem. Soc. 72 (1950) 1942-1952.

[12] C.R. Boss, J.C.W. Chien, Oxygen diffusion limitation in autoxidation of polypropylene, J. Polym. Sci.: Part A: Polym. Chem. 4(6) (1966) 1543-1551.
[13] G.C. Furneaux, K.J. Ledbury, A. Davis, Photo-oxidation of thick polymer samples. Part I: The variation of photo-oxidation with depth in naturally and artificially weathered low density polyethylene, Polym. Degrad. Stab. 3(6) (1981) 431-442.

[14] T. Seguchi, S. Hashimoto, K. Arakawa, N. Hayakawa, W. Kawakami, I. Kuriyama, Radiation induced oxidative degradation of polymers. Part I: Oxidation region in polymeric films irradiated in oxygen under pressure, Radiat. Phys. Chem. 17(4) (1981) 195-201.

[15] A.V. Cunliffe, A. Davis, Photo-oxidation of thick polymer samples, Polym. Degrad. Stab. 4(1) (1982) 17-37.

[16] S.G. Kiryushkin, Y.A. Shlyapnikov, Diffusion-controlled polymer oxidation, Polym. Degrad. Stab. 23 (1989) 185-192.

[17] K.M.B. Jansen, Analytical approximation of degradation profiles in polymer products, Polym. Eng. Sci. 34(21) (1994) 1619-1627.

[18] X. Colin, C. Marais, J. Verdu, A new method for predicting the thermal oxidation of thermosets. Application to an amine crosslinked epoxy, Polym. Test. 20(7) (2001) 795-803.

[19] X. Colin, C. Marais, J. Verdu, Thermal oxidation kinetics for a poly-(bismaleimide), J. Appl. Polym. Sci. 82(14) (2001) 3418-3430.

[20] X. Colin, J. Verdu, Strategy for studying thermal oxidation of organic matrix composites, Compos. Sci. Technol. 65 (2005) 411-419.

[21] X. Colin, G. Teyssèdre, M. Fois, Ageing and degradation of multiphase polymer systems, dans Handbook of Multiphase Polymer Systems, Vol. 2/2, A. Boudenne, L. Ibos, Y. Candau and S. Thomas eds, John Wiley & Sons Ltd, Chichester, Chap. 21, pp. 797-841, 2011.
[22]X. Colin, F. Essatbi, J. Delozanne, G. Moreau, Towards a general kinetic model for the thermal oxidation of epoxy diamine networks. Effect of the molecular mobility around the glass transition temperature », Polym. Degrad. Stab. 181, 109314, 2020.

[23]J. Crank, G.S. Park, Diffusion in Polymers, 4th Edition, Academic Press Inc., London, 1981.

[24]D.W. Van Krevelen, K. Te Nijenhuis, Properties of Polymers: Their correlation with chemical structure; their numerical estimation and prediction from additive group contributions, 4th Edition, Elsevier, Amsterdam, 2009.

[25]X. Colin, F. Essatbi, J. Delozanne, G. Moreau, A new analytical model for predicting the thermal oxidation kinetics of composite organic matrices. Application to diamine cross-linked epoxy, Polym. Degrad. Stab. 186 (2021) 109513.

[26]X. Colin, J. Delozanne, G. Moreau, New advances in the kinetic modelling of thermal oxidation of epoxy-diamine networks, Front. Mater. 8 (2021) 720455.

[27]M.C. Celina, A. Quintana, Oxygen diffusivity and permeation through polymers at elevated temperature, Polymer 150 (2018) 326-342.

[28]S. Terekhina, M. Mille, B. Fayolle, X. Colin, Oxidation induced changes in viscoelastic properties of a thermostable epoxy matrix, Polym. Sci. Ser. A 55(10) (2013) 614-624.

[29]J.D. Nam, J.C. Seferis, Anisotropic thermo-oxidative stability of carbon fiber reinforced polymeric composites, SAMPE Quart. 24(1) (1992) 10-18.

[30]I.M. Salin, J.C. Seferis, Anisotropic effects in thermogravimetry of polymeric composites, J. Polym. Sci. 31(8) (1993) 1019-1027.

[31]G.P. Tandon, K.V. Pochiraju, G.A. Schoeppner, Modeling the oxidative development in PMR-15 resin, Polym. Degrad. Stab. 91(8) (2006) 1861-1869.
[32] G.P. Tandon, K.V. Pochiraju, G.A. Schoeppner, Evolution of stress and deformations in high-temperature polymer matrix composites during thermo-oxidative aging, Mech. Time-Depend. Mater. (12) (2008) 45-68.

[33] X. Busch, M.E.R. Shanahan, Influence of the gaseous environment on the thermal degradation of a structural adhesive, J. Appl. Polym. Sci. 76(7) (2000) 987-992.

[34] T. Devanne, Radiochemical ageing of an epoxy network, PhD dissertation, ENSAM, Paris, France, 2003.

[35] N.Q. Ho, Damage of carbon/epoxy composites induced by thermo-oxidation of the matrix in the presence of mechanical stresses, PhD dissertation, ENSMA Poitiers, France, 2006.

[36] W. Trabelsi, Aging of carbon/epoxy composite materials for aeronautical applications, PhD dissertation, ENSAM, Paris, France, 2006.

[37] I. Ammar-Khodja, Thermal and thermo-oxidative aging of carbon/epoxy composites: relationship between chemical, physical and mechanical properties, PhD dissertation, Université du Havre, Le Havre, France, 2007.

[38] L. Olivier, N.Q. Ho, J.C. Grandidier, M.C. Lafarie-Frenot, Characterization by ultra-micro indentation of an oxidized epoxy polymer: Correlation with the predictions of a kinetic model of oxidation, Polym. Degrad. Stab. 93 (2008) 489-497.

[39] M.C. Lafarie, J.C. Grandidier, M. Gigliotti, L. Olivier, X. Colin, J. Verdu, J. Cinquin. Thermo-oxidation behaviour of composite materials at high temperatures: A review of research activities carried out within the COMEDI program. Polym. Degrad. Stab. 95(6) (2010) 965-974.

[40] M. Gigliotti, M. Minervino, M.C. Lafarie-Frenot, J.C. Grandidier, Effect of thermo-oxidation on the local mechanical behaviour of epoxy polymer materials for high temperature applications, Mech. Mater. 101 (2016) 118-135.
[41] E. Ernault, Thermo-oxidation of epoxy/amine networks, PhD dissertation, ENSAM, Paris, France, 2016.

[42] J. Cinquin, X. Colin, B. Fayolle, M. Mille, S. Terekhina, L. Chocinski-Arnault, M. Gigliotti, J.-C. Grandidier, M.-C. Lafarie-Frenot, M. Minervino, C. Cluzel, F. Daghia, P. Ladeveze, F. Zhang, Thermo-oxidation behavior of organic matrix composite materials at high temperatures", Adv. Aircr. Spacecr. Sci. 3(2) (2016), 171-195.

[43] M. Zhang, C. Zuo, B. Sun, B. Gu, Thermal ageing degradation mechanisms on compressive behavior of 3-D braided composites in experimental and numerical study, Compos. Struct. 140 (2016) 180-191.

[44] M. Zhang, B. Sun, B. Gu, Accelerated thermal ageing of epoxy resin and 3-D carbon fiber/epoxy braided composites, Compos. A 85 (2016) 163-171.

[45] M. Pecora, Y. Pannier, M.-C. Lafarie-Frenot, M. Gigliotti, Effect of thermo-oxidation on the failure properties of an epoxy resin, Polym. Test. 52 (2016) 209-217.

[46] M. Zhang, B. Sun, B. Gu, Meso-structure ageing mechanism of 3-D braided composite's compressive behaviors under accelerated thermo-oxidative ageing environment, Mech. Mater. 115 (2017) 47-63.

[47] J. Delozanne, Durability of epoxies: Application to aeronautical structural bonding, PhD dissertation, ENSAM, Paris, France, 2018.

[48] M. Pecora, Development of a cyclic indentation method for the characterisation of materials gradients in polymers and polymer composites due to thermal aging, PhD dissertation, ENSMA, Poitiers, France, 2018.

[49] B. Passilly, R. Delannoy, Characterization of the ageing of an epoxy resin using high temperature nanoindentation, Mater. Techn. 107 (2019) 206.
[50] E. Ernault, J. Dirrenberger, E. Richaud, B. Fayolle, Prediction of stress induced by heterogeneous oxidation: Case of epoxy/amine networks, Polym. Degrad. Stab. 162 (2019) 112-121.

[51] M. Pecora, O. Smerdova, M. Gigliotti, Gradients of cyclic indentation mechanical properties in PR520 epoxy and its 3D carbon fiber composite induced by aging at 150°C, Polym. Degrad. Stab. 193 (2021) 109720.

[52] A. Inamdar, Y.-H. Yang, A. Prisacaru, P. Gromala, B. Han, High temperature aging of epoxy-based molding compound and its effect on mechanical behavior of molded electronic package, Polym. Degrad. Stab. 188 (2021) 109572.

[53] A. Inamdar, M. van Soestbergen, A. Mavinkurve, W.D. van Driel, G.Q. Zhang, Modelling thermomechanical degradation of moulded electronic packages using physics-based digital twin, Microelectron. Reliab. 157 (2024) 115416.

[54] J.-P. Pascault, H. Sautereau, J. Verdu, R.J.J. Williams, Basic physical properties of networks, in: Thermosetting Polymers, Marcel Dekker Inc., New York, Chap 10, pp. 282-322, 2002.