# Simulation of thermal degradation in a composite material using phase field method

Mhadji Abdoussalam, Azdine Naït-Ali, Benjamin Batiot, Mathieu Calvat, Damien Halm

► To cite this version:

Mhadji Abdoussalam, Azdine Naït-Ali, Benjamin Batiot, Mathieu Calvat, Damien Halm. Simulation of thermal degradation in a composite material using phase field method. Composites Science and Technology, 2025, 261, pp.111015. 10.1016/j.compscitech.2024.111015. hal-04911977

HAL Id: hal-04911977

https://isae-ensma.hal.science/hal-04911977v1

Submitted on 24 Feb 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
Simulation of thermal degradation in a composite material using phase field method

M.Abdoussalam<sup>a</sup>, A.Nait-Ali<sup>a</sup>, B.Batiot<sup>b</sup>, M.Calvat<sup>c</sup>, D.Halm<sup>a,*</sup>

<sup>a</sup>Institut P', Physics and Mechanics of Materials Department, CNRS / ISAE-ENSMA / Université de Poitiers, 1 Av. Clément Ader, 86360 Chasseneuil-du-Poitou, France.

<sup>b</sup>Institut P', Fluids, Thermal Science and Combustion Department, CNRS / ISAE-ENSMA / Université de Poitiers, 1 Av. Clément Ader, 86360 Chasseneuil-du-Poitou, France.

<sup>c</sup>University of Illinois at Urbana-Champaign Materials Science and Engineering Department, Urbana, IL, USA

**Abstract**

Carbon fibers/epoxy resin composite laminates decompose by pyrolysis when submitted to high heat fluxes under inert atmosphere. A rigorous thermodynamic approach with internal variables has been adopted to better capture the phenomenon under study. In the literature, the models used to describe this degradation generally do not take into account the influence of microstructural heterogeneity on the decomposition rate, on the degradation kinetics and, consequently, on the propagation of the thermal front decomposition. To consider the variability of the composite microstructure, simulations at the microscale were conducted, involving a strong coupling between the evolution of the thermal degradation rate (which follows an Arrhenius law) and the temperature evolution using the phase-field method. An experimental approach with cone calorimeter has been undertaken to achieve two major objectives : first, to calibrate model parameters, and then to compare the numerical results with experimental data for the purpose of model validation. This comparison will focus on the analysis of degradation kinetics as well as the evolution of mass loss.

**Keywords:** Composite materials, microstructure, pyrolysis, thermal degradation, interface, coupling, gradient-based model, phase field

**1. Introduction**

Carbon fiber/epoxy resin composite materials are increasingly being used in aeronautical structures and compressed hydrogen storage due to their high stiffness, lightweight, and corrosion resistance properties Mori and Hirose (2009), Villalonga et al. (2009), Barral and Barthélémy (2006). However, their low fire resistance remains a significant limitation. Understanding the decomposition of these composite materials under fire, supported by relevant numerical models, could greatly enhance the design of composite structures in terms of safety, environmental impact, and cost reduction Halm et al. (2017). When a solid fuel is exposed to heat, it may experience mass loss due to thermal stress, combustion and heat transfer. In the case of carbonised matter, pyrolysis results in the production of combustion gases and a solid material in the form of 'char'. This new solid phase is primarily carbonised (burnt) and may undergo further decomposition. The degradation of carbon fiber reinforced polymers (CFRP) under fire results from several combined phenomena. Heat transfer is primarily conductive for limited fluxes, typically below 20

kW/m² Biasi et al. (2014), Milke and Vizzini (1991). When thermal fluxes exceed this threshold, pyrolysis or oxidation reactions occur, altering the material's composition. The composite's effective properties change during decomposition, depending on temperature and degradation state. These transformations produce solid residues and decomposition gases. The gas phase is carried within and released from the material's surface, potentially interacting with the surrounding atmosphere, including the possibility of ignition Mouritz and Gibson (2007). The objective is to produce a model that simulates transformation phenomena to observe the transition from the un-degraded solid material to the carbonized material at the microstructure scale. The model will be based on the condensed phases during pyrolysis, which will capture the transformation from the virgin material to the char material.

Thermal degradation phenomena have been extensively documented in the literature, particularly in the case of homogeneous materials. Various models have been developed over time to simulate the degradation of composites subjected to significant heat fluxes, primarily based on the work of Henderson et al. (1985), Henderson and Wiecek (1987). The models consider a decomposing composite material as a blend of virgin and charred materials. These methods have been greatly revised and enhanced by Florio Jr et al. (1991), Sullivan and Salamon (1992), Dimitrienko (1997), Galgano et al. (2009), who consider essential information from the literature.

<footnote>* Corresponding Author</footnote>

<footnote>Email addresses: mhadji.abdoussalam@ensma.fr
(M.Abdoussalam), azdine.nait-ali@ensma.fr (A.Nait-Ali), benjamin.batiot@univ-poitiers.fr (B.Batiot), mcalvat@illinois.edu (M.Calvat), damien.halm@ensma.fr
(D.Halm)</footnote>
mation on non-thermal equilibrium effects, thermal expansion, and effective stress modifications. While these models can describe temperature fields within the material with accuracy, they are constrained to one-dimensional configurations, presuming constant and well-understood heat fluxes. Significant advancements have been made in the research conducted by Biasi and colleagues Biasi et al. (2014), who introduced a model with the objective to improve the accuracy of predicting thermo-chemical degradation of polymer and fibre-based composites under non-uniform heating. However, the effect of heterogeneity on the evolution of degradation, material decomposition and char front propagation is not considered in the proposed models. Indeed, in this configuration, potential interface issues may be observed, and it is crucial to carefully consider these in the energy balance. The objective of this study is to present a model for simulating the decomposition of a composite material subjected to a heat flux at the microstructure scale, using the phase field method. Unlike conventional approaches, this model will incorporate the influence of microstructural heterogeneities on the propagation of the thermal degradation front. In fact, in the literature, the modeling of composite material degradation under thermal assault is predominantly based on approaches utilizing homogenized properties, where the composite is treated as an equivalent homogeneous material (Henderson et al. (1985), Henderson and Wiecek (1987), Florio Jr et al. (1991), Galgano et al. (2009), Biasi et al. (2014)). These macroscopic properties are derived from the aggregation of contributions from the various constituent phases. However, the influence of microstructural heterogeneity on the progression of thermal degradation has been rarely explored in the literature, despite its potential impact on the material's char formation. The model established in this work will focus on condensed phases during pyrolysis, capturing the transition from the virgin material to char. By accounting for microstructural variability, this model aims to provide key insights into how heterogeneities affect not only the spread of thermal degradation but also the formation and distribution of char zones, thereby influencing the composite material's properties. Moreover, it aims to manage the interfaces between degraded and non-degraded zones effectively.

A theory based on the thermodynamic approach will be introduced to describe the material's decomposition process under thermal loading, incorporating the phase field variable. This method is highly applicable for investigating models with a diffuse interface and non-negligible depth, where the local state of a heterogeneous microstructure is described by conservative and non-conservative fields associated with the nature of the phases (undamaged part of the material or degraded part) Chen (2002), Choudhury (2017), Yamanaka (2023), Gomez et al. (2019), Choudhury (2013). This approach will enhance our comprehension and management of the interface between degraded and undegraded regions, while taking into account the temporal evolution of properties with respect to degradation

through a strong coupling between temperature and thermal degradation.

To achieve this, we will initiate an experimental investigation to gather essential input and output data, which will be utilized for calibrating the model we will develop. The output parameters will include tracking the advancement of the char front and quantifying mass loss.

In the second phase of our work, we will formulate a theoretical model based on the principles of thermodynamics, thus introducing the phase field method into the resolution of this thermal degradation problem. Finally, our approach will culminate in a comparison between the numerical results derived from the model and the experimental result, aiming to validate the robustness and accuracy of the model.

## 2. Thermal degradation in a composite material

The objective is to modelling the thermal decomposition of the composite material during the pyrolysis phenomenon. The focus will be on the decomposition of the solid state, transforms the sound material into degraded material Koga et al. (2012), Rouquerol (1973), Gérard and Ono (2005), Brown et al. (1980). When a heat flux is applied to a solid's surface, a temperature front travels through the solid fuel. At first, only a small area is affected. However, as the thermal front moves through the material, a larger portion of the solid is heated. The speed of this thermal front decreases over time as it moves away from the heat source towards the cold back surface Colombiano et al. (2022). If the material being heated is a fuel, it will degrade and turn into a new solid when it reaches a specific temperature, known as the critical value. This new solid has disparate thermal and mechanical properties Torero (2016).

Models used for simulating such phenomena in heterogeneous materials, like composite materials, often disregard the effect of material heterogeneity and the interface between degraded and sound zones Biasi et al. (2014), Wang (2023), Dugast et al. (2020), Perron et al. (2021). The aim is to describe the thermal degradation of solids, taking into account the effects of microstructure, using the phase field method. The phase field method is a suitable approach for the study of diffuse interfaces. It has been applied mainly to solidification dynamics Boettinger et al. (2002), Gránásy et al. (2004), Echebarria et al. (2004), but also to other situations such as viscous fingering Folch et al. (1999), Casademunt (2011), fracture mechanics Bourdin et al. (2000), Amor et al. (2009), Bleyer and Alessi (2018), hydrogen embrittlement Martínez-Pañeda et al. (2018), Kristensen et al. (2020), Ngiam et al. (2022), and bubble dynamics Valizadeh and Rabczuk (2019), Ashour et al. (2021), Valizadeh and Rabczuk (2022), among others.
<table><caption>TABLE 1: Nomenclature and notations used</caption><thead><tr><th colspan="2">Nomenclature</th></tr></thead><tbody><tr><td><i>E</i></td><td>Internal energy (J)</td></tr><tr><td><i>Q</i></td><td>Received calorific power (J.s<sup>-1</sup>)</td></tr><tr><td><i>e</i></td><td>Specific internal energy (J.kg<sup>-1</sup>)</td></tr><tr><td><i>ε</i></td><td>Energy density (J.m<sup>-3</sup>)</td></tr><tr><td><i>S</i></td><td>Total entropy of the system (J.K<sup>-1</sup>)</td></tr><tr><td><i>E</i><sub>i</sub></td><td>Activation energy (J.mol<sup>-1</sup>)</td></tr><tr><td><i>r</i></td><td>Volumetric calorific power (J.m<sup>-3</sup>.s<sup>-1</sup>)</td></tr><tr><td><i>m</i></td><td>Mass (kg)</td></tr><tr><td><i>q</i></td><td>Imposed heat flux (W.m<sup>-2</sup>)</td></tr><tr><td><i>q&#773;</i> = <i>q&#773;</i></td><td><i>q</i> vector</td></tr><tr><td><i>&epsilon;</i></td><td>Interface energy (J.m<sup>-1</sup>)</td></tr></tbody></table><table><tbody><tr><td><i>A</i></td><td>Pre-exponential factor (s<sup>-1</sup>)</td></tr><tr><td><i>T</i></td><td>Temperature (K)</td></tr><tr><td><i>w</i></td><td>Volumetric energy density (Jm<sup>-3</sup>)</td></tr><tr><td><i>t</i></td><td>Time (s)</td></tr><tr><td><i>&xi;</i></td><td>Entropy per unit volume (J.K<sup>-1</sup>.m<sup>-3</sup>)</td></tr><tr><td><i>&alpha;</i></td><td>Thermal degradation rate</td></tr><tr><td><i>&rho;</i></td><td>Density (kg.m<sup>-3</sup>)</td></tr><tr><td><i>C</i></td><td>Heat capacity (J.kg<sup>-1</sup>.K<sup>-1</sup>)</td></tr><tr><td><i>&lambda;</i></td><td>Thermal conductivity (W.m<sup>-1</sup>.K<sup>-1</sup>)</td></tr><tr><td><i>&nabla;</i> and <i>&nabla;</i></td><td>Gradient and divergence operator</td></tr><tr><td><i>&Delta;H</i></td><td>Reaction enthalpy (J.kg<sup>-1</sup>)</td></tr></tbody></table>

A partial differential equation replaces boundary conditions at the interface. This equation governs the evolution of a non-conservative variable, known as the phase field, which acts as an order parameter. The phase field takes two distinct values (for example, 0 and 1) in each of the phases, with a smooth transition between both values in the zone around the interface, which then diffuses with a finite width. For our study, we use the phase-field variable $\alpha$ to describe whether the material is degraded (turns into char) or remains sound. $\alpha$ is a function of position and time, and its rate is governed by an equation that is coupled to equations for pyrolysis. These equations allow us to estimate the temperature field $T(x, t)$. The interfaces between char and sound material are described by smooth but highly localised changes of the variable. With that, all of the field (temperature and the thermal degradation rate) parameters in a phase field model have a smooth transition across the interface and thus yield diffuse (finite width) interfaces, making this method particularly well-suited for simulating such phenomena Ansari et al. (2021) Ammar et al. (2009).

## 2.1. Experimental study

To establish a robust and accurate numerical model, an experimental study has been carried out to observe and understand the physical mechanisms, on the one hand, and to obtain data to validate the degradation model, on the other hand. The focus of this experimental study has been put on two main aspects: firstly, degradation kinetics, involving the tracking of the char front evolution and, secondly, quantification of mass loss. Two experimental test campaigns were carried out. The first campaign aimed to analyze and understand the thermal decomposition phenomenon of a composite material and to calibrate the parameters of the numerical model. A "Cal" prefix will be added to all samples used in this campaign. In order to calibrate certain model parameters, micrographic images of the "Cal" samples will be used in numerical simulations. The second test campaign involves samples labelled with a

"Val" prefix and focuses on model validation experiments. The microstructures of the 'Val' samples will be used in order to compare the numerical and experimental results and to validate the model.

### 2.1.1. Material and methods

The studied material is a carbon fiber /epoxy matrix composite containing 58% of unidirectional (UD) Toray T700 carbon fibres. It is a material commonly used in the manufacture of hyperbaric hydrogen tanks. The unidirectional composite lay up allows to easily observe the different phases on a sample edge under a microscope, including the fibres, matrix, and porous regions. Rectangular samples of the material with fibers perpendicular to the sample axis, measuring 100m x 25mm x 5mm (Figure 1), are obtained by cutting large plates by water jet. Before and after conducting the experimental tests, the microstructures of each sample are reconstructed using an optical microscope (see Figure 1).

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Sample with its microstructure visualised using an optical microscope.</figcaption></figure>

This figure shows the heterogeneity of the distribution of the different phases in the material. This makes it an ideal object for testing a model dedicated to the degradation of inhomogeneous materials.
The samples were exposed to 35 kW/m² heat flux for different time periods (100s, 120s) under nitrogen conditions (oxygen concentration less than 5%). Afterwards, the propagation rates of the char front and the extent of mass loss have been measured.

The thermal aggression is performed by using a cone calorimeter apparatus (ISO 5660) Dao et al. (2013) Brauskas (1984) ISO (2002).

The specimen was positioned on a hollow silicate holder that had an identical cavity size to the specimen for hassle-free insertion (Figure 2). The use of silicate is preferred as it has low conductivity and can limit heat transfer to the surrounding environment. The weight of each sample was measured before and after decomposition.

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Sample holder with sample.</figcaption></figure>

The specimen is positioned at first in a furnace maintained at a temperature between 500-600°C. After the oxygen content within the furnace falls below 5%, a heat flux of 35kW/m² is applied to the top surface of the specimen for different time periods. The selection of this flux value is driven by the objective of monitoring the progression of the char front in samples with a depth of 5 mm. Therefore, it was essential to choose a flux capable of inducing sufficient decomposition reactions throughout the entire depth of the sample, while avoiding a rapid and complete decomposition across the full depth.

### 2.1.2. Results

In order to monitor the propagation of the char front through the depth of the sample, the most degraded zone, usually the central area, was examined under the optical microscope. This allows to observe the interface between degraded and undegraded areas of the sample.

From the microstructure visualization shown in Figure 3, a gradient of thermal degradation is observed throughout the depth. It is apparent that the char front has expanded up to 2.7 mm within the sample's depth. However, the significant porosity in the degraded zone poses difficulties in visualizing the interface between the degraded and undegraded regions of the sample. The results of the initial experimental test campaign are summarised in Table 2 and Table 3.

<figure><img src="image_4.png" /><figcaption>FIGURE 3: Sample exposed to a heat flux of 35kW/m² for 120s</figcaption></figure>

<table><caption>TABLE 2: Remaining mass of calibration sample for 100s exposition</caption><thead><tr><td>Samples</td><td>Experimental remaining mass (%)</td><td>Maximum degraded depth (mm)</td></tr></thead><tbody><tr><td>CalB30</td><td>97</td><td>1.6</td></tr><tr><td>CalB31</td><td>96</td><td>1.8</td></tr><tr><td>CalB32</td><td>98,6</td><td>1.2</td></tr></tbody></table>

<table><caption>TABLE 3: Remaining mass of calibration sample for 120s exposition</caption><thead><tr><td>Samples</td><td>Experimental remaining mass (%)</td><td>Maximum degraded depth (mm)</td></tr></thead><tbody><tr><td>CalC41</td><td>94.3</td><td>2.5</td></tr><tr><td>CalC42</td><td>92.4</td><td>2.7</td></tr><tr><td>CalC43</td><td>96</td><td>2</td></tr></tbody></table>

It is important to note that during all the exposure period, none of the samples underwent full transformation into char throughout their entire depth. During the decomposition phase, only the matrix undergoes degradation. The maximum temperature reached during the tests was around 600°C, which is insufficient to trigger decomposition phenomena in the fibers Feih and Mouritz (2012). Therefore, the mass losses recorded in the tests are mainly due to matrix decomposition.

The results of this experimental campaign allow us to focus on setting up a model that can accurately reproduce the thermal degradation phenomenon in thermodynamic framework.

### 2.2. Development of thermal degradation model

This section adopts a rigorous thermodynamic approach with internal variables to better capture the phenomenon of thermal degradation. The goal is to obtain a consistent model capable of simulating thermal decomposition using the phase field method. It is important to emphasize that this study focuses exclusively on thermal decomposition.
and its impact on thermal properties. Mechanical properties will not be addressed within the scope of this article.

### 2.2.1. Fundamental principles of thermodynamics and introduction of phase fields

The local form of the energy conservation and the Clausius-Duhem inequality for the studied problem of temperature-driven thermal decomposition is as follows Lemaitre et al. (2020), Maitournam (2010) :

$$ \frac{\partial \mathcal{E}}{\partial t} = -\nabla \underline{q} + r : \text{energy conservation} \quad (1) $$

$$ \mathcal{D} = \underbrace{\dot{\xi}T - \dot{\mathcal{E}}}_{\mathcal{D}_1} - \underbrace{-\frac{1}{T}\nabla T \cdot \underline{q}}_{\mathcal{D}_2} \geq 0 : \text{Clausius Duhem inequality} \quad (2) $$

with $\mathcal{D}$ : Total dissipation; $\mathcal{E}$ = $\rho e$ : internal energy per unit volume; $r$ : volumetric calorific power (inductive heating) $\xi$ = $\rho s$ : Specific entropy per unit volume; $T$ : Temperature; $q$ : imposed heat flux received by conduction through the boundary; $\underline{\nabla}$ : gradient operator.

This dissipation can be split into two components : an intrinsic dissipation component $\mathcal{D}_1$ and a thermal component $\mathcal{D}_2$, assuming that each term is positive :

$$ \begin{cases} \mathcal{D}_1 = \dot{\xi}T - \dot{\mathcal{E}} \geq 0 \\ \mathcal{D}_2 = -\frac{1}{T}\nabla T \cdot \underline{q} \geq 0 \end{cases} \quad (3) $$

### 2.2.2. Generalized Standard Materials framework (GSM)

In order to describe the behavior of a given medium accurately, it is essential to determine the variables that measure the state of the medium. According to the local state method, the thermodynamic state of the medium at a specific moment and point is defined entirely by a group of state variables, which depend solely on the considered location Lemaitre et al. (2020) Maitournam (2010). Given that the time derivatives of these variables are immaterial to the state definition, this assumption infers that any evolution can be envisaged as a sequence of successive equilibrium states. The particular choice of the type and number of state variables allows a more or less comprehensive representation of physical phenomena. The state variables, also known as thermodynamic or independent variables, include both observable variables and internal variables of the system. In the current case, which deals with the transformation of a material into char, the total energy of the system at a particular point depends on temperature (T) and the phase field variable $\alpha$ which represents the structural degradation variable ($\alpha$ = 0 when the material is sound and $\alpha$ approaches 1 as the material degrades and transforms into char). The GSM framework

assumes the existence of a thermodynamic potential from which the state laws are derived. In this way, we define a specific free energy potential, denoted $w(T, \alpha)$, which depends on these state variables :

$$ \begin{aligned} \mathcal{E}(\xi, \alpha) &= w(T, \alpha) + T\xi \\ \dot{\mathcal{E}}(\xi, \alpha) &= \dot{w}(T, \alpha) + \dot{\xi}T + \xi\dot{T} \end{aligned} \quad (4) $$

By substituting the elements from (4) into (2), we obtain :

$$ -\left(\frac{\partial w}{\partial T} + \xi\right) \dot{T} - \frac{\partial w}{\partial \alpha} \dot{\alpha} - \frac{1}{T} \underline{\nabla} T \cdot \underline{q} \geq 0 \quad (5) $$

To guarantee that dissipation remains positive, we find the state laws, which can be expressed in terms that express the irreversibility of the phenomenon under study. : Lemaitre et al. (2020), Maitournam (2010) :

$$ \begin{cases} \xi = -\frac{\partial w}{\partial T} \\ A_{\alpha} = -\frac{\partial w}{\partial \alpha} \end{cases} \quad (6) $$

Thus, the Clausius-Duhem inequality can be expressed as follows :

$$ A_{\alpha} \dot{\alpha} + \underline{G} \underline{q} \geq 0 \quad (7) $$

with

$$ \underline{G} = T \underline{\nabla}(\frac{1}{T}) = -\frac{1}{T} \underline{\nabla} T \quad (8) $$

To establish the evolution law that governs the dissipation of energy, we postulate the existence of a dissipation potential $\Omega^*$ composed of two terms, one associated with thermal conduction and the other with the degradation phenomenon :

$$ \Omega^*(A_{\alpha}, \underline{G}; T, \alpha) = \Phi^*(A_{\alpha}) + \chi^*(\underline{G}) \quad (9) $$

where $\Omega^*$ is a convex function that is positive at all points and has a value of zero at the origin.

A large body of work deals with the evolution of the thermal degradation of solids subjected to a heat flux Batiot et al. (2016), Torero (2016), Batiot et al. (2021). From this literature, it is possible to propose an expression for the dissipation potential that is consistent with experimental measurements :

$$ \Phi^*(A_{\alpha}; \alpha, T) = A(1 - \alpha)e^{-\frac{E}{RT}} A_{\alpha} \quad (10) $$

The normality rule $\dot{\alpha} \in \partial\Phi^*(A_{\alpha})$, together with the mathematical aforementioned properties of the dissipation potential, enables us to infer the degradation internal variable rate as follows, which satisfies the principles of thermodynamics outlined above (in particular, the positivity of dissipation) Germain et al. (1983) :
$$ \dot{\alpha} = \frac{\partial \Phi^*(A_{\alpha}; \alpha, T)}{\partial A_{\alpha}} = A(1 - \alpha)e^{-\frac{E}{RT}} \quad (11) $$

Regarding the dissipation associated with thermal conduction, the following expression can be used Lemaitre et al. (2020) :

$$ \chi(\underline{G}; \alpha) = \frac{1}{2} \underline{G} F(\alpha) \underline{G} \quad (12) $$

where F is the tensor associated with the material's thermal diffusion properties.

A similar normality rule leads to :

$$ \underline{q} = \frac{\partial \chi(\underline{G}, \alpha)}{\partial \underline{G}} = F(\alpha)\underline{G} = -\frac{F(\alpha)}{T}\underline{\nabla}T \quad (13) $$

Assuming isotropic heat conduction in the direction of degradation front propagation, $\frac{F(\alpha)}{T}$ is reduced to a scalar $\lambda$, which is identified as the effective thermal conductivity of the material, whose value will vary depending on the degradation variable $\alpha$ :

$$ \underline{q} = -\lambda(\alpha)\underline{\nabla}T \quad (14) $$

Revisiting Equation (1), integrating Equation (4) into it, and using (6), while assuming there is no source term, leads at the following resulting equations :

$$ T\dot{\xi} = -\nabla\underline{q} + A_{\alpha}\dot{\alpha} \quad (15) $$

and

$$ \dot{\xi} = -\frac{\partial^2 w}{\partial \alpha \partial T} \dot{\alpha} - \frac{\partial^2 w}{\partial T^2} \dot{T} = \frac{\partial A_{\alpha}}{\partial T} \dot{\alpha} + \frac{\partial \xi}{\partial T} \dot{T} \quad (16) $$

The change in mass density due to temperature is assumed negligible compared to its change due to degradation, $\alpha$. This assumption, together with the definition of the entropy per unit volume $\xi = \rho s$, allows to write the specific heat $C = T \frac{\partial s}{\partial T}$ as follows :

$$ \rho C \approx T \frac{\partial \xi}{\partial T} = -T \frac{\partial^2 w}{\partial T^2} \quad (17) $$

By substituting (17) into (15) :

$$ \rho C \dot{T} = (A_{\alpha} - T \frac{\partial A_{\alpha}}{\partial T}) \dot{\alpha} - \nabla \underline{q} \quad (18) $$

As it can be seen in the previous lines, the equation governing the temperature evolution needs the expression of the free energy $w$ whose role is to reflect the underlying physics of the thermal decomposition phenomenon. The phase field method, which is commonly used to simulate diffusion phenomena, frequently employs a Ginzburg-Landau type

approach. The Ginzburg-Landau approach is based on a total free energy functional. This functional is expressed as the sum of a homogeneous free energy density contribution and a quadratic gradient term related to the interface energy Penrose and Fife (1990) Ammar et al. (2009). In the case we are dealing with here, this decomposition is written as :

$$ w(T, \alpha) = w_0(T, \alpha) + \frac{\epsilon}{2} \nabla \alpha \cdot \nabla \alpha \quad (19) $$

with the assumption that $\alpha$ and $\nabla \alpha$ are non-independent variables and $\epsilon$ is a positive constant related to the interface energy between degraded and sound areas. The following expression is assumed for $w_0$ :

$$ w_0(T, \alpha) = -\rho_0 \left[ C(\alpha) \left( T \ln \left( \frac{T}{T_0} \right) - T + T_0 \right) + (1-\alpha) \Delta H \left\langle \frac{T}{T_0} - 1 \right\rangle \right] \quad (20) $$

where : $C$ is the heat capacity, $\rho_0$ the density of the virgin material, $\Delta H$ the degradation reaction enthalpy, $T_0$ a reference temperature chosen as initial temperature and $\langle \cdot \rangle$ the McAulay brackets.

$$ \left\{ \begin{array}{ll} \langle x \rangle = x & \text{if } x > 0 \\ \langle x \rangle = 0 & \text{else} \end{array} \right. \quad (21) $$

The first term of $w_0$ represents the classical internal energy commonly employed in thermal analysis Maitournam (2010), which is dependent on temperature and thus on the degradation rate $\alpha$. The second term of $w_0$ accounts for the energy associated with the decomposition reaction. It is crucial to express the energy function in a way that can accommodate the transition between different material states (from virgin material to charred material), capturing all transition phenomena accurately. This energy formulation offers a particularly insightful approach. When $T \le T_0$ and thus $\alpha = 0$, indicating no decomposition reaction, the material retains only the energy associated with its heat capacity. However, as $T$ exceeds $T_0$ and $\alpha$ increases, thermal degradation initiates, resulting in a progressive alteration of the material properties until the material transforms into char. When $\alpha = 1$, the material is completely converted into char, signifying the completion of the reaction. At this stage, the remaining energy in the material is solely attributed to the heat capacity of the obtained char.

The quadratic gradient term in the free energy, associated with interface energy, captures the effects of microstructural heterogeneity. This aspect represents a key innovation and highlights the originality of this new model. When the $\epsilon_c$ term is set to zero, the model reduces to $w = w_0$, thus recovering the classical model from the literature commonly used to simulate material degradation under heat flux or elevated temperatures.
By using 21, equation 20 can be expressed in two forms :

$$ w_0(T, \alpha) \begin{cases} = -\rho_0 \left[ C(\alpha) \left( T \ln\left(\frac{T}{T_0}\right) - T + T_0 \right) \dots \right. \\ \left. \quad + (1-\alpha)\Delta H \left( \frac{T}{T_0} - 1 \right) \right] \\ \text{or} \\ = -\rho_0 \left[ C(\alpha) \left( T \ln\left(\frac{T}{T_0}\right) - T + T_0 \right) \dots \right] \end{cases} \quad \text{if } \frac{T}{T_0} - 1 > 0 \\ \text{else} \tag{22} $$

The case of the application of a flux to the edge of the sample dealt with in the experimental part corresponds to $\frac{T}{T_0} - 1 > 0$. From this expression for $w$ and $w_0$, it is possible to calculate the set of equations describing the evolution of degradation and temperature :

$$ \left\{ \begin{aligned} A_{\alpha} &= \rho_0 \left[ \frac{\partial C}{\partial \alpha} \left( T \ln \left( \frac{T}{T_0} \right) - T + T_0 \right) - \Delta H \left( \frac{T}{T_0} - 1 \right) \right] + \epsilon \Delta \alpha \\ \frac{\partial A_{\alpha}}{\partial T} &= \rho_0 \left[ \frac{\partial C}{\partial \alpha} \ln \left( \frac{T}{T_0} \right) - \frac{\Delta H}{T_0} \right] \\ \dot{\alpha} &= A(1-\alpha)e^{-\frac{E_i}{RT}} \end{aligned} \right. \tag{23} $$

The Gâteaux derivative of $\nabla\alpha \cdot \nabla\alpha$ which give the term of $\Delta\alpha$ (Laplace operator) is a common result which can be found in works related to phase-field modelling, f.ex. Nguyen et al. (2015). By substituting $A_\alpha$, $\frac{\partial A_\alpha}{\partial T}$, $\dot{\alpha}$, and $\nabla q$ into (18) and by using the classical assumption that the dissipation resulting from specific heat variation is negligible compared to that arising from the transformation into char (see Appendix for further details), Equation 18 leads to :

$$ \rho C \dot{T} - \rho_0 \Delta H A (1-\alpha) e^{-\frac{E_i}{RT}} - \epsilon A (1-\alpha) e^{-\frac{E_i}{RT}} \Delta \alpha = \lambda(\alpha) \Delta T \tag{24} $$

Equation 24 is the usual heat equation involving pyrolysis with an additional gradient term. Henceforth, we will refer to this model as the **Gradient Model (GM)**. The simulation of the degradation phenomena induced by a heat flux in a material, needs to solve the following system of equations for the **Gradient Model (GM)** :

$$ \left\{ \begin{aligned} \rho C \frac{\partial T}{\partial t} - \lambda(\alpha) \Delta T - \rho_0 \Delta H A (1 - \alpha) e^{-\frac{E_i}{RT}} - \epsilon A (1 - \alpha) e^{-\frac{E_i}{RT}} \Delta \alpha &= 0 \\ \dot{\alpha} &= A(1 - \alpha) e^{-\frac{E_i}{RT}} \end{aligned} \right. \tag{25} $$

When $\epsilon = 0$, the classical heat equation with pyrolysis utilized in literature in the absence of oxidation is returned Biasi et al. (2014)Batiot et al. (2016)Torero (2016).

$$ \left\{ \begin{aligned} \rho C \frac{\partial T}{\partial t} - \lambda(\alpha) \Delta T - \rho_0 \Delta H A (1 - \alpha) e^{-\frac{E_i}{RT}} &= 0 \\ \dot{\alpha} &= A(1 - \alpha) e^{-\frac{E_i}{RT}} \end{aligned} \right. \tag{26} $$

This model will be referred to as **no gradient model (NGM)**.

This work focuses on simulating thermal degradation at the microstructural scale, where significant heterogeneity is present. The NGM model, being a specific case of the GM model ($\epsilon = 0$), cannot take into account heterogeneity effects and is therefore not suitable for this purpose. A brief comparison between the GM and NGM models is included in the sensitivity analysis of the parameter $\epsilon$, highlighting its influence on the propagation of the degradation front.

### 2.2.3. Material properties and boundary conditions

A T700 carbon fibre reinforced epoxy matrix composite is used in this study. This material was characterised as part of the FireComp project Blanc-Vannet et al. (2019) on the fire resistance of hyperbaric hydrogen storage tanks Quach et al. (2016). To conduct microstructure scale numerical simulations, it is crucial to know the properties of each composite phase (matrix, fibres and pores). The data for the fibres is widely available in literature Biasi et al. (2014), Feih and Mouritz (2012). For reasons of industrial confidentiality, the exact composition of the matrix cannot be published. For this reason, literature data on an epoxy similar to the resin in this material was used. A technique based on the mixture law has been established to estimate the properties of the epoxy matrix and its char (Table 4). It is important to distinguish between the properties of the virgin composite and the char of the composite, as well as between the matrix property and its char. In the applications considered here, the temperatures range from 50°C to 600°C and only the matrix is capable of decomposing and transforming into char. As the microstructure studied here (see Figure 1) is composed of fibers oriented in the direction transverse to thermal loading, an inverse mixture law combining the properties of the sound material and the degraded one is used to determine thermal conductivity Kulkarni and Brady (1997), Rivière et al. (2016). As regard the density and the heat capacity their expression is found to follow a classical mixing law as explained in Torre et al. (1998), Rivière et al. (2016). The kinetic data ($A$, $\Delta H$, $E_i$) are obtained through thermogravimetric analysis (TGA) experiments that have been carried out with the same methodology as Quach et al. (2016). The parameters are categorized into two groups : the known thermal properties for each phase of the composite, summarized in Table 4, and the estimated properties of the matrix, which are approximated using a mixing law based on the known thermal properties. Additionally, the kinetic parameters are detailed in Table 5.

As previously explained, only the matrix undergoes decomposition and turns into char. Therefore, we assume that only the thermal properties of the matrix evolve with the degradation parameter $\alpha$ Tranchard et al. (2017).
<table><caption>TABLE 4: Established thermal properties of each material phase</caption>
  <tbody>
    <tr>
      <th colspan="2" style="text-align: center; font-weight: bold;">Thermal parameters known<br>at reference temperature T<sub>0</sub></th>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Composite properties Halm et al. (2017)</td>
    </tr>
    <tr>
      <td>C<sub>comp</sub></td>
      <td>952.75 J.kg<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&lambda;<sub>comp</sub></td>
      <td>0.48 W.m<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&rho;<sub>comp</sub></td>
      <td>1350 kg.m<sup>-3</sup></td>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Char composite properties Halm et al. (2017)</td>
    </tr>
    <tr>
      <td>C<sub>comp-char</sub></td>
      <td>4.15 &times; 10<sup>3</sup> J.kg<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&lambda;<sub>comp-char</sub></td>
      <td>0.48 W.m<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&rho;<sub>comp-char</sub></td>
      <td>1015 kg.m<sup>-3</sup></td>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Fiber properties Kim et al. (2003)</td>
    </tr>
    <tr>
      <td>C<sub>fib</sub></td>
      <td><sup>1</sup>&frasl;<sub>&rho;<sub>fib</sub></sub> (1.1995 &times; 10<sup>6</sup> + 7811.8T<sub>0</sub> - 7.200T<sub>0</sub><sup>2</sup>) J.kg<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&lambda;<sub>fib</sub></td>
      <td>3.569 + 0.0110T<sub>0</sub> W.m<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&rho;<sub>fib</sub></td>
      <td>1700 kg.m<sup>-3</sup> Wang (2023)</td>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Pores properties Gerzhova et al. (2019)</td>
    </tr>
    <tr>
      <td>C<sub>p</sub></td>
      <td>1.075 &times; 10<sup>3</sup> J.kg<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&lambda;<sub>p</sub></td>
      <td>0.0466 W.m<sup>-1</sup>.K<sup>-1</sup></td>
    </tr>
    <tr>
      <td>&rho;<sub>p</sub></td>
      <td>0.588 kg.m<sup>-3</sup></td>
    </tr>
  </tbody>
</table>

<table><caption>TABLE 5: Matrix properties and kinetic parameters: $F_{fib}$, $F_{mat}$, and $F_p$ represent the volume fractions of fibers, matrix, and pores, respectively.</caption>
  <tbody>
    <tr>
      <th colspan="2" style="text-align: center; font-weight: bold;">Estimated thermal properties</th>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Estimated properties of the matrix</td>
    </tr>
    <tr>
      <td>C<sub>mat</sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(C<sub>comp</sub> - C<sub>fib</sub>F<sub>fib</sub> - C<sub>p</sub>F<sub>p</sub>)</td>
    </tr>
    <tr>
      <td><sup>1</sup>&frasl;<sub>&lambda;<sub>mat</sub></sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(<sup>1</sup>&frasl;<sub>&lambda;<sub>comp</sub></sub> - <sup>F<sub>fib</sub></sup>&frasl;<sub>&lambda;<sub>fib</sub></sub> - <sup>F<sub>p</sub></sup>&frasl;<sub>&lambda;<sub>p</sub></sub>)</td>
    </tr>
    <tr>
      <td>&rho;<sub>mat</sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(&rho;<sub>comp</sub> - &rho;<sub>fib</sub>F<sub>fib</sub> - &rho;<sub>p</sub>F<sub>p</sub>)</td>
    </tr>
    <tr>
      <td colspan="2" style="text-align: center; font-style: italic;">Estimated properties of the matrix char</td>
    </tr>
    <tr>
      <td>C<sub>mat-char</sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(C<sub>comp-char</sub> - C<sub>fib</sub>F<sub>fib</sub> - C<sub>p</sub>F<sub>p</sub>)</td>
    </tr>
    <tr>
      <td><sup>1</sup>&frasl;<sub>&lambda;<sub>mat-char</sub></sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(<sup>1</sup>&frasl;<sub>&lambda;<sub>comp-char</sub></sub> - <sup>F<sub>fib</sub></sup>&frasl;<sub>&lambda;<sub>fib</sub></sub> - <sup>F<sub>p</sub></sup>&frasl;<sub>&lambda;<sub>p</sub></sub>)</td>
    </tr>
    <tr>
      <td>&rho;<sub>mat-char</sub></td>
      <td><sup>1</sup>&frasl;<sub>F<sub>mat</sub></sub>(&rho;<sub>comp-char</sub> - &rho;<sub>fib</sub>F<sub>fib</sub> - &rho;<sub>p</sub>F<sub>p</sub>)</td>
    </tr>
    <tr>
      <th colspan="2" style="text-align: center; font-weight: bold;">Kinetic parameters</th>
    </tr>
    <tr>
      <td>&Delta;H</td>
      <td>317 &times; 10<sup>2</sup> kJ.kg<sup>-1</sup></td>
    </tr>
    <tr>
      <td>E<sub>i</sub></td>
      <td>130 &times; 10<sup>3</sup> kJ.mol<sup>-1</sup></td>
    </tr>
    <tr>
      <td>A</td>
      <td>9.9 &times; 10<sup>8</sup> s<sup>-1</sup></td>
    </tr>
  </tbody>
</table>

$$
\begin{align}
\rho(\alpha) &= (1 - \alpha)\rho_{\text{mat}} + \alpha\rho_{\text{mat-char}}, \nonumber \\
C(\alpha) &= (1 - \alpha)C_{\text{mat}} + \alpha C_{\text{mat-char}}, \tag{27} \\
\frac{1}{\lambda(\alpha)} &= \frac{1}{\lambda_{\text{mat}}}\alpha + \frac{1}{\lambda_{\text{mat-char}}}(1 - \alpha) \nonumber
\end{align}
$$

with: $\rho_{mat}$, $C_{mat}$, $\lambda_{mat}$ respectively the density, the specific heat, and the conductivity of the virgin matrix. $\rho_{mat-char}$, $C_{mat-char}$, $\lambda_{mat-char}$ respectively the density, the specific heat, and the conductivity of matrix fully transformed into char. Properties values of virgin matrix and come from Table (see Table 4).

In addition to the matrix properties evolution, the conductivity and the thermal capacity of the fibers are assumed to follow Kim et al. (2003).

$$
\begin{equation}
\begin{aligned}
C_{fib}(T) &= \frac{1}{\rho_{fib}} (1.1995 * 10^6 + 7811.8T - 7.200T^2) && (28) \\
\lambda_{fib}(T) &= 3.569 + 0.0110T
\end{aligned}
\end{equation}
$$

where the temperature is in Kelvin (K)

## 2.3. Boundary conditions and numerical resolution

With a view to simulating sample exposure tests in a calorimeter cone, we consider a composite material that is exposed at one edge by heat flux $q = 35 \text{ kW/m}^2$ absorbed to the surface, while at the other edges, heat exchange with the external environment is prevented (zero flux) in order to replicate the experimental conditions (Figure 4). The selection of this flux value aims to facilitate monitoring of the char front progression within a moderate timeframe (80-120 s) in samples with a depth of 5 mm. Therefore, it was essential to choose a flux capable of inducing sufficient decomposition reactions throughout the entire depth of the sample, while avoiding a rapid and complete decomposition across the full depth. Furthermore, the work of Quach Quach et al. (2017) demonstrates that precise characterization of thermal degradation depends not on flux intensity but on the total energy delivered (E=flux*exposure time), i.e., the thermal flux exposure duration. Consequently, a single heat flux value is sufficient to effectively quantify the progression of thermal degradation within the material.

<figure><img src="image_5.png" /><figcaption>FIGURE 4: Scheme of the problem with boundary conditions.</figcaption></figure>

A finite difference method is used to solve the equations. The time and spatial discretization procedure is carried out as follows:

Time discretization :

$$
\frac{\partial T}{\partial t} = \frac{T_{i,j}^{n+1} - T_{i,j}^{n}}{\Delta t} \qquad (29)
$$

$$
\frac{\partial \alpha}{\partial t} = \frac{\alpha_{i,j}^{n+1} - \alpha_{i,j}^{n}}{\Delta t} \qquad (30)
$$
Spatial discretization :

$$
\Delta T = \frac{T_{i+1,j}^{n+1} - 2T_{i,j}^{n+1} + T_{i-1,j}^{n+1}}{\Delta x^2} + \frac{T_{i,j+1}^{n+1} - 2T_{i,j}^{n+1} + T_{i,j-1}^{n+1}}{\Delta y^2} \quad (31)
$$

$$
\Delta\alpha = \frac{\alpha_{i+1,j}^{n+1} - 2\alpha_{i,j}^{n+1} + \alpha_{i-1,j}^{n+1}}{\Delta x^2} + \frac{\alpha_{i,j+1}^{n+1} - 2\alpha_{i,j}^{n+1} + \alpha_{i,j-1}^{n+1}}{\Delta y^2} \quad (32)
$$

$n$ and $i, j$ are respectively the time and space discretisation indexes, $\Delta t$ and $\Delta x, \Delta y$ are respectively the time and space path.. By substituting equations (31) and (32) into the governing equations (25), the discretized form of the temperature equation is obtained.

$$
\begin{equation}
\begin{aligned}
T_{i,j}^n = & \left(\frac{-\lambda \Delta t}{\Delta x^2 \rho C}\right) T_{i-1,j}^{n+1} + \left(\frac{-\lambda \Delta t}{\Delta x^2 \rho C}\right) T_{i+1,j}^{n+1} \\
& + \left(1 + \frac{2\lambda \Delta t}{\rho C}\left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right)\right) T_{i,j}^{n+1} \\
& + \left(\frac{-\lambda \Delta t}{\Delta y^2 \rho C}\right) T_{i,j-1}^{n+1} + \left(\frac{-\lambda \Delta t}{\Delta y^2 \rho C}\right) T_{i,j+1}^{n+1} \\
& + \rho_0 \Delta H A(1-\alpha)e^{\frac{-E_i}{RT_{i,j}^{old}}} \\
& + \epsilon A(1-\alpha^n)e^{\frac{-E_i}{RT_{old}}} \left[ \frac{\alpha_{i+1,j}^n - 2\alpha_{i,j}^n + \alpha_{i-1,j}^n}{\Delta x^2} \right. \\
& \qquad \left. + \frac{\alpha_{i,j+1}^n - 2\alpha_{i,j}^n + \alpha_{i,j-1}^n}{\Delta y^2} \right]
\end{aligned}
\tag{33}
\end{equation}
$$

To simplify the solution and remain within the fra-
mework of a linear problem, the dissipation term is im-
plemented as a source term for the numerical resolution,
using $T_{old}$, which represents the temperature at the pre-
vious time step. This will allow to solve a matrix problem
of the form :

$$
T^n = KT^{n+1} \tag{34}
$$

where $T^{n+1}$ will be obtained by inverting the matrix $K$. $K$ is the finite difference matrix involving the material property and the discretisation parameters (see Table 6).

On the other hand, since thermal degradation is repre-
sented by an update function related to temperature, we
can discretize it and solve it explicitly. As long as the tem-
perature approximation at a node is accurate, the value of
the thermal degradation will be directly determined. Thus,
the discretized form of the thermal degradation equation
can be expressed as :

$$
\alpha_{i,j}^{n+1} = \Delta t A e^{\frac{-E}{RT_{i,j}^{n}}} + \left(1 - \Delta t A e^{\frac{-E}{RT_{i,j}^{n}}}\right) \alpha_{i,j}^{n} \quad (35)
$$

The goal is to solve the MG model equation, specified
as (25). The resolution procedure is carried out as follows :

- **Thermal Problem Resolution**: The first step involves solving the thermal equation to compute the temperature distribution across the domain at a specific time $t$. This calculation determines the temperature at each node of the microstructure.

- **Thermal Degradation Equation Update**: Using the computed temperature at time $t$, the degradation equation is solved step-by-step to evaluate the degradation rate $\alpha$ at each node.

- **Thermal Properties Update**: Once the thermal degradation is determined at each point, the thermal properties, which vary according to the degradation level, are then estimated locally.

This resolution methodology offers the advantage of relatively simple implementation but requires relatively short time steps to ensure result reliability. All the resolution procedure is summarized in Table 6

Find $T(\underline{x}, t)$, $\alpha(\underline{x}, t)$ that satisfy :

At $t=0: \forall \underline{x}$

$T(\underline{x}, 0) = T_{old} = 40^{\circ}C$

$\alpha(\underline{x}, 0) = 0$

At any time $t > 0$ :

Step 1 : Solving the temperature problem.

$$
\begin{equation}
\begin{aligned}
T_{i,j}^n = & \left(\frac{-\lambda \Delta t}{\Delta x^2 \rho C}\right) T_{i-1,j}^{n+1} + \left(\frac{-\lambda \Delta t}{\Delta x^2 \rho C}\right) T_{i+1,j}^{n+1} \\
& + \left(1 + \frac{2\lambda \Delta t}{\rho C}\left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right)\right) T_{i,j}^{n+1} \\
& + \left(\frac{-\lambda \Delta t}{\Delta y^2 \rho C}\right) T_{i,j-1}^{n+1} + \left(\frac{-\lambda \Delta t}{\Delta y^2 \rho C}\right) T_{i,j+1}^{n+1} \\
& + \rho_0 \Delta H A(1-\alpha)e^{\frac{-E_i}{RT_{old}}} \\
& + \epsilon A(1-\alpha)e^{\frac{-E_i}{RT_{old}}} \left[ \frac{\alpha_{i+1,j} - 2\alpha_{i,j} + \alpha_{i-1,j}}{\Delta x^2} \right. \\
& \qquad \left. + \frac{\alpha_{i,j+1} - 2\alpha_{i,j} + \alpha_{i,j-1}}{\Delta y^2} \right]
\end{aligned}
\tag{36}
\end{equation}
$$

$T^n = KT^{n+1}$

Step 2 : Solving the thermal degradation problem.

$$
\alpha_{i,j}^{n+1} = \Delta t A e^{\frac{-E}{RT_{i,j}^{n}}} + \left(1 - \Delta t A e^{\frac{-E}{RT_{i,j}^{n}}}\right) \alpha_{i,j}^{n} \quad (37)
$$

Step 3 : Update the thermal properties.

$$
\rho(\alpha)_{i,j}^{n} = (1 - \alpha_{i,j}^{n})\rho_{\text{mat}} + \alpha_{i,j}^{n}\rho_{\text{mat-char}}
$$

$$
C(\alpha)_{i,j}^{n} = (1 - \alpha_{i,j}^{n}) C_{\text{mat}} + \alpha_{i,j}^{n} C_{\text{mat-char}}
$$

$$
\frac{1}{\lambda(\alpha)}_{i,j}^{n} = \frac{1}{\lambda_{\text{mat}}} \alpha_{i,j}^{n} + \frac{1}{\lambda_{\text{mat-char}}} (1 - \alpha_{i,j}^{n})
$$

Step 4 : Update $T = T_{old}$ and $\alpha$

<caption>TABLE 6: Solve procedure of the equations</caption>

It should be highlighted that in this study only the
matrix is degraded during the degradation process, hence
$α = 0$ at the fibre and porosity level.

## 2.4. *Structure studied and simulation results*

The simulations are performed at the microstructure scale. To obtain a full depth microstructure map, a real microstructure was visualised under an optical microscope. Post-processing was then performed using the Fiji software to separate the different phases that make up the composite (Figure 1). The matrix is shown in red, the fibres in green and the pores in purple. This microstructural mapping helped us to better understand the morphology and spatial distribution of the fibres in the material. To ensure reliable simulation results, it is essential to identify all model parameters and understand their contributions to the studied phenomenon. Consequently, the results section of this study is organized into three parts :

### 1. Influence and Methodology for Approximating Interface Energy

This first paragraph outlines the methodology used to determine the value of the interface energy, $\epsilon_c$. The impact of this parameter on the propagation of the degradation front is analyzed. Based on literature recommendations suggesting very small regularization values, simulations are conducted on the calibration samples ("Cal") using various random values of $\epsilon_c$ within a narrow range. This approach identifies the optimal values of $\epsilon_c$ that accurately capture the char front progression, consistent with experimental results. These findings will be used to establish a law describing the evolution of interfacial energy as a function of microstructural parameters, enabling its prediction for any given microstructure.

### 2. Determination of the Representative Elementary Volume (REV)

This part describes the procedure for determining the REV, a critical step to ensure that the results are representative of the entire material, thereby enhancing the reliability of predictions at a larger scale.

### 3. Application to Validation Sample Microstructures

Once the REV and $\epsilon_c$ are determined, the methodology is applied to new microstructures from the validation sample category ("Val"). This phase involves performing simulations on these microstructures and comparing numerical results with experimental data from these samples, thereby validating the robustness of the methodology across diverse microstructures.

#### 2.4.1. *Influence and Methodology for Approximating Interface Energy* $\epsilon$

The $\epsilon$ value can be seen as a regularization parameter that defines the interface energy between the damaged and undamaged areas. The chosen value of $\epsilon$ can either accelerate or decelerate the propagation of the degradation front. The use of phase fields in a heterogeneous medium with strong microstructure gradients remains relatively uncommon in the scientific literature. In the literature, terms similar to those in $\epsilon$ have a regularising function and affect the thickness of the degradation front (or more generally the interface between two media) Ammar et al. (2009)Ansari et al. (2021)Chen (2002)Boettinger et al. (2002)Gránásy et al. (2004)Echebarria et al. (2004). In the context of cracking, this regularisation parameter takes the form of a characteristic length corresponding to the crack width Nguyen et al. (2015)Bourdin et al. (2000)Amor et al. (2009). Thus, increasing this parameter leads to a widening of the crack, while decreasing it leads to a narrower crack with an increasingly pronounced localisation of the phenomenon.

In order to understand the influence of $\epsilon$ in our case, simulations were conducted on a representative volume of a microstructure encompassing the entire depth of the sample. Details regarding the approach employed to identify this representative volume for each microstructure will be explained in Paragraph 2.4.2. The analysis of the simulation results has shown that an adaptation of the value of $\epsilon$ based on the degree of microstructural heterogeneity is necessary to obtain results consistent with experimental observations Khachaturyan (2013)Chen (2002). This dependence is directly related to the presence of a significant gradient of the variable $\alpha$ within the microstructure. In the temperature range considered here, the fibers do not undergo degradation ($\alpha = 0$) whereas the matrix transformation rate may vary between 0 (for the sound side) and 1 (for the degraded side). Regarding porosity, $\alpha$ is 0 in this phase, whatever the temperature. This situation results in a significant gradient of the variable $\alpha$ from one phase to another, leading to a rapid increase in the temperature field, thereby accelerating the degradation process within the microstructure. It is therefore essential to determine the appropriate regularisation energy values associated with each microstructure to ensure accurate simulations despite abrupt fluctuations in $\alpha$ values. Simulations were carried out on the same microstructure with a heat flux of $35\,\text{kW/m}^2$ for an exposure duration of 100 seconds while varying the values of $\epsilon$ to more clearly illustrate the effect of the regularisation energy value. The results of these simulations are shown in Figure 5 in which it is possible to observe the progression of the degradation field through the depth of the sample for different values of $\epsilon$:

The curves of Figure 5 illustrate the average value of $\alpha$ in the matrix along the height throughout the depth of
<figure><img src="image_6.png" /><figcaption>FIGURE 5: Char front propagation for different $\epsilon$ values in $Jm^{-1}$ for 100s exposure</figcaption></figure>

the microstructure. It is clear from this result that reducing the value of $\epsilon$ down to 0 results in a notable increase in the depth of the degradation zone. Precisely, at $\epsilon = 0$, which represents the NMG model (26), along the depth of the microstructure, the value of $\alpha$ is significantly higher at each point compared to other cases where $\epsilon$ is non-zero. A criterion is established whereby areas with $\alpha \ge 0.8$ are considered degraded. This allows for the measurement of the total depth of the affected zone. The establishment of the $\alpha$ criterion will be discussed in Section 3. To avoid excessive propagation of the propagation front, which is not evidenced by the experiment (see paragraph 2.1.2), the regularization energy penalises the total energy, slo wing down the degradation kinetics. Each microstructure has its own unique intrinsic value for $\epsilon$ which has to be determined by linking it to the heterogeneity. It is important to note that, during this first testing campaign with the 'Cal' samples, for those exposed to 100 seconds, the size of the degraded area varied between 1.2 mm and 1.8 mm 2. However, in the simulations, when $\epsilon = 0$, the size of the degraded area could reach up to 3.6 mm (see Figure 5), which does not align with the experimental results. This discrepancy highlights the limitations of the NGM model. Therefore, to accurately simulate the evolution of thermal degradation at the microstructural scale, it will be neces sary to adopt the GM model (25) for future simulations.

Identification is based on observations made on "Cal-Sample" samples from the first experimental campaign (see Section 2.1.1). A correlation has been found between the $\epsilon$ parameter and the average distance $D_m$ that separates the fibres and the porosities (see Figure 6) :

$$D_m = \frac{1}{n} \sum_{i=1}^{n} D_i = \frac{1}{n} (D_1 + \dots + D_n) \quad (38)$$

with $D_i$ is the shortest matrix length between two ad-

jacent fibres, two adjacent pores or a pore and an adjacent fibre, from edge to edge in the representative volume.

<figure><img src="image_7.png" /><figcaption>FIGURE 6: Determination of the average distance $D_m$</figcaption></figure>

In summary, the interfacial energy identification was performed in 5 steps :

- Selection of a representative volume (RV) for each sample.

- Calculation of $D_m$ by image processing.

- Simulation of char front propagation.

- Determination of the $\epsilon$ value required to obtain the correct depth of char obtained in Table 2,3.

- Plot the point $(D_m, \epsilon)$ on the graph Figure 7.

After estimating the average distance $D_m$ of a given microstructure, the interfacial energy can be determined using the linear equation derived from the Figure 7.

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Interface energy approximation.</figcaption></figure>

This provides a general estimation of the regularization energy that establishes a direct correlation between the extent of the phase separation zone and the value of the regularization energy. In accordance with the principle that a narrow phase separation and therefore, a steep gradient (i.e., a reduced $D_m$ resulting in a high $\nabla\alpha$), requires a higher regularization energy to penalize the energy (19). Accurate identification of this term is crucial for the reliability of the obtained results in microstructure-scale simulations.

### 2.4.2. Identification of a representative volume (RV)

The aforementioned procedure to identify the interface energy is based on image processing of a representative
volume. To ensure that reliable results are obtained at the scale of the real structure, the selection of an appropriate representative volume (RV) for each sample is essential. To carry out this identification, the entire depth (5mm) for each RV associated with the microstructure of each sample is considered. The aim is to find the minimum height required so that the simulation results for the propagation of the char front and of the temperature front in depth no longer depend on *h*. The height of the sample (the edge on which the heat flux is applied) is determined according to the following method : first, a window size is selected that encompasses the microstructure and covers the entire depth. Then the height *h* of this window is increased gradually, and simulations are run for each window under conditions similar to those shown in Figure 4 (except for the boundary conditions on lateral edges which are taken as periodic and not zero flux Figure 8).

<figure><img src="image_9.png" /><figcaption>FIGURE 8: Identification VR method.</figcaption></figure>

The aim is to determine the minimum window size (*h* minimum value) from which the temperature and degradation results stabilise (Figure 10, Figure 9).

Thus, the following results (Figure 9 and Figure10) are obtained for different values of height *h*.

It is noticeable from Figure 9 and 10 that an increase in window height leads to a gradual convergence of temperature and degradation profiles towards similar results. As shown for this particular microstructure, the minimum height required to obtain representative results is *h* = 0.65mm. It is important to emphasise that this value may vary depending on the microstructure under investigation. Therefore, a similar approach has been taken for each of microstructures to ensure representative results for all the simulations.

### 3. Comparison of numerical and experimental results

The methodologies for determining a Representative Volume (RV) and the interfacial energy values have been applied to the validation samples called "ValSample" (see 2.1) in the following. Although their microstructures are different, these validation samples are made of the same material as the calibration ones. The objective is to conduct a series of cone exposure tests and to compare the experimental observations with the corresponding simulations.

<figure><img src="image_10.png" /><figcaption>FIGURE 9: Thermal degradation profile in space (120s)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10: Temperature profile in space (exposure duration : 120s)</figcaption></figure>

This process aims to evaluate the accuracy and robustness of the developed numerical model. The goal is to compare the measured degradation kinetics and mass loss with the simulated ones under the same loading conditions (heat flux, exposure time). Tables 7 and 8 provide morphological data for the two samples considered here.

**Remark :** It is important to note that the simulation framework does not explicitly consider porosity growth resulting from decomposition. Instead, pore growth is taken implicitly into account in the degradation variable *α* which characterizes the overall thermal degradation process, including pore growth.

The profiles results of the degraded samples from simulation and experimental are superimposed on figures 11 and figure 12.

<table><caption>TABLE 7: Microstructure properties of "ValST1" sample</caption><tr><td>$F_{fib}$</td><td>= 0.61;</td><td>$F_{mat}$</td><td>= 0.29;</td><td>$F_{pores}$</td><td>= 0.10</td></tr><tr><td>$D_m$</td><td>= $1.04 * 10^{-6}m$;</td><td>$\epsilon$</td><td>= $8.6 * 10^{-4}J/m$</td></tr><tr><td colspan="4">number of discretization point</td><td>= 9 090 896 nodes</td></tr></table>
<figure><img src="image_12.png" /><figcaption>FIGURE 11: Experimental and numerical thermal degradation of "ValST1" sample exposed to a heat flux of 35kW/m² during 100s a) enlargement of the 3 windows which show the numerical degradation state α in each window; b) simulation of the degradation c) microstructure of the exposed sample; d) enlargement of the 3 windows which show the experimental degradation state in each window</figcaption></figure>

The curve plotting alpha against depth is determined by averaging, for a given absciss, the value of the degradation over the entire height of the RV. The same applies to the temperature curve. In the experimental case, the degradation rate (α) is determined by analyzing the microstructure of the degraded sample using optical microscopy. The α curve is reconstructed through a grayscale thresholding method applied to the degraded microstructure. Specifically, the parameter α represents the material's degradation rate, which is experimentally correlated with the grayscale intensity of the degraded sample. This correlation stems from the contrast variations observed in the images of degraded samples. For a 8 bit figure, grayscale values range from 0 (white) to 255 (black), where lower values typically indicate non-degraded region intermediate values correspond to moderately degraded areas, and higher values represent highly degraded regions (e.g., charred zones). By normalizing these data between 0 and 1, the experimental α curve will be extracted and superimposed on the numerical evolution of α.

It is noteworthy to observe in Figure 11 and Figure 12 that the propagation velocities of the char front in the

<table><caption>TABLE 8: Microstructure properties of "ValS01" sample</caption><tr><td>F<sub>fib</sub> = 0.59; F<sub>mat</sub> = 0.30; F<sub>pores</sub> = 0.11</td></tr><tr><td>D<sub>m</sub> = 1.04 * 10<sup>-6</sup>m; ε = 8.6 * 10<sup>-4</sup>J/m</td></tr><tr><td>number of discretization point = 12 014 400 nodes</td></tr></table>

<figure><img src="image_13.png" /><figcaption>FIGURE 12: Experimental and numerical thermal degradation of "ValS01" sample exposed to a heat flux of 35kW/m² during 120s a) enlargement of the 3 windows which show the numerical degradation state α in each window; b) simulation of the degradation c) microstructure of the exposed sample; d) enlargement of the 3 windows which show the experimental degradation state in each window</figcaption></figure>

numerical model are highly consistent with the experimental results. Indeed, the size of the region transformed into char is very similar. Qualitatively, the comparison of Figures 11b and 12b on the one side and Figure 11c and 12c on the other side proves the model is able to fairly capture the extent of the decomposed area. For a more quantitative comparison, it is first necessary to set a degradation threshold above which the material can be considered to have been transformed into char. To do this, we use the work of Quach et al. Quach et al. (2016), who carried out ATG tests on the same material and showed that the decomposition reaction of the epoxy matrix in
nitrogen is largely completed at around 450°C. Thus, if we consider a temperature of 400°C, close to the total decomposition temperature, over 80% of the material is transformed into carbon. At this stage, there is no visual difference between a degradation rate of 80% and 100%. Consequently, regions where the temperature exceeds 400°C correspond to a degradation extent of $α ≥ 0.8$, and these areas can be reasonably considered fully degraded. By using this criterion, the observed curves in 11 and 12 provide a more detailed insight into the propagation of the char front along the depth. In both exposure time configurations (100s and 120s), the simulated spatial evolution of the degradation closely matches the experimental observations throughout the material depth. As illustrated in Figure 11, for a 100s exposure, the size of the completely degraded zone, where $α ≥ 0.8$, is approximately 1 mm, in agreement with both numerical simulations and experimental findings. Similarly, for a 120s exposure duration (Figure 12), this degraded zone extends to approximately 3.5 mm. It should be noted that there are some disturbances in the experimental degradation curve (see Figure 12). These fluctuations can be attributed to the significant presence of pores within the microstructure (both pre- and post-degradation), which induce substantial fluctuations during grayscale thresholding. In contrast, the simulation directly plots the average $α$ value across the matrix along the height throughout the depth of the microstructure, resulting in a more stable curve.

In addition to the simulation of the char depth, by using the hypothesis that the mass variation is caused by the decomposition of epoxy matrix, it is possible to calculate the mass loss as a function of temperature and time. The remaining mass (expressed as a percentage of the initial mass) is determined from the following mixture law 39, taking into account that the material is a mixture of pores, fibres and matrix in a state of varying degrees of degradation:

$$ \%m = 100 * \frac{\rho_f V_f F_f + \rho_m (1 - \bar{\alpha}) F_m V_m + \rho_{char} \bar{\alpha} F_m V_{char}}{\rho_f V_f F_f + \rho_m V_m F_m} \quad (39) $$

where $ρ_f$, $ρ_m$ and $ρ_{char}$ are the densities of the fibres, matrix and char formed during degradation, respectively. $V_m$, $V_f$ and $V_{char}$ are the volumes occupied by matrix, fibres and char formed during degradation, respectively, and $F_m$ and $F_f$ are the volume fractions of matrix and fibres. $\bar{\alpha}$ is the average degradation of the material at each time which characterises the volume fraction of the matrix that is transformed into char throughout the structure. This average is calculated as follows (40):

$$ \bar{\alpha} = \frac{1}{V} \int_V \alpha(x, t) dV \quad (40) $$

with V the volume occupied by the matrix in the microstructure

The results of the mass loss are summarised in Tables 9 and 10:

<table><caption>TABLE 9: Remaining mass of validation sample for 100s exposure</caption><thead><tr><th>Samples</th><th>Experimental remaining mass (%)</th><th>Simulated remaining mass (%)</th></tr></thead><tbody><tr><td>ValST1</td><td>96</td><td>97.8</td></tr><tr><td>ValST4</td><td>96</td><td>95.4</td></tr><tr><td>ValB31</td><td>94</td><td>97.2</td></tr></tbody></table>

<table><caption>TABLE 10: Remaining mass of validation sample for 120s exposure</caption><thead><tr><th>Samples</th><th>Experimental remaining mass (%)</th><th>Simulated remaining mass (%)</th></tr></thead><tbody><tr><td>ValS01</td><td>93</td><td>92.2</td></tr><tr><td>ValS04</td><td>94</td><td>94.6</td></tr><tr><td>ValB4</td><td>91</td><td>96</td></tr></tbody></table>

Simulated mass losses are very close to the experimental measures, which highlights the model's ability to provide a reliable representation of material degradation.

## 4. Conclusion

The thermal degradation phenomena observed in composite materials are highly complex and result from the interaction of thermal and chemical processes. A model has been developed in this study to simulate these phenomena at the microstructure scale under a nitrogen atmosphere using the phase field method. This model provides, the possibility to take into account the microstructural heterogeneity for more accurate simulation of thermal decomposition. As observed when $ϵ = 0$, representing the NGM model (a classic model used to simulate material pyrolysis), the model fails to faithfully reproduce correctly the thermal decomposition at the microstructure scale: the degradation kinetics does not align with experimental observations. The NGM model, which is widely used in the literature, yields satisfactory results that agree with experimental observations for homogeneous or homogenized materials. Indeed, the effects of material heterogeneity are directly accounted for in the effective properties of the homogenized material, reducing the impact of strong property gradients within the material, which are significant at the microstructure scale. The model presented here offers a robust coupling methodology, considering direct interactions between the temporal evolution of temperature and degradation rate.

Nevertheless, it should be noted that certain situations may demonstrate a slight delay or advancement in degradation compared to the model's forecast, which can be
explained by several factors : the established numerical model does not take into account all the complex physical phenomena involved in actual material degradation, such as an explicit consideration of pore formation and growth. These factors may impact the propagation of the char front and contribute to discrepancies between experimental and numerical outcomes. Furthermore, it is difficult to precisely replicate the exact experimental conditions within a numerical model due to environmental variations. This includes pressure, and atmospheric composition that can affect material degradation. As a result, minor variations may be observed between the experimental and numerical results. Although some slight differences were observed in our results, the comparison between the numerical and experimental results confirms the ability of the model to accurately predict char front propagation and composite degradation.

Moreover, in this study, several simplifying hypotheses were adopted, including assumptions on material properties considered as isotropic, and idealized boundary conditions, to facilitate computational implementation. The properties of the matrix were determined using mixing laws, which, while practical, could be further refined through direct measurement on the material itself. Decomposition was considered to be associated with a single reaction, which simplifies the analysis but limits the representation of complex degradation processes. This approach can be extended by introducing multiple successive or competing reactions to better capture the underlying physics of degradation. From a numerical standpoint, it is important to highlight the limitations associated with the use of an alternating resolution scheme, which necessitates small time steps. While this approach is straightforward to implement and provides a strong degree of robustness, it would be beneficial to enhance the current code by parallelization or integration of adaptive time-stepping, could be considered in future developments to increase efficiency and reduce computational times.

It would be interesting to further explore the approximation of the interface energy in order to find an expression directly linking it to other microstructural and kinetic parameters. This would enable an adaptable expression of this parameter for any microstructure from different types of composite materials. To enrich the model further, the explicit integration of oxidation and pore growth phenomena would be an interesting step to validate the model's ability to account for other physical phenomena that may occur during the degradation of a composite material subjected to a heat flux. This approach would allow the simulation of more realistic and common scenarios under specific atmospheric conditions.

Finally, incorporating multiphysical coupling into the model would be a significant advancement, enabling the simultaneous simulation of various physical phenomena. For instance, coupling thermal decomposition with mechanical damage would allow the model to capture the reciprocal effects : the influence of thermal degradation on the mate-

rial's mechanical properties and the impact of mechanical damage on decomposition kinetics and thermal behavior. The ultimate objective is to employ this thermal degradation simulation tool at the microstructural scale to develop a methodology for optimizing the manufacturing processes of these composites. This approach would facilitate the selection of optimal microstructures, enhancing their resistance to fire-induced damage.

**Acknowledgments** - This work was partially funded by the French Government programs "Investissements d'Avenir" LABEX INTERACTIFS (reference ANR-11-LABX-0017-01) - Pprime Institute gratefully acknowledges "Région Nouvelle-Aquitaine" for their financial support to the reported work within the framework of agreement AAPR2021-2021-12236310. Computation were performed on the supercomputer facilities at the mésocentre de calcul SPIN hosted by Poitiers University (P-2019-BAFE-154)

We extend our deepest gratitude to Djimédo Kondo for his invaluable discussions and expertise on thermodynamics. We are sincerely thankful to Amélie Caradec for her contributions to the image process analysis. We also wish to express our appreciation to Marc Poisson for his assistance in conducting the Cone Calorimeter tests.

## Références

Ammar, K., Appolaire, B., Cailletaud, G., Feyel, F., Forest, S., 2009. Formulation éléments finis des modèles de champ de phases basée sur la théorie de l'équilibre des microforces. PlastOx 2007-Mécanismes et Mécanique des Interactions Plasticité-Environnement 2007, 117-129.

Amor, H., Marigo, J.J., Maurini, C., 2009. Regularized formulation of the variational brittle fracture with unilateral contact : Numerical experiments. Journal of the Mechanics and Physics of Solids 57, 1209-1229. doi:10.1016/j.jmps.2009.04.011.

Ansari, T.Q., Huang, H., Shi, S.Q., 2021. Phase field modeling for the morphological and microstructural evolution of metallic materials under environmental attack. Computational Materials 7, 143.

Ashour, M., Valizadeh, N., Rabczuk, T., 2021. Isogeometric analysis for a phase-field constrained optimization problem of morphological evolution of vesicles in electrical fields. Computer Methods in Applied Mechanics and Engineering 377, 113669. doi:https://doi.org/10.1016/j.cma.2021.113669.

Atkins, P.W., De Paula, J., Keeler, J., 2023. Atkins' physical chemistry. Oxford university press.

Babrauskas, V., 1984. Development of the cone calorimeter—a bench-scale heat release rate apparatus based on oxygen consumption. Fire and Materials 8, 81-95.

Barral, K., Barthélémy, H., 2006. Hydrogen high pressure tanks storages : overview and new trends due to h2 energy specifications and constraints. WHEC16 : 16 World Hydrogen Energy Conference, France.

Batiot, B., Rogaume, T., Collin, A., Richard, F., Luche, J., 2016. Sensitivity and uncertainty analysis of Arrhenius parameters in order to describe the kinetic of solid thermal degradation during fire phenomena. Fire Safety Journal 82, 76-90. doi:10.1016/j.firesaf.2016.03.007.

Batiot, B., Rogaume, T., Richard, F., Luche, J., Collin, A., Guillaume, E., Torero, J.L., 2021. Origin and justification of the use of the arrhenius relation to represent the reaction rate of the thermal decomposition of a solid. Applied Sciences 11, 4075.
Biasi, V., Leplat, G., Feyel, F., Beauchêne, P., 2014. Heat and mass transfers within decomposing carbon fibers/epoxy resin composite materials, in: 11th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, p. 2678.

Blanc-Vannet, P., Jallais, S., Fuster, B., Fouillen, F., Halm, D., van Eekelen, T., Welch, S., Breuer, P., Hawksworth, S., 2019. Fire tests carried out in fch ju firecomp project, recommendations and application to safety of gas storage systems. International journal of hydrogen energy 44, 9100–9109.

Bleyer, J., Alessi, R., 2018. Phase-field modeling of anisotropic brittle fracture including several damage mechanisms. Computer Methods in Applied Mechanics and Engineering 336, 213–236. doi: https://doi.org/10.1016/j.cma.2018.03.012.

Boettinger, W.J., Warren, J.A., Beckermann, C., Karma, A., 2002. Phase-field simulation of solidification. Annual review of materials research 32, 163–194.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48, 797–826.

Brown, M.E., Dollimore, D., Galwey, A.K., 1980. Reactions in the solid state. Elsevier.

Casademunt, J., 2011. Exploring topological singularities with phase-field methods: Centrifugal viscous fingering as a case study. The European Physical Journal Plus 126, 1–11.

Chen, L.Q., 2002. Phase-field models for microstructure evolution. Annual review of materials research 32, 113–140.

Choudhury, A., 2017. Phase-field modeling as a method relevant for modeling phase transformation during interdiffusion, in: Handbook of Solid State Diffusion, Volume 1. Elsevier, pp. 363–389.

Choudhury, A.N., 2013. Quantitative phase-field model for phase transformations in multi-component alloys. Ph.D. thesis. Karlsruhe, Karlsruher Institut für Technologie (KIT), Diss., 2012.

Colombiano, J., Batiot, B., Dréan, V., Richard, F., Guillaume, E., Rogaume, T., 2022. Numerical analysis of the characteristics of the decomposition zone of a burning wood sample under cone calorimeter and evaluation of the limiting process. Journal of Analytical and Applied Pyrolysis 168, 105752.

Dao, D.Q., Luche, J., Richard, F., Rogaume, T., Bourhy-Weber, C., Ruban, S., 2013. Determination of characteristic parameters for the thermal decomposition of epoxy resin/carbon fibre composites in cone calorimeter. International journal of hydrogen energy 38, 8167–8178.

Dimitrienko, Y.I., 1997. Thermomechanical behaviour of composite materials and structures under high temperatures: 1. materials. Composites Part A: Applied Science and Manufacturing 28, 453–461.

Dugast, G., Settar, A., Chetehouna, K., Gascoin, N., Marceau, J.L., Bouchez, M., De Bats, M., 2020. Experimental and numerical analysis on the thermal degradation of reinforced silicone-based composites: Effect of carbon fibres and silicon carbide powder contents. Thermochimica Acta 686, 178563.

Echebarria, B., Folch, R., Karma, A., Plapp, M., 2004. Quantitative phase-field model of alloy solidification. Physical review E 70, 061604.

Feih, S., Mouritz, A., 2012. Tensile properties of carbon fibres and carbon fibre–polymer composites in fire. Composites Part A: Applied Science and Manufacturing 43, 765–772. doi: https://doi.org/10.1016/j.compositesa.2011.06.016.

Florio Jr, J., Henderson, J.B., Test, F.L., Hariharan, R., 1991. A study of the effects of the assumption of local-thermal equilibrium on the overall thermally-induced response of a decomposing, glass-filled polymer composite. International Journal of Heat and Mass Transfer 34, 135–147.

Folch, R., Casademunt, J., Hernández-Machado, A., Ramirez-Piscina, L., 1999. Phase-field model for hele-shaw flows with arbitrary viscosity contrast. i. theoretical approach. Physical Review E 60, 1724.

Galgano, A., Di Blasi, C., Branca, C., Milella, E., 2009. Thermal response to fire of a fibre-reinforced sandwich panel: Model formulation, selection of intrinsic properties and experimental validation. Polymer degradation and stability 94, 1267–1280.

Gérard, N., Ono, S., 2005. Hydride formation and decomposition kinetics. Hydrogen in Intermetallic Compounds II: Surface and Dynamic Properties, Applications, 165–195.

Germain, P., Suquet, P., Nguyen, Q.S., 1983. Continuum thermodynamics. ASME Journal of Applied Mechanics 50, 1010–1020.

Gerzhova, N., Côté, J., Blanchet, P., Dagenais, C., Ménard, S., 2019. A conceptual framework for modelling the thermal conductivity of dry green roof substrates. BioResources 14, 8573–8599.

Gomez, H., Bures, M., Moure, A., 2019. A review on computational modelling of phase-transition problems. Philosophical Transactions of the Royal Society A 377, 20180203.

Gránásy, L., Pusztai, T., Warren, J.A., 2004. Modelling polycrystalline solidification using phase field theory. Journal of Physics: Condensed Matter 16, R1205.

Halm, D., Fouillen, F., Lainé, E., Gueguen, M., Bertheau, D., van Eekelen, T., 2017. Composite pressure vessels for hydrogen storage in fire conditions: Fire tests and burst simulation. International Journal of Hydrogen Energy 42, 20056–20070. doi:10.1016/j.ijhydene.2017.06.088.

Henderson, J., Wiecek, T., 1987. A mathematical model to predict the thermal response of decomposing, expanding polymer composites. Journal of composite materials 21, 373–393.

Henderson, J.B., Wiebelt, J., Tant, M., 1985. A model for the thermal response of polymer composite materials with experimental verification. Journal of composite materials 19, 579–595.

ISO, I., 2002. 5660-1: "reaction-to-fire tests-heat release, smoke production and mass loss rate-part 1 (cone calorimeter method)". International Organization for Standardization, International Standard Organization, .

Khachaturyan, A.G., 2013. Theory of structural transformations in solids. Courier Corporation.

Kim, S.K., Sung Jung, B., June Kim, H., Il Lee, W., 2003. Inverse estimation of thermophysical properties for anisotropic composite. Experimental Thermal and Fluid Science 27, 697–704. doi: https://doi.org/10.1016/S0894-1777(02)00309-6.

Koga, N., Šesták, J., Šimon, P., 2012. Some fundamental and historical aspects of phenomenological kinetics in the solid state studied by thermal analysis, in: Thermal analysis of Micro, Nano-and Non-Crystalline Materials: Transformation, Crystallization, Kinetics and Thermodynamics. Springer, pp. 1–28.

Kristensen, P.K., Niordson, C.F., Martínez-Pañeda, E., 2020. A phase field model for elastic-gradient-plastic solids undergoing hydrogen embrittlement. Journal of the Mechanics and Physics of Solids 143, 104093.

Kulkarni, M., Brady, R., 1997. A model of global thermal conductivity in laminated carbon/carbon composites. Composites science and Technology 57, 277–285.

Lemaitre, J., Chaboche, J.L., Benallal, A., Desmorat, R., 2020. Mécanique des matériaux solides-3e éd. Dunod.

Maitournam, M., 2010. Mécanique des structures anélastiques: Comportement, chargements cycliques et fatigue. Laboratoire de Mécanique des Solides (LMS)-École Polytechnique, polytechnique édition.

Martínez-Pañeda, E., Golahmar, A., Niordson, C.F., 2018. A phase field formulation for hydrogen assisted cracking. Computer Methods in Applied Mechanics and Engineering 342, 742–761. doi: https://doi.org/10.1016/j.cma.2018.07.021.

Milke, J.A., Vizzini, A.J., 1991. Thermal response of fire-exposed composites. Composites Technology and Research 13, 145–151.

Mori, D., Hirose, K., 2009. Recent challenges of hydrogen storage technologies for fuel cell vehicles. International journal of hydrogen energy 34, 4569–4574.

Mouritz, A.P., Gibson, A.G., 2007. Fire properties of polymer composite materials. volume 143. Springer Science & Business Media.

Ngiam, Y., Cao, Z., Huang, M., 2022. Understanding hydrogen embrittlement in press-hardened steel by coupling phase field and hydrogen diffusion modeling. Materials Science and Engineering: A 834, 142523.

Nguyen, T., Yvonnet, J., Zhu, Q.Z., Bornert, M., Chateau, C., 2015. A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of
their microstructure. Engineering Fracture Mechanics 139, 18–39. doi:10.1016/j.engfracmech.2015.03.045.

Penrose, O., Fife, P.C., 1990. Thermodynamically consistent models of phase-field type for the kinetic of phase transitions. Physica D: Nonlinear Phenomena 43, 44–62. doi:https://doi.org/10.1016/0167-2789(90)90015-H.

Perron, N., Prévereaud, Y., Dellinger, N., Balat-Pichelin, M., Annaloro, J., 2021. 3d modeling and simulation of the thermal degradation of composite material during space debris atmospheric reentry, in: 8th European Conference on Space Debris, ESA/ESOC.

Quach, T., Benelfellah, A., Batiot, B., Halm, D., Rogaume, T., Luche, J., Bertheau, D., 2016. Determination of the tensile residual properties of a wound carbon/epoxy composite first exposed to fire. Journal of Composite Materials 51. doi:10.1177/0021998316637419.

Quach, T.H.Y., Benelfellah, A., Batiot, B., Halm, D., Rogaume, T., Luche, J., Bertheau, D., 2017. Determination of the tensile residual properties of a wound carbon/epoxy composite first exposed to fire. Journal of Composite Materials 51, 17–29.

Rivière, L., Caussé, N., Lonjon, A., Dantras, É., Lacabanne, C., 2016. Specific heat capacity and thermal conductivity of peek/ag nanoparticles composites determined by modulated-temperature differential scanning calorimetry. Polymer Degradation and Stability 127, 98–104.

Rouquerol, J., 1973. Critical examination of several problems typically found in the kinetic study of thermal decomposition under vacuum. Journal of thermal analysis 5, 203–216.

Sullivan, R., Salamon, N., 1992. A finite element method for the thermochemical decomposition of polymeric materials—i. theory. International Journal of Engineering Science 30, 431–441.

Thiry-Muller, A., 2018. Modélisation de la décomposition thermique des solides. Ph.D. thesis. Université de Lorraine.

Torero, J., 2016. Flaming ignition of solid fuels. SFPE handbook of fire protection engineering, 633–661.

Torre, L., Kenny, J.M., Maffezzoli, A., 1998. Degradation behaviour of a composite material for thermal protection systemspart ii process simulation. Journal of materials science 33, 3145–3149.

Tranchard, P., Samyn, F., Duquesne, S., Estèbe, B., Bourbigot, S., 2017. Modelling behaviour of a carbon epoxy composite exposed to fire: Part i—characterisation of thermophysical properties. Materials 10, 494.

Valizadeh, N., Rabczuk, T., 2019. Isogeometric analysis for phase-field models of geometric pdes and high-order pdes on stationary and evolving surfaces. Computer Methods in Applied Mechanics and Engineering 351, 599–642. doi:https://doi.org/10.1016/j.cma.2019.03.043.

Valizadeh, N., Rabczuk, T., 2022. Isogeometric analysis of hydrodynamics of vesicles using a monolithic phase-field approach. Computer Methods in Applied Mechanics and Engineering 388, 114191. doi:https://doi.org/10.1016/j.cma.2021.114191.

Villalonga, S., Nony, F., Magnier, C., Yvernes, J., Thomas, C., Delmas, B., Mazabraud, P., 2009. Composite 700 bar-vessel for on-board compressed gaseous hydrogen storage, in: Proc. of 17th International Conference on Composite Materials, Edinburgh, UK.

Wang, X., 2023. Modelisation of the thermal degradation of composites in order to predict their fire reaction when they are exposed to flames. Ph.D. thesis. Normandie Université.

Yamanaka, A., 2023. Phase-field modeling and simulation of solid-state phase transformations in steels. ISIJ International 63, 395–406.

# Appendix

The aim of this appendix is to specify the hypothesis outlined in paragraph 2.2.2, acknowledging that the dissipation originating from specific heat variation is insignificant in comparison to that arising from the transformation into char.

From the relation (see Atkins et al. (2023)):

$$dH = C_p(T) dT \quad \text{with } C_p \text{ the heat capacity, } dH \text{ the free enthalpy} \quad (41)$$

when a medium is heated from $T_0$ to $T$, its enthalpy changes from $H(T_0)$ to

$$H(T) = H(T_0) + \int_{T_0}^{T} C_p(T) dT \quad (42)$$

It is assumed that no phase transition takes place in the temperature range of interest. Because this equation applies to each medium in the reaction, the standard reaction enthalpy at temperature $T$ is:

$$\Delta_r H^\theta(T) = \Delta_r H^\theta(T_0) + \int_{T_0}^{T} \Delta_r C_p^\theta(T) dT \quad \text{Kirchhoff's law} \quad (43)$$

where $\Delta_r C_p^\theta(T)$ is the difference of the molar heat capacities of products and reactants under standard conditions. During a transformations without a change of state (liquid-solid or solid-liquid), the relative weight of the two terms ($\Delta_r^\theta H(T_0)$ and $\Delta_r C_p^\theta(T)$) in the enthalpy expression is on the order of 1 per 1000, with $H^\theta(T_0)$ being on the order of $10^7$ to $10^8$ J kg$^{-1}$ and $\Delta_r C_p^\theta(T_0)$ ranging from 10 to $100$ J kg$^{-1}$ K$^{-1}$ (considering a reaction over a temperature range of about hundred degrees Kelvin) Thiry-Muller (2018).

This observation leads to the neglect of the term $\Delta_r C_p^\theta(T_0)$ in favor of the term $\Delta_r^\theta H(T_0)$. This results in the formulation of a widely used approximation, known as the Ellingham approximation, which reads as follows, for any $T$:

$$\int_{T_0}^{T} \Delta_r C_p^\theta(T) dT \ll \Delta_r H^\theta(T_0) \quad (44)$$

On the basis of the relationship between temperature and degradation, the variation of $C_p$ as a function of thermal degradation can be identified here by the variation of $C_p$ as a function of temperature: $\Delta_r C_p^\theta(T) = \frac{\partial C}{\partial \alpha}$.

By linking this justification and demonstration to the model established throughout this study, an analogy can be drawn between the variation of $C_p$ with temperature
in this appendix and the variation of $C_p$ with the thermal degradation parameter shown in Equation (23). This relationship is expressed as :

$$ \Delta_r C_p^0(T) \approx \frac{\partial C}{\partial \alpha} \qquad (45) $$

This reasoning introduced in this appendix provides justification for the hypothesis presented in 2.2.2, leading to the Equation 24.