Xavier Colin
email: xavier.colin@ensam.eu

Fatima Essatbi

Justine Delozanne

Gurvan Moreau

A new analytical model for predicting the thermal oxidation kinetics of composite organic m atrices. Application to diamine cross-linked epoxy

Keywords:

kinetic model Carbonyl build-up Accelerated aging conditions

Introduction

The vast majority of composite material structures used in the civil aeronautical sector are made of epoxy or polyimide matrix reinforced with carbon fibers. Since the early 1980s, numerous studies clearly showed that, during their in-service use just below their glass transition temperature, these composite materials can perish by matrix embrittlement induced by thermal oxidation. This latter is restricted to a superficial layer owing to its kinetic control by oxygen diffusion [ 22 , 23 ]. "Spontaneous" cracking can then occur in this superficial layer (without any external loading) due to the development of tensile stress gradients and the catastrophic fall of fracture properties [ 16 , 21 , 24 ], thus allowing oxygen penetration into deeper layers [ 16 , 25 , 26 ]. The repetition of this sequential scenario allows damage to propagate until the core of the composite material [START_REF] Colin | Interaction between cracking and oxidation in organic matrix composites[END_REF] . Thenceforth, failure can occur untimely.

The first attempts at kinetic modeling of the thermal oxidation process in composite materials were initiated in the early 1980s by Nelson [START_REF] Nelson | Thermal ageing of graphite/polyimide composites[END_REF] , and continued for only a few years by Bowles and coll. [ 28 , 29 ]. For these authors, oxidation was a purely surface phenomenon leading to the formation of organic volatile compounds according to a one-step polymer/oxygen reaction, of which the apparent reaction order had to be determined experimentally, in particular from weight loss measurements. In the mid-1990s, Seferis and coll. [ 6 , 7 , 30 , 31 ] perfected this first kinetic model by introducing differential equations derived from solid/gas reactions in order to describe the propagation of the oxidation front from the external surfaces towards the core of the composite material. Thus, for these authors, thermal ageing transformed the composite material into a sandwich structure consisting in a completely oxidized superficial layer surrounding a still intact core. The resulting weight losses were expressed as a function of the square root of time if the thermal oxidation kinetics was diffusion-controlled, otherwise they were written as a linear function of time. However, this approach was rapidly abandoned because it led to inconsistent results. As an example, it was several times reported that the weight loss measurements did not obey the two behavioral laws mentioned above. In addition, this approach did not give access to the oxidation profiles in the sample thickness. It was thus impossible to precisely locate, in the sample thickness, the critical oxidation events responsible for its embrittlement and subsequently, to fully apply the methodology currently used for lifetime prediction.

Undoubtedly, the most realistic approach would have been to write a real coupling between oxygen diffusion and its chemical consumption (by the oxidation reaction), as already done since the 1950s for unfilled hydrocarbon polymers [START_REF] Gillen | Techniques for monitoring heterogeneous oxidation of polymers[END_REF][START_REF] Audouin | Review: Role of oxygen diffusion in polymer ageing: Kinetic and mechanical aspects[END_REF][START_REF] Tobolsky | Low temperature oxidation of hydrocarbons. The phenomenon of maximum rates[END_REF][START_REF] Boss | Oxygen diffusion limitation in autoxidation of polypropylene[END_REF][START_REF] Furneaux | Photo-oxidation of thick polymer samples. Part I: The variation of photo-oxidation with depth in naturally and artificially weathered low density polyethylene[END_REF][START_REF] Seguchi | Radiation induced oxidative degradation of polymers. Part I: Oxidation region in polymeric films irradiated in oxygen under pressure[END_REF][START_REF] Cunliffe | Photo-oxidation of thick polymer samples[END_REF][START_REF] Kiryushkin | Diffusion-controlled polymer oxidation[END_REF][START_REF] Jansen | Analytical approximation of degradation profiles in polymer products[END_REF]. In the case of unidirectional diffusion (e.g. in the sample thickness), this coupling can be written as follows:

∂C ∂t = D ∂ 2 C ∂ z 2 -r ( C ) (1)
where C is the local oxygen concentration, r(C) is the local rate of the oxygen chemical consumption (commonly called the "oxidation rate"), D is the coefficient of oxygen diffusion into the sample thickness, t is the exposure time and z is the spatial co-ordinate (i.e. the depth beneath the sample surface).

From this point of view, the kinetic model developed in the late 1990s by McManus and coll. [START_REF] Cunningham | Coupled diffusion-reaction model for predicting the distribution of degradation in polymer matrix composites[END_REF][START_REF] Mcmanus | Materials and mechanics analyses of durability tests for high-temperature polymer matrix composites[END_REF][START_REF] Crews | Modelling the high temperature degradation of graphite/epoxy[END_REF][START_REF] Mcmanus | Mechanism-based modeling of long-term degradation[END_REF] was much more satisfying. However, the oxidation rate was always expressed in a completely empirical form: r ( C ) = r 0 ∂a ( C ) ∂t

(2)

∂a ∂t = k ( 1 -a ) n C p ( 3
)
where r 0 is a constant of proportionality, a is the local conversion ratio of the oxidation process generally expressed as a function of the local variation of the polymer density, k is a rate constant obeying an Arrhenius law, and n and p are apparent reaction orders. Indeed, parameters r 0 , n and p were very difficult to justify from mechanistic considerations. In fact, they were purely adjustable parameters. Their values were found to considerably vary with the thermal aging conditions, but the physical meaning of these variations was not well understood.

At the turn of the last century, our research group at Arts et Métiers proposed to benefit from the advances realized during the past half century by the scientific community on the oxidation of unfilled hydrocarbon polymers [ 22 , 23 , 32-38 ]. Thanks to a careful kinetic analysis of the so-called "closed-loop" mechanistic scheme [START_REF] Audouin | Close-loop mechanistic schemes for hydrocarbon polymer oxidation[END_REF] , it was shown that a non-empirical analytical solution can be determined for the oxidation rate without making the usual simplifying assumptions [START_REF] Colin | Determination of thermal oxidation rate constants by an inverse method. Application to polyethylene[END_REF] that seriously degraded the reliability of all the analytical kinetic models established until then. In particular, the following general expression was obtained only if considering the steady-state assumption for hydroperoxides and radical species [ 45 , 46 ]:

r ( C ) = 2 r 0 βC 1 + βC 1 - βC 2 ( 1 + βC ) (4)
where r 0 and β are kinetic parameters which express as a function of the elementary rate constants and the substrate concentration. Both parameters are easily accessible from conventional laboratory techniques. Indeed, as schematized in Fig. 1 , β -1 corresponds the critical value of the oxygen concentration C C above which oxygen excess is reached, and r 0 is the maximal value of the oxidation rate reached in oxygen excess.

The reliability of this analytical solution was successfully demonstrated in a wide range of temperatures and oxygen partial pressures for both epoxy and polyimide matrices [ 45 , 46 ]. Subsequently, Equation ( 4) was introduced into Equation (1) by several authors in order to determine the oxidation profiles, but also to predict the consequences of oxidation on mechanical properties, first of all in epoxy [ 45 , 47 ] and polyimide matrices [ 46 , 4 8 , 4 9 ], then in their corresponding composite materials [START_REF] Olivier | Development of experimental, theoretical and numerical tools for studying thermo-oxidation of CFRP composites[END_REF][START_REF] Barjasteh | Thermal ageing of fiberglass/carbon-fiber hybrid composites[END_REF][START_REF] Barjasteh | Effect of filler on thermal aging of composites for next-generation power lines[END_REF][START_REF] Pochiraju | Evolution of stress and deformations in high temperature polymer matrix composites during thermo-oxidative aging[END_REF][START_REF] Pochiraju | Interaction of oxidation and damage in high temperature polymeric matrix composites[END_REF][START_REF] Tandon | Heterogeneous thermo-oxidative behavior of multidirectional laminated composites[END_REF][START_REF] Pochiraju | Modeling thermo-oxidative aging and degradation of composites[END_REF] .

It is important to remember here that, due to the starting (steady-state) assumption, Equation ( 4) is not valid at low conversion ratios of the oxidation reaction. Consequently, this equation cannot be used for predicting the thermal oxidation kinetics during the induction period, but also the auto-acceleration of the oxidation kinetics at the end of the induction period. This is not very penalizing at high temperatures (typically for T > 180 °C) at which accelerated aging tests in laboratory are usually performed, because the induction period is extremely short and therefore, these two first oxidation stages can be reasonably neglected. In contrast, the invalidity of Equation ( 4) prevents access to the effects of the very first oxidation events over relatively long time scales at lower temperatures (in particular, at the temperature of the industrial application), which then poses an serious problem for determining the lifetime of composite materials.

Aware of this weakness, the decision was taken to privilege the numerical approach (over the previous analytical approach) in order to eradicate all simplifying assumptions. Thus, the oxygen diffusion term (i.e. the Fick's second law) was directly introduced into the system of differential equations derived from the closed-loop mechanistic scheme, and this system was solved simultaneously in space (z) and time (t) using the numerical algorithms especially recommended for "stiff problems" of chemical kinetics [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF] . It was found that the numerical solution gives very satisfying results for both epoxy and polyimide matrices, but also for their corresponding composite materials [START_REF] Colin | Strategy for studying thermal oxidation of organic matrix composites[END_REF][START_REF] Lafarie-Frenot | Thermo-oxidation behaviour of composite materials at high temperatures: A review of research activities carried out within the COMEDI program[END_REF][START_REF] Cinquin | Thermo-oxidation behavior of organic matrix composite materials at high temperatures[END_REF] . Subsequently, this new approach was generalized by several authors to other families of unfilled polymers and organic composite materials, e.g. in reference [START_REF] Konica | A thermodynamically consistent chemo-mechanically coupled large deformation model for polymer oxidation[END_REF] .

Despite this recognized success, the numerical approach poses several practical problems. In addition to being relatively complex to understand for non-specialists in chemistry and chemical kinetics, it is very cumbersome to deploy in the common mechanical calculation codes. Indeed, the calculation time increases almost exponentially with the number of freedom degrees and the dimension of the numerical problem to be solved, so that the computing possibilities are rapidly exceeded [START_REF] Plantet | Résolution d'un problème raide de cinétique chimique[END_REF] . If there is no doubt that a numerical solution is essential for an advanced understanding of the thermal oxidation mechanisms and kinetics, it also appears clearly that a more rigorous analytical solution is lacking to allow mechanical researchers to take into account the durability of composite structures from the moment of their design and sizing. In addition, this analytical solution could be used to predict the development of damage in composite materials and their lifetime under service conditions.

The purpose of the present study is twice. At first, a new method for solving analytically the system of differential equations will be presented and the reliability of the corresponding analytical solution obtained for the oxidation rate will be checked nu- merically. This new solution will be used to accurately identify the thermal oxidation behavior of two epoxy-diamine matrices (EPO-DA) considered for composite structure applications at low to moderate temperatures (typically between 70 and 150 °C) in the civil aeronautical sector, but also to position them within the whole EPO-DA family from a durability point of view. Then, the practical interest of having an analytical solution will be demonstrated through the determination of the accelerated aging conditions inducing a simple acceleration of the natural aging kinetics without any deformation of the structural degradation state. Indeed, the research of accelerated aging conditions perfectly representative of the natural aging remains a completely open challenge in industry even today. Knowing them would allow significantly reducing the duration and cost of experimental campaigns currently performed in industry in order to qualify a material for a given application.

Theory

Closed-loop mechanistic scheme

In EPO-DA matrices, the main oxidation sites are oxy-methylene (-O-CH 2 -), amino-methylene ( > N-CH 2 -) and methanol groups ( > CH-OH) because they contain the most labile C-H bonds. Indeed, aliphatic C-H bonds located in α of an heteroatom (O or N) are characterized by a much lower dissociation energy (E D ≈ 376 kJ.mol -1 ) than aliphatic C-H bonds in polymethylenic sequences (E D ≈ 393 kJ.mol -1 ) and aromatic C-H bonds (E D ≈ 465 kJ.mol -1 ) [START_REF] Colin | Ageing and degradation of multiphase polymer systems[END_REF] .

If oxy-methylene (-O-CH 2 -) and amino-methylene ( > N-CH 2 -) groups can be already present in the starting epoxy monomers, in contrast, methanol groups ( > CH-OH) are formed during the polymerization reaction between the epoxy monomer and diamine hardener. All these oxidation sites are reported in Table 1 , as well as the corresponding name of the starting epoxy monomer.

In this study, no difference will be made between these three types of oxidation sites. They will be all noted PH for a sake of simplicity.

In a recent publication [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] , it was shown that the thermal oxidation of six EPO-DA matrices with very different glass transition temperatures can be accurately described using a single numerical kinetic model derived from the so-called closed-loop mechanistic scheme, provided that the impact of molecular mobility is taken into account in the Arrhenius laws of the different rate constants. This mechanistic scheme can be summarized by the following six reactions:

Initiation:

1) δPOOH → λP • + μPO 2 • (k 1 )
Propagation: 2)

P • + O 2 → PO 2 • (k 2 ) 3) PO 2 • + PH → POOH + P • (k 3 ) Terminations: 4) P • + P • → Inactive products (k 4 ) 5) P • + PO 2 • → Inactive products (k 5 ) 6) PO 2 • + PO 2 • → Inactive products + O 2 (k 6 )
where PH, POOH, P • and PO 2

• designate an oxidation site, an hydroperoxide, alkyl and peroxy radicals, respectively. δ, λ, and μ are stoichiometric coefficients and k i are reaction rate constants. This general form of mechanistic scheme allows describing a large variety of thermal oxidative agings differing essentially in their initiation step. In most cases of thermal aging, POOH decomposition is exclusively unimolecular at high temperature (i.e. δ = 1, λ = 2, μ = 0 and k 1 = k 1u ) or mainly bimolecular at low temperature ( δ = 2, λ = μ = 1, and k 1 = k 1b ), but rarely a combination of both modes. The boundary delimiting the domains in which each initiation predominates over each other displays the shape of Fig. 2 . This boundary corresponds to a critical hydroperoxide concentration [POOH] C reached when the respective rates of the two initiation modes are equal:

k 1u [ POOH ] = k 1b [ POOH ] 2 (5) i.e. [ POOH ] C = k 1u k 1b (6)
Since both rate constants obey an Arrhenius' law, it comes:

[ POOH ] C = A 0 Exp - E C RT (7) with A 0 = k 1u0 k 1b0 and E C = E 1u -E 1b
where k 1u0 , E 1u , k 1b0 and E 1b are the Arrhenius' parameters of the rate constants k 1u and k 1b , respectively.

For EPO-DA matrices, typical values for these Arrhenius' parameters are [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] : k 1u0 = 1.3 × 10 13 s -1 , E 1u = 130 kJ.mol -1 , k 1b0 = 1.9 × 10 9 L.mol -1 .s -1 and E 1b = 90 kJ.mol -1 .

Since E 1u > E 1b , then E C > 0 and [POOH] C is an increasing function of temperature. As an example, the numerical application of Equation ( 7) for EPO-DA matrices gives: [POOH] C = 1.7 × 10 -2 mol.L -1 at 100 °C, 7.9 × 10 -2 mol.L -1 at 150 °C, and 2.6 × 10 -1 mol.L -1 at 200 °C.

Different chemical (e.g. FTIR spectroscopy and iodometry) or calorimetric techniques (DSC) can be used to measure the initial hydroperoxide concentration [POOH] 0 . For weakly pre-oxidized polymers, [POOH] 0 is often close to 10 -2 mol.l -1 [ 65 , 66 ]. This value is of the same order of magnitude as [POOH] C at the two exposure temperatures under investigation (i.e. 120 and 150 °C). That is the reason why, in the present article, it will be considered that thermal oxidation is directly initiated in bimolecular mode:

1b

) 2POOH → P • + PO 2 • (k 1b )

System of differential equations

At the low conversion ratios of the thermal oxidation reaction (i.e. for [PH] ≈ [PH] 0 ≈ constant), the following system of differential equations can be derived from the closed-loop mechanistic scheme:

d [ P • ] dt = k 1b [ POOH ] 2 -k 2 C [ P • ] + k 3 [ PO 2 • ] [ PH ] -2k 4 [ P • ] 2 -k 5 [ P • ] [ PO 2 • ] (8) d [ PO 2 • ] dt = k 1b [ POOH ] 2 + k 2 C [ P • ] -k 3 [ PO 2 • ] [ PH ] -k 5 [ P • ] [ PO 2 • ] -2k 6 [ PO 2 • ] 2 (9) d [ POOH ] dt = -2k 1b [ POOH ] 2 + k 3 [ PO 2 • ] [ PH ] (10)
where the oxygen concentration C is related to the oxygen partial pressure P O2 in the exposure environment by the classical Henry's law: C = S . P O2 [START_REF] Tsotsis | Long-term thermo-oxidative aging in composite materials: Experimental methods[END_REF] and S is the coefficient of oxygen solubility into the polymer matrix.

There is no rigorous analytical solution for this system of differential equations, which can be solved numerically without any simplifying assumption. However, an approached analytical solution would be very useful because it would allow appreciating the behavioral trends of the reactive species

[P • ], [PO 2 • ], [POOH]
with key exposure parameters, in particular: temperature T and oxygen partial pressure P O2 . In addition, this analytical solution would then allow calculating key physico-chemical properties, from a practical point of view, because they can be checked experimentally, such as: oxygen consumption Q, carbonyl build-up [P = O], weight changes m / m 0 , etc. (see section 2.5 ). Analytical solving is only possible through the use of simplifying assumptions. The most common assumption is the so-called "steady state" assumption according to which the radical concentrations remain constant throughout the exposure duration. This assumption can only be justified in the steady-state regime where the decomposition of hydroperoxides equilibrates their formation. According to Equation [START_REF] Parvatareddy | Evaluation of chemical aging/oxidation in high performance composites via the Vickers micro-indentation technique[END_REF] , it can be thus written:

2k 1b [ POOH ] 2 ∞ = k 3 [ PO 2 • ] ∞ [ PH ] (12)
where [ POOH ] ∞ and [ PO 2 • ] ∞ are the steady concentrations of hydroperoxides and peroxy radicals, respectively. In summary, the application of the steady-state assumption to hydroperoxides is equivalent to considering a constant initiation rate for the thermal oxidation. Therefore, the concentrations of the two radical species should rapidly tend towards their steady values:

[P • ] ∞ and [ PO 2 • ] ∞ .
However, it is inconceivable to consider that the radical concentrations are constant during the induction period because this would mean that the radical species have absolutely no influence on the thermal oxidation kinetics. In oxygen excess, Tobolsky et al. [START_REF] Tobolsky | Low temperature oxidation of hydrocarbons. The phenomenon of maximum rates[END_REF] used an artefact to round this difficulty. They assumed that the sum of the radical concentrations (i.e.

[P • ] + [ PO 2 • ]
) is also a constant whereas each individual radical concentration can vary freely, which leads finally to:

d [ P • ] dt + d [ PO 2 • ] dt = 0 (13) i.e. 2k 1b [ POOH ] 2 = 2k 6 [ PO 2 • ] 2 (14)
i.e.

[ PO 2

• ] = k 1b k 6 1 / 2 [ POOH ] (15)
These results are physically reasonable insofar as only hydroperoxides and peroxy radicals are considered. Indeed, it was shown that the corresponding concentrations of these two reactive species increase in a pseudo-exponential manner with exposure time [ 67 , 68 ], which is in perfect agreement with the available experimental data. In contrast, the concentration changes in alkyl radicals are completely unrealistic because, according to Equation ( 13) , [P • ] would be maximum at the beginning of exposure and would decrease continuously to finally vanish when reaching the steadystate regime, while the thermal oxidation reaction would be in its auto-acceleration phase:

[ P • ] = [ PO 2 • ] ∞ -[ PO 2 • ] (16)
It seemed to us that a sound analytical kinetic model, aimed at describing accurately the thermal oxidation kinetics during the induction period and the auto-acceleration of the oxidation kinetics at the end of the induction period, and at predicting the induction time and the lifetime of EPO-DA matrices, cannot be based on the steady-state assumption. In other words, this latter must be replaced by a more physically reasonable and mathematically consistent assumption. How to choose such an assumption?

First novelty for solving the system of differential equations

In the previous section, it was shown that the Equation (15) obtained in oxygen excess can be used as a starting assumption. Indeed, the idea that, during the entire course of thermal oxidation, [ PO 2

• ] is proportional to [ POOH ] is quite realistic whatever the oxygen concentration [START_REF] Colin | The classical kinetic model for radical chain oxidation of hydrocarbon substrates initiated by bimolecular hydroperoxide decomposition[END_REF] . Consequently, it was decided to generalize this assumption to the entire range of oxygen concentrations:

[ PO 2 • ] = m [ POOH ] ( 17
)
where m is a function of the oxygen concentration C. The mathematical expression of m can be established by applying the steady-state assumption to radical species. In the general case, Equation ( 13) leads to:

2k 1b [ POOH ] 2 = 2k 4 [ P • ] 2 + 2k 5 [ P • ] [ PO 2 • ] + 2k 6 [ PO 2 • ] 2 ( 18
)
The introduction of Equation ( 12) into Equation ( 18) leads to:

k 4 [ P • ] 2 + k 5 [ P • ] [ PO 2 • ] - k 3 [ PH ] 2 [ PO 2 • ] + k 6 [ PO 2 • ] 2 = 0 (19)
whose the positive root is:

[ P • ] = k 5 [ PO 2 • ] 2k 4 -1 + 1 + ψ [ PO 2 • ] 0 [ PO 2 • ] -1 1 / 2 (20)
where is a kinetic parameter quantifying the competition between the different terminations of radical species:

ψ = 4 k 4 k 6 k 2 5 ( 21
)
and [ PO 2 • ] 0 is the steady concentration of peroxy radicals in oxygen excess [ 67 , 68 ]:

[ PO 2 • ] 0 = k 3 [ PH ] 2k 6 (22)
If, as shown by Gillen et al. [START_REF] Gillen | General solution for the basic autoxidation scheme[END_REF] , ψ << 1, then there is a more or less wide range of oxygen concentrations for which:

1 ≥ [ PO 2 • ] [ PO 2 • ] 0 > 5 ψ 1 + 5 ψ (23)
In this range of oxygen concentrations, Equation ( 20) can be simplified as follows:

[ P • ] = k 6 k 5 ( [ PO 2 • ] 0 -[ PO 2 • ] ) (24)
If introducing Equation ( 24) into Equation ( 9) and if applying the steady-state assumption to peroxy radicals, it comes finally: [START_REF] Bowles | Isothermal aging effect on PMR15[END_REF] where β -1 corresponds to the critical value of the oxygen concentration C C above which oxygen excess is reached:

[ PO 2 • ] ∞ = k 3 [ PH ] 2k 6 βC 1 + βC
β = 1 C C = k 6 k 2 k 5 k 3 [ PH ] (26)
The corresponding critical value of the oxygen partial pressure P O2C in the exposure environment is deduced from the classical Henry's law (Eq. 11):

P O2C = C C S = 1 β S = k 5 k 3 [ PH ] k 6 k 2 S (27)
Finally, Equation ( 12) leads to: [START_REF] Bowles | Thermo-oxidative stability studies of Celion 60 0 0/PMR-15 unidirectional composites[END_REF] It is now possible to integrate Equation [START_REF] Parvatareddy | Evaluation of chemical aging/oxidation in high performance composites via the Vickers micro-indentation technique[END_REF] . The introduction of Equation [START_REF] Bowles | Thermo-oxidative stability studies of Celion 60 0 0/PMR-15 unidirectional composites[END_REF] into Equation [START_REF] Parvatareddy | Evaluation of chemical aging/oxidation in high performance composites via the Vickers micro-indentation technique[END_REF]

[ POOH ] ∞ = k 3 [ PH ] 2k 1b [ PO 2 • ] ∞ 1 / 2 = k 3 [ PH ] 2 ( k 1b k 6 ) 1 / 2 βC 1 + βC 1 / 2 (28) and m = [ PO 2 • ] [ POOH ] = [ PO 2 • ] ∞ [ POOH ] ∞ = k 1b k 6 1 / 2 βC 1 + βC 1 / 2
leads to: d [ POOH ] dt = -2 k 1b [ POOH ] 2 + K [ POOH ] ( 30
)
where K is a rate constant expressed by: [START_REF] Salin | Mass transfer effects in degradation of bismaleimide matrix composite[END_REF] Equation ( 30) can be rewritten such as:

K = m k 3 [ PH ] = k 3 [ PH ] k 1b k 6 1 / 2 βC 1 + βC 1 / 2
d [ POOH ] dt = -K [ POOH ] 2 [ POOH ] ∞ -[ POOH ] (32)
i.e.

d [ POOH ] [ POOH ] ∞ [ POOH ] 2 [ POOH ] 2 ∞ -[ POOH ] [ POOH ] ∞ = -K dt (33)
i.e.

[ POOH ] [ POOH ]

∞ -1 - 1 [ POOH ] [ POOH ] ∞ d [ POOH ] [ POOH ] ∞ = -K dt ( 34
)
The integration of Equation ( 34) with respect to time t leads to:

Ln 1 - [ POOH ] ∞ [ POOH ] = -K t + a ( 35
)
where a is a constant to be determined by considering the initial conditions, i.e.

[ POOH ] = [ POOH ] 0 when t = 0 . Rearranging Equation ( 35) in a more adequate form leads finally to:

[ POOH ] ∞

[ POOH ]

= 1 + b Exp ( -K t ) (36)
i.e.

[ POOH ] =

[ POOH ] ∞ 1 + b Exp ( -K t ) ( 37
) with b = [ POOH ] ∞ -[ POOH ] 0 [ POOH ] 0 (38)
For weakly pre-oxidized samples, [ POOH ] ∞

[ POOH ] 0 . Experimental measurements [ 65 , 66 ] show that: [ POOH ] ∞ > 10 [ POOH ] 0 , therefore it can be reasonably considered that: b >> 1. In a first approach, a value of b = 10 was chosen in this study.

Then, the introduction of Equation (37) into Equation [START_REF] Bowles | Thermo-oxidative stability studies of Celion 60 0 0/PMR-15 unidirectional composites[END_REF] gives:

[ PO 2 • ] = [ PO 2 • ] ∞ 1 + b Exp ( -K t ) (39)
Similarly, the introduction of Equation (37) into Equation [START_REF] Colin | Mechanisms and kinetics of organic matrix thermal oxidation, In: Long-term durability of polymer matrix composites[END_REF] gives:

[ P • ] = k 3 [ PH ] 2k 5 1 - βC 1 + βC 1 1 + b Exp ( -K t ) (40)
Undoubtedly, the proportionality assumption between [ PO 2

• ] and [ POOH ] leads to sounder analytical solutions for [ PO 2 • ] and

[ POOH ] which are valid in almost the entire range of oxygen concentrations, except at the very low oxygen concentrations not satisfying Equation ( 23) . This validity range can be rewritten such as:

βC > 5 ψ (41)
In contrast, this analytical method used for solving the system of differential equations is not completely satisfying because Equation [START_REF] Mcmanus | Materials and mechanics analyses of durability tests for high-temperature polymer matrix composites[END_REF] , obtained by applying the steady-state assumption to radical species ( Eq. 24 ), gives always unrealistic changes of [P • ] with exposure time. As in Equation ( 16) , [P • ] decreases with exposure time, but fortunately here, this concentration remains very low and does not vanish when reaching the steady-state regime:

[ P • ] ∞ = k 3 [ PH ] 2k 5 1 1 + βC ( 42
)
Despite this undeniable improvement, it was decided to look for a more realistic equation.

Table 2

Result of the analytical solving of the system of differential equations without using the classical simplifying assumptions. Initial and steady values.

Species

Analytical solution Initial value Steady value General case Steady value Oxygen excess

[ POOH ] [ POOH ] ∞ 1+ b Exp ( -K t ) [ POOH ] 0 k 3 [ PH] 2 ( k 1b k 6 ) 1 / 2 ( βC 1+ βC ) 1 / 2 k 3 [ PH] 2 ( k 1b k 6 ) 1 / 2 [ PO 2 • ] [ PO 2 • ] ∞ 1+b Exp (-K t) [ PO 2 • ] ∞ [ POOH ] 0 [ POOH ] ∞ k 3 [ PH] 2 k 6 βC 1+ βC k 3 [ PH] 2 k 6 [ P ] [P • ] ∞ ( 1 1+b Exp (-K t) ) 2 [P • ] ∞ [ POOH ] 2 0 [ POOH ] 2 ∞ k 3 [ PH] 2 k 5 1 1+ βC k 3 [ PH] 2 k 5 1 βC

Second novelty for solving the system of differential equations

As in the previous section, it seemed to us interesting to start from the observation that, during the entire course of thermal oxidation, [P • ] is proportional to [ POOH ] 2 [START_REF] Colin | The classical kinetic model for radical chain oxidation of hydrocarbon substrates initiated by bimolecular hydroperoxide decomposition[END_REF] . From the Equations ( 28) and ( 42) established in the steady-state regime, it can be written:

m ' = [ P • ] ∞ [ POOH ] 2 ∞ = 2k 1b k 2 C ( 43
)
It was decided to generalize this assumption to the entire course of thermal oxidation:

[ P • ] = m ' [ POOH ] 2 = [ P • ] ∞ 1 1 + b Exp ( -K t ) 2 (44)
Undoubtedly, the proportionality assumption between [P • ] and [ POOH ] allows avoiding the last contradictions induced by the application of the steady-state assumption to radical species. Indeed, in Equation ( 44) , [P • ] increases with exposure time up to its maximal value in the regime. It should be recalled here that this sounder analytical solution for [P • ] is also valid in almost the entire range of oxygen concentrations, except at very low oxygen concentrations (see Eq. 41 ).

All these new analytical solutions, obtained for [ POOH ] , [ PO 2 • ] and [P • ] , have been compiled in Table 2 , as well as their initial and steady values (both in the general case and in oxygen excess).

Calculation of the physico-chemical properties

From these analytical solutions, it is now possible to calculate several key physico-chemical properties, from a practical point of view, because these properties can be checked experimentally. Undoubtedly, hydroperoxide concentration [POOH] and oxygen consumption Q are the most relevant properties because their mathematical expressions can be derived from the closed-loop mechanistic scheme without using any additional adjustable parameter. Thus, they can be used to accurately determine the different rate constants from experiments.

Recall that the rate of oxygen consumption writes:

r ( C ) = - d [ O 2 ] dt = k 2 C [ P • ] -k 6 [ PO 2 • ] 2 (45)
i.e.

r ( C ) = 2k 1b [ POOH ] 2 -k 6 [ PO 2 • ] 2 (46)

If replacing [ POOH ] and [ PO 2

• ] by their analytical solutions in Equation ( 46) , it comes:

r ( C ) = k 2 3 [ PH ] 2 2 k 6 βC 1 + βC 1 - βC 2 ( 1 + βC ) 1 1 + b Exp ( -K t ) 2 (47) i.e. r ( C ) = r ∞ 1 1 + b Exp ( -K t ) 2 (48) with r ∞ = 2 r 0 βC 1 + βC 1 - βC 2 ( 1 + βC ) (49) and r 0 = k 2 3 [ PH ] 2 4 k 6 ( 50
)
where r ∞ and r 0 are the respective steady rates of oxygen consumption in the general case and in oxygen excess.

It should be noted that the rate constant K ( Eq. 31 ) can be also expressed such as:

K = 2 ( r 0 k 1b ) 1 / 2 βC 1 + βC 1 / 2 (51)
The integration of Equation ( 48) with respect to time t gives access to the concentration of oxygen consumed by the chemical reaction Q:

Q = t ∫ 0 r ( C ) dt (52) i.e. Q = r ∞ t ∫ 0 dt [ 1 + b Exp ( -K t ) ] 2 (53) i.e. Q = r ∞ t ∫ 0 dt - t ∫ 0 b Exp ( -K t ) dt 1 + b Exp ( -K t ) - t ∫ 0 b Exp ( -K t ) dt ( 1 + b Exp ( -K t ) ) 2 (54) i.e. Q = r ∞ K Kt + Ln 1 + b Exp ( -K t ) 1 + b - 1 1 + b Exp ( -K t ) + 1 1 + b (55)
This curve admits as asymptotic straight-line Q ∞ in the steadystate regime (i.e. when t → ∞ ):

Q ∞ = r ∞ K Kt -Ln ( 1 + b ) - b 1 + b ( 56
)
It is now possible to deduce the mathematical expression of the induction time corresponding to the intersection point of this straight-line with the time axis:

t i = 1 K Ln ( 1 + b ) + b 1 + b ( 57
)
As b >> 1, a correct approximation is:

t i ≈ 1 -Ln ( Y 0 ) K ( 58
) with Y 0 = [ POOH ] 0 [ POOH ] ∞ (59)
It should be recalled here that K in a hyperbolic function of oxygen concentration Eq. 51 ). Thus, Equations ( 47) and ( (58) allows highlighting a very important feature of the thermal oxidation kinetics, already put in evidence on many unfilled polymers in the literature, e.g. in PP [START_REF] François-Heude | Impact of oxygen transport properties on polypropylene thermal oxidation. Part 1: Effect of oxygen solubility[END_REF]: Oxidation kinetics does not only depends on temperature, but also on the oxygen partial pressure in the exposure environment. Indeed, increasing the oxygen concentration C within the sample not only reduces the induction time ( Eq. 58 ), but also increases the oxidation rate ( Eq. 47 ). Consequently, oxygen partial pressure can be used in combination with temperature to accelerate the thermal oxidation kinetics.

Another key property, largely used in the literature to estimate the orders of magnitude of the different rate constants, is the carbonyl concentration [P = O]. In the closed-loop mechanism, these oxidation products can only be formed in initiation (1) and termination reactions (6) through specific chemical events (e.g. β scission, disproportionation, etc.) generally in competition with many other chemical events, in particular with hydrogen abstraction which gives alcohols. Consequently, the calculation of the formation rate of these oxidation products requires the use of adjustable parameters:

d [ P = O ] dt = γ 1 CO k 1b [ POOH ] 2 + γ 6 CO k 6 [ PO 2 • ] 2 (60)
where γ 1 CO and γ 6 CO are the respective formation yields of carbonyl products in initiation (1) and termination [START_REF] Nam | Anisotropic thermo-oxidative stability of carbon fiber reinforced polymeric composites[END_REF].

In a recent publication [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] , it has been shown that at 120 and 150 °C for the two EPO-DA matrices under study, γ 1 CO = γ 6 CO = 0 . 30 ± 0 . 05 . This average value will be kept in this study.

If replacing [ POOH ] and [ PO 2

• ] by their analytical solutions in Equation ( 60) , it comes:

d [ P = O ] dt = k 2 3 [ PH ] 2 4 k 6 βC 1 + βC γ 1 CO + γ 6 CO βC 1 + βC × 1 1 + b Exp ( -K t ) 2 (61)
i.e.

d [ P

= O ] dt = r CO ∞ 1 1 + b Exp ( -K t ) 2 (62) with r CO ∞ = r 0 βC 1 + βC γ 1 CO + γ 6 CO βC 1 + βC (63)
The integration of Equation ( 62) with respect to time t gives a mathematical expression similar in form to Q, thus admitting the same induction time t i :

[ P = O ] = r CO ∞ K Kt + Ln 1 + b Exp ( -K t ) 1 + b - 1 1 + b Exp ( -K t ) + 1 1 + b ( 64
)

Experimental

Materials

The two EPO-DA matrices under study result from the reaction of a common bi-functional epoxy monomer: diglycidyl ether of bisphenol F (DGEBF) or diglycidyl ether of bisphenol A (DGEBA), with an aromatic diamine hardener: 9,9-bis(3-chloro-4aminophenyl)fluorene (CAF). Films of unfilled matrices, with thicknesses typically ranging between 25 and 100 μm, were produced by compression molding then post-cured under primary vacuum (i.e. 10 -3 bar) in accordance with the recommended industrial cure cycle, with the objective of reaching the maximum crosslinking density while avoiding any undesired pre-oxidation before exposure to thermal ageing. These precautions allowed minimizing the concentration in structural defects (in particular, in hydroperoxides POOH) in both EPO-DA matrices. That is the reason why, as explained in previous in section 2.3 (see Equation ( 38) ), a high value was chosen for parameter b in this study (typically b = 10).

After processing, the films were characterized by conventional laboratory techniques. In particular, their glass transition temperature (T g ) was determined by differential scanning calorimetry (DSC, with a 20 °C.min -1 heating rate, under nitrogen) and confirmed by mechanical spectrometry (DMA, in tensile mode, with a 1 Hz frequency and a 2 °C.min -1 heating rate, under nitrogen). These measurements were found very close to the theoretical values of T g determined with the Di Marzo's equation for perfect networks (i.e. without dangling chains) [START_REF] Terekhina | Oxidation induced changes in viscoelastic properties of a thermostable epoxy matrix[END_REF] : T g ≈ 158 °C for DGEBF-CAF and 182 °C for DGEBA-CAF.

In addition, the concentration in oxidation sites [PH] was estimated from the knowledge of the repetitive monomer unit of perfect network. These concentrations were found to be roughly the same for both matrices: [PH] ≈ 14.4 mol.L -1 for DGEBF-CAF and 13.7 mol.L -1 for DGEBA-CAF. All these physico-chemical characteristics have been compiled in Table 3 for the two EPO-DA matrices under study.

Experimental procedures

The thermal oxidation kinetics of DGEBF-CAF and DGEBA-CAF matrices was studied at 120 and 150 °C under an oxygen partial pressure ranged between 0.21 bar (in air-ventilated ovens) and 10 bars (in autoclaves) for a minimum duration of 3 months. All the films were periodically removed from the ageing chambers and cooled to room temperature in a desiccator containing silica gel for preventing any moisture recovery prior to be characterized by FTIR spectrophotometry. The FTIR spectra were recorded in a transmission mode between 400 and 4000 cm -1 with a Perkin Elmer Frontier apparatus, after having averaged the 16 scans obtained with a minimum resolution of 4 cm -1 . As already shown in a recent publication [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] , the main structural modifications were observed in the carbonyl region where two new wide absorption bands appeared and grew rapidly with exposure time: one centered around 1680-1690 cm -1 , and the other around 1720-1730 cm -1 . As an example, Fig. 3 shows the changes over time in the FTIR spectrum of DGEBF-CAF matrix at 150 °C under 0.21 bar of oxygen (i.e. ambient air). These two bands were assigned to amides and other types of carbonyl products, respectively. Unfortunately, the great variety of these latter products (aldehydes, carboxylic acids, phenyl formates, etc.), often resulting from oxidation induced chain scissions in the hydroxyl propyl ether segment [START_REF] Bellenger | Oxidative skeleton breaking in epoxy-amine networks[END_REF] , did not allow identifying them precisely. Their average concentration throughout the film thickness [P = O] was determined by applying the classical Beer-Lambert's law:

[ P = O ] = OD ep ε ( 65
)
where OD is the optical density of the IR absorption band centered at 1720-1730 cm -1 (dimensionless), ε is the corresponding molar extinction coefficient (expressed in L.mol -1 .cm -1 ), and ep is the film thickness (in cm).

For carbonyl products, typical values of ε are ranged between 150 L.mol -1 .cm -1 (for ketones) and 850 L.mol -1 .cm -1 (for carboxylic acids) [ 65 , 73-76 ]. In a first approximation, an average value of 500 L.mol -1 .cm -1 was chosen for ε in the present study.

Checking the validity of the analytical solution

The numerical model recently developed for predicting the thermal oxidation kinetics of the whole family of EPO-DA matrices [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] was naturally chosen to check the reliability of the analytical Equation ( 64) in this study. Examples of numerical simulations of experimental data between 0.21 and 10 bars of oxygen at 150 °C for DGEBF-CAF and at 120 °C for DGEBA-CAF are recalled in Figures 4 and5 , respectively. It can be clearly seen that the numerical model is capable of accurately accounting for the four main successive stages commonly observed throughout the course of thermal oxidation, i.e.: the period of induction, the sharp auto-acceleration of oxidation kinetics at the end of the induction period, the steadystate regime and finally, the sudden slow-down of oxidation kinetics when the concentration of oxidation sites vanishing.

The same simulations were performed with analytical Equation [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] . Examples of comparisons between the simulations with Equation ( 64) and the numerical model at 150 °C for DGEBF-CAF and at 120 °C for DGEBA-CAF are shown in Figures 6 and7 , respectively. It can be clearly seen that Equation ( 64) is also capable of accurately accounting for all stages of the ther- mal oxidation kinetics, except the last one. Indeed, recall that Equation (64) was established by assuming low conversion ratios (i.e. [PH] ≈ [PH] 0 ≈ constant). For this reason, it cannot account for the sudden slow-down of oxidation kinetics at long-term. This behavioral deviation between the analytical and numerical models appears when the concentration in oxidation sites typically decreases by about ten percent, and it progressively increases with exposure time. However, since the embrittlement of EPO-DA matrices is generally observed at low conversion ratios, i.e. during the induction period or the auto-acceleration of the oxidation

γ 6 CO (%) 0.3 0.3 0.3 0.3 r 0 (mol.L -1 .s -1 )
4.5 × 10 -5 2.5 × 10 -4 6.0 × 10 -5 6.5 × 10 -4

β (L.mol -1 ) 13 40 15 18 k 1b (L.mol -1 .s -1 ) 2.0 × 10 -3 2.0 × 10 -2 2.0 × 10 -3 8.0 × 10 -3

Table 5

Critical values of oxygen concentration and oxygen partial pressure in DGEBF-CAF and DGEBA-CAF matrices at 150 and 120 °C.

Matrix DGEBF-CAF DGEBA-CAF T ( °C) 120 150 120 150 C C (mol.L -1 ) 2.3 × 10 -1 7.5 × 10 -2 2.0 × 10 -1 1.7 × 10 -1 P O2C (bar)
15.9 5.2 13.8 11.5

kinetics, the domain of validity of the analytical model (and thus, Eq. 64 ) is largely sufficient in order to apply the current methodology for lifetime prediction thereafter. From these comparisons, it can be thus concluded that the reliability of Equation ( 64) is demonstrated.

The values of the seven adjustable parameters used for simulating the carbonyl build-up with Equation ( 64) in Figures 6 and7 are reported in Table 4 . Do not forget that four of these parameters were already fixed in previous sections from theoretical considerations (b), or determined in a previous publication [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] (S, γ 1 CO and γ 6 CO ). Finally, the simulations allowed only identifying the three remaining parameters characterizing the thermal oxidation behavior of the two EPO-DA matrices: r 0 , β and K (i.e. k 1b ). Their values determined at 120 and 150 °C call for the following comments: i) k 1b takes almost the same values as those recently determined with the numerical model in reference [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] . It is thus confirmed that k 1b obeys an Arrhenius law with the same values of pre-exponential factor and activation energy given just below Equation [START_REF] Salin | Anisotropic effects in thermogravimetry of polymeric composites[END_REF] .

ii) In contrast, β and r 0 show much more complicated variations with temperature because these two parameters depend on rate constants k 3 and k 6 (see Eq. ( 26) and Eq. ( 50) ) which, as shown in reference [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] , are affected by the molecular mobility and, for this reason, no longer obey an Arrhenius' law in the action zone of the glass transition. Of course, the impact of the molecular mobility on β and r 0 is expected to be more pronounced for the DGEBF-CAF matrix because its glass transition temperature (T g = 158 °C) is much closer to the two exposure temperatures under investigation. This is exactly what can be seen in Table 4: When increasing the temperature from 120 to 150 °C, β increases faster whereas r 0 increases slower for this epoxy matrix, in accordance with the temperature variations of k 3 and k 6 previously reported in reference [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] .

The critical values of the oxygen concentration C C and oxygen partial pressure P O2C above which oxygen excess is reached were deduced from the values of β using Equations ( 26) and ( 27) , respectively. As expected, C C and P O2C also show complicated variations with temperature but they vary in the opposite direction of β. In the glass transition zone (i.e. at 150 °C for DGEBF-CAF), it is found that P O2C is about 5 bars, i.e. a value often reported for EPO-DA networks in the literature, e.g. in reference [START_REF] Olivier | Prévision du vieillissement thermos-oxydant de composites à matrice organique dédiés à l'aéronautique: Prise en compte des couplages multiphysiques[END_REF] . In the glassy domain, however, P O2C is found to be between two and three times higher, which requires realizing thermal ageing experiments under very high oxygen partial pressures (typically 50 bars) for accessing the thermal oxidation kinetics in oxygen excess.

Using the analytical solution for determining the accelerated aging conditions

The analytical kinetic model developed in previous sections can be easily implemented into common mechanical calculation codes for determining the consequences of oxidation on mechanical properties. In addition, this model can be used to determine the accelerated aging conditions inducing a simple acceleration of the natural aging kinetics without any deformation of the structural degradation state.

To accelerate the thermal aging, it is usual to increase the exposure temperature T. In general, care is taken not to change of physical state (e.g. not to pass from glassy to rubbery state) because the literature reports examples of the impact of the molecular mobility on the oxidation kinetics [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] . However, as repeatedly shown in the literature, e.g. in references [ 16 , 45 ], the rise in temperature causes a reduction in the thickness of the superficial oxidized layer (TOL), which will affect the global mechanical properties of the sample, but can also completely modify its mode of rupture when it is mechanically stressed. Fortunately, it is possible to play on another environmental factor to keep TOL constant: the oxygen partial pressure P O2 .

According to Audouin et al. [START_REF] Audouin | Review: Role of oxygen diffusion in polymer ageing: Kinetic and mechanical aspects[END_REF] , the steady value of TOL can be estimated by using a simple scaling law:

TOL ≈ D C r ∞ 1 / 2 (66)
where C and r ∞ are the oxygen concentration (Eq. 11) and the steady value of oxidation rate ( Eq. 49 ) at the sample surface, respectively.

As already shown by Audouin et al. [START_REF] Audouin | Review: Role of oxygen diffusion in polymer ageing: Kinetic and mechanical aspects[END_REF] , two extreme kinetic regimes can be distinguished (see Fig. 1 ): i) When the oxygen concentration C is much higher than its critical value C C above which oxygen is in excess (i.e. for βC 1 ), r ∞ ≈ r 0 , thus Equation [START_REF] Huang | Thermal oxidation of poly(dicyclopentadiene). Decomposition of hydroperoxides[END_REF]

Remember that, in a given physical state (i.e. in rubbery or glassy state), r 0 and D obey an Arrhenius law. Their respective activation energies will be noted E r and E D thereafter. In contrast, S erratically varies with temperature, so that it is generally considered independent of temperature [START_REF] Van Krevelen | Properties of Polymers: Their correlation with chemical structure; their numerical estimation and prediction from additive group contributions[END_REF] . Finally, the introduction of these temperature dependences into Equation (68) leads to: with

TOL ∝ P 1 / 2 O2 Exp - E TOL RT (69)
E TOL = 1 2 ( E D -E r ) ( 70
)
It should be pointed out that Equation ( 69) presents a real physical meaning as it allows accounting for the general trends of the diffusion-limited oxidation effects experimentally observed for many polymers by several authors. First, it well predicts that TOL increases according to the square root of the oxygen partial pressure P O2 , as observed for instance for PE and its EPR copolymer [ 23 , 35 , 38 ]. In addition, as in most cases: E D < E r , i.e. E TOL < 0 , it also predicts that TOL is a decreasing function of temperature, as already reported between 150 and 200 °C in air for epoxy matrices [ 16 , 45 ]. In this latter case, it was also shown that Equation [START_REF] Gillen | General solution for the basic autoxidation scheme[END_REF] gives the correct order of magnitude of TOL. ii) In contrast, when C is much lower than C C (i.e. for βC 1 ), r ∞ ≈ 2r 0 βC , thus Equation (68) writes as follows:

TOL ≈ D 2r 0 β 1 / 2 (71)
In this case, TOL only depends on temperature. Finally, the introduction of the temperature dependences in a given physical state (i.e. in rubbery or glassy state) of r 0 , β and D into Equation (71) leads to:

TOL ∝ Exp - E TOL RT (72)
with

E TOL = 1 2 E D -E r -E β ( 73
)
Equation ( 72) (or a very similar formulation) was shown to be particularly suitable for capturing the general trends of the diffusion-limited oxidation in PP [START_REF] Kiryushkin | Diffusion-controlled polymer oxidation[END_REF] .

At intermediary oxygen concentrations, i.e. when the increase in r ∞ is strongly slowed down with the increase in C (see Fig. 1 ), r ∞ can be approached by a simple power law in each elementary interval C of the entire concentration range under investigation (see Fig. As expected, in the general case, TOL depends both on temperature and oxygen partial pressure. Finally, the introduction of the temperature dependences in a given physical state (i.e. in rubbery or glassy state) of r 0 , β, D and S into Equation ( 76) leads to:

TOL ∝ P ( 1 -n ) / 2 O2 Exp - E TOL RT ( 77
) with E TOL = 1 2 E D -E r -n E β ( 78
)
It can be noticed that Equations ( 77) and ( 78) can be considered valid over the whole range of concentrations under investigation, because they allow finding the equations of the two extreme kinetic regimes if considering the corresponding extreme values of exponent n, i.e. Equations ( 69) and ( 70) if taking n = 0 , and Equations ( 72) and ( 73) if taking n = 1 . That is the reason why this general Equation (77) will be used, in the next paragraphs, to predict the combined effects of oxygen partial pressure and temperature.

Let us now consider a natural aging, denoted "N", characterized by a couple of mild environmental conditions ( P O2N , T N ). An accelerated aging "A" will be representative of this natural aging only if TOL remains unchanged:

TO L N TO L A = 1 (79)
i.e. if the accelerated environmental conditions ( P O2A , T A ) satisfy the following equation:

P O 2 N P O 2 A 1 -n 2 Exp - E TOL R 1 T N - 1 T A = 1 (80)
i.e.

P O2 A

= P O2 N Exp - 2E TOL ( 1 -n ) R 1 T N - 1 T A (81)
Equation (81) predicts that the conditions for accelerating natural ageing can be represented by a curve of exponential shape in the ( P O2 , T ) map.

Let us now consider the cruising flight conditions of a commercial airliner: P O2 N = 0 . 21 bar (i.e. for ambient air) and T N = 70 • C . In these conditions, TOL is typically 1.2 mm for EPO-DA matrices [START_REF] Colin | A new method for predicting the thermal oxidation of thermosets. Application to an amine crosslinked epoxy[END_REF] . Let us only focus on accelerated ageing conditions allowing to significantly reduce the duration of the ageing experiments, in particular high values of P O2 for which it can be assumed that n 1 (see Fig. 8 ). In this case, Equation (81) can be simplified as follows:

P O2 A = P O2 N Exp - 2E TOL R 1 T N - 1 T A ( 82
)
An example of accelerated ageing conditions predicted with Equation (82) has been plotted in Fig. 9 for the DGEBA-CAF matrix in glassy state (T g ≈ 182 °C). The value of E D recently identified between 10 and 110 °C for the EPO-DA matrices cross-linked with an aromatic diamine hardener [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] was used for this calculation. In addition, E r was estimated from the activation energies of the rate constants k 3 and k 6 in glassy state determined for the whole family of EPO-DA matrices in the same publication [START_REF] Colin | Towards a general kinetic model for the thermal oxidation of epoxy-diamine networks. Effect of the molecular mobility around the glass transition temperature[END_REF] . Indeed, from Equation [START_REF] Olivier | Development of experimental, theoretical and numerical tools for studying thermo-oxidation of CFRP composites[END_REF] , it can shown that: E r = 2 E 3 -E 6 (82)

Finally, E TOL was deduced from Equation ( 70) . The corresponding values of E D , E r and E TOL are compiled in Table 6 . In Fig. 9 , it can be observed that the oxygen partial pressure should be considerably increased, typically up to 527 bars for DGEBA-CAF, in order to keep TOL at a constant value when the temperature is raised up to 180 °C. This result is not surprising because when βC 1 , the oxidation rate r ∞ is almost no longer dependent on P O2 . It is therefore increasingly difficult to counterbalance the effects of an increase in temperature by an increase in oxygen partial pressure. In other words, to have an effect of P O2 comparable to that of temperature on TOL, P O2 must be increased much more than temperature.

Such extreme environmental conditions are necessary to avoid a deformation of the structural degradation state and thus, to guarantee the reliability of the results of mechanical tests. Unfortunately, they may be very difficult to access with the current experimental tools available in research laboratories.

Conclusion

The system of differential equations derived from the so-called closed-loop mechanistic scheme was solved analytically without making the usual simplifying assumptions that seriously degraded the reliability of all the analytical kinetic models established until then, but only generalizing the proportionalities observed between the different steady concentrations of reactive species: [ POOH ] , [ PO 2

• ] and [P • ] , to the entire course of thermal oxidation. These proportionality assumptions were found quite realistic. They led to a sounder analytical kinetic model capable of accurately describing the three first stages of the thermal oxidation kinetics: the induction period, the auto-acceleration of the oxidation kinetics at the end of the induction period, and the steady-state regime. This new kinetic model was used to identify the thermal oxidation behavior at 120 and 150 °C in a large range of oxygen partial pressures (between 0.21 and 10 bars) of two epoxy-diamine matrices considered for composite structure applications at low to moderate temperatures in the civil aeronautical sector. It can be now easily implemented into commercial mechanical calculation codes for determining the consequences of oxidation on mechanical properties.

In addition, the new analytical model was used to predict the accelerated aging conditions inducing a simple acceleration of the natural aging kinetics without any deformation of the structural degradation state. These accelerated aging conditions correspond to a curve of exponential shape in the ( P O2 , T ) map. It was found that the oxygen partial pressure must be increased much more than the temperature to avoid any deformation of the structural degradation state in the two EPO-DA matrices under study. However, accessing such extreme environmental conditions could pose almost insurmountable problems in practice.

Fig. 1 .

 1 Fig. 1. Shape of the curve: oxidation rate versus oxygen concentration.

  Abbreviations: Diglycidyl ether of bisphenol A (DGEBA), Diglycidyl ether of bisphenol F (DGEBF), Diglycidyl ether of bisphenol D (DGEBD), Triglycidyl ether of para-amino phenol (TGPAP), Triglycidyl ether of triphenyl methane (TGTPM), Tetraglycidyl ether of methylene dianiline (TGMDA).

Fig. 2 .

 2 Fig. 2. Shape of the boundary delimitating thermal oxidation initiated by bimolecular (1b) and unimolecular POOH decompositions (1u).

Fig. 3 .

 3 Fig. 3. Changes in the FTIR spectrum of DGEBF-CAF matrix during its thermal ageing at 150 °C under 0.21 bar of oxygen (i.e. ambient air). The final spectrum (in sky blue) was recorded after 156 hours of thermal exposure.

Fig. 4 .

 4 Fig. 4. Carbonyl build-up in 25 μm thick films of DGEBF-CAF matrix at 150 °C between 0.21 (ambient air) and 10 bars of oxygen. Comparison between simulations with the numerical model (solid lines) and experimental data (points) [64] .

Fig. 5 .

 5 Fig. 5. Carbonyl build-up in 25 μm thick films of DGEBA-CAF matrix at 120 °C between 0.21 (ambient air) and 10 bars of oxygen. Comparison between simulations with the numerical model (solid lines) and experimental data (points) [64] .

Fig. 6 .

 6 Fig. 6. Simulation of carbonyl build-up in the DGEBF-CAF matrix at 150 °C under 0.21 (ambient air), 3 and 10 bars of oxygen. Comparison between the analytical Equation (64) (dashed black lines) and the numerical model (solid colored lines).

Fig. 7 .

 7 Fig. 7. Simulation of carbonyl build-up in the DGEBA-CAF matrix at 120 °C under 0.21 (ambient air), 1, 5 and 10 bars of oxygen. Comparison between the analytical Equation (64) (dashed black lines) and the numerical model (solid colored lines).

Fig. 8 .

 8 Fig. 8. Decomposition of the oxidation rate (black line) into a series of power functions (colored lines). The equations of these different functions are specified nearby their domain of validity.

  8 ):r ∞ ≈ A r 0 ( βC ) n[START_REF] Domke | Oxidative structures in polyolefins: FT-IR method of quantitative determination[END_REF] with A = 0 . 85 ± 0 . 10 .Exponent n is close to unity when βC 1 , but decreases rapidly with C to finally almost cancel when βC 1 .The introduction of Equation[START_REF] Domke | Oxidative structures in polyolefins: FT-IR method of quantitative determination[END_REF] into Equation (66) leads to:

Fig. 9 .

 9 Fig. 9. Accelerated aging conditions representative of a natural aging of DGEBA-CAF matrix under 0.21 bar of oxygen at 70 °C predicted with Equation (82) .

Table 1

 1 Main oxidation sites in EPO-DA networks.

Table 3

 3 Molar mass of the repetitive monomer unit (mUCR), density (r) and concentration in oxidation sites (PH) for perfect EPO-DA networks.

	Matrix	m UCR (g.mol -1 )	ρ	[PH] (mol.L -1 )	T g ( °C)
	DGEBF-CAF	1041	1.25	14.4	158
	DGEBA-CAF	1097	1.25	13.7	182
	Abbreviations: Diglycidyl ether of bisphenol F (DGEBF), Diglycidyl ether
	of bisphenol A (DGEBA) and 9,9-bis(3-chloro-4-aminophenyl)fluorene
	(CAF).

Table 4

 4 Values of parameters used for simulating the carbonyl build-up in DGEBF-CAF and DGEBA-CAF matrices at 120 and 150 °C with Equation (64) .

	Matrix	DGEBF-CAF		DGEBA-CAF
	T ( °C)	120	150	120	150
	b (mol.L -1 )	10	10	10	10
	S (mol.L -1 .Pa -1 )	1.45 × 10 -7	1.45 × 10 -7	1.45 × 10 -7	1.45 × 10 -7
	γ 1 CO (%)	0.3	0.3	0.3	0.3

Table 6

 6 Activation energies of oxygen diffusion, oxidation rate and thickness of oxidized layer for DGEBF-CAF and DGEBA-CAF matrices in glassy state. E D (kJ.mol -1 ) 20 E r (kJ.mol -1 ) 112 E TOL (kJ.mol -1 ) -46

Declaration of Competing Interest

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

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests.

CRediT authorship contribution statement


Ruben Ibanez
email: ruben.ibanez-pinillo@ensam.eu

Pierre Gilormini
email: pierre.gilormini@ensam.eu

Elias Cueto
email: ecueto@unizar.es

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Numerical experiments on unsupervised manifold learning applied to mechanical modeling of materials and structures

Keywords: Nonsupervised manifold learning, State variables, Dimensionality reduction, k-PCA, Structural analysis, Material constitutive equations

The present work aims at analyzing issues related to the data manifold dimensionality. The interest of the study is twofold: (i) first, when too many measurable variables are considered, manifold learning is expected to extract useless variables; (ii) second, and more important, the same technique, manifold learning, could be utilized for identifying the necessity of employing latent extra variables able to recover singlevalued outputs. Both aspects are discussed in the modeling of materials and structural systems by using unsupervised manifold learning strategies.

Introduction

Recently, data-driven description of materials has been gaining popularity. Many complex material behaviors resisting traditional modeling procedures, or that are too complex from a microstructural viewpoint, are approached by using data-based descriptions. Different approaches are being considered. Among them include those based exclusively on measured data, others that extract the manifolds related to data, and others that attempt to enforce thermodynamic and thermomechanical consistency. The interested reader can refer to [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Bessa | A framework for data-driven analysis of materials under uncertainty: countering the curse of dimensionality[END_REF][START_REF] Liu | Microstructural material database for self-consistent clustering analysis of elastoplastic strain softening materials[END_REF][START_REF] Gonzalez | Thermodynamically consistent data-driven computational mechanics[END_REF][START_REF] Ibanez | A manifold learning approach to datadriven computational elasticity and inelasticity[END_REF][START_REF] Latorre | What-you-prescribe-is-what-you-get orthotropic hyperelasticity[END_REF][START_REF] Ladeveze | Data-driven computation for history-dependent materials[END_REF] and the numerous references therein. This work focuses on techniques based on the use of manifolds and their associated manifold learning procedures for extracting them from the available data. In general, data involve many dimensions. Consider first a sequence of three-dimensional (3D) fields defined in a domain Ω ⊂ R 3 partitioned into D voxels. Each of these fields contains many data, one datum at each voxel. Each field can be represented as a point in a vector space of dimension D (the number of voxels), where we can presume each of the D coordinate axes as reporting the value that the field of interest takes in the associated voxel. Thus, each field becomes a point in that high-dimensional space of dimension D, R D . If important correlations exist among the different fields, these points are expected to be distributed on a low-dimensional subspace embedded in the D-dimensional space. Techniques aiming at extracting these reduced subspaces, the so-called slow manifolds, sketched in Figure 1, are key tools for manipulating data and extracting their hidden information.

Thus, data define in general slow manifolds embedded in very large vector spaces due to the significant hidden correlations among them. The number of uncorrelated explicative dimensions usually becomes much smaller than the a priori assumed dimension of the space for accommodating the data. The extraction of these slow manifolds can be successfully accomplished by using linear and nonlinear dimensionality reduction techniques such as principal component analysis (PCA) in the linear case and its nonlinear counterparts ( PCA, kernel-based PCA [k-PCA], LLE, tSNE, etc.) [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF][START_REF] Maaten | Visualizing data using t-SNE[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Kambhatla | Dimension reduction by local principal component analysis[END_REF][START_REF] Zhang | Principal manifolds and nonlinear dimensionality reduction via tangent space alignment[END_REF].

These techniques can be applied to several physical systems. Moreover, when the slow manifold is determined, the solution at any point on it can be computed very accurately from a simple interpolation of neighboring data on the manifold [START_REF] Badias | An augmented reality platform for interactive aerodynamic design and analysis[END_REF], enabling almost real-time predictions and the associated real-time decision-making.

However, extracting knowledge from data associated with an existing but hidden model requires the following:

• identifying the manifold intrinsic dimension, • discovering hidden parameters, • discarding useless parameters, and • discovering the models originating the data. These questions are addressed, illustrated, and discussed in the present work in a purely methodological manner, aiming at illustrating the key concepts that could open numerous future possibilities in the field of mechanics of materials, processes, structures, and systems.

In our previous works, we addressed problems involving thousands of dimensions [START_REF] Badias | An augmented reality platform for interactive aerodynamic design and analysis[END_REF][START_REF] Gonzalez | kPCA-based parametric solutions within the PGD framework[END_REF][START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF][START_REF] Lopez | On the effect of phase transition on the manifold dimensionality: application to the Ising model[END_REF], proving that despite the apparent richness in many cases, the embedding regards a relatively low-dimensional space. However, in most of the problems that we have treated until now, their complexity prevented fine analyses of their solutions. The present paper, which is purely methodological, considers simple problems defined in low dimensions, with known solutions being easily visualizable, for facilitating an analysis and discussion. Of course, and as proved in the works just referred to, all the methodologies apply to multidimensional settings.

Unsupervised manifold learning

Let us consider a vector y ∈ R D containing experimental or synthetic data from measurements or numerical simulation. These results are often referred to as snapshots. If they are obtained by numerical simulation, they consist of nodal values of the essential variable. Therefore, these variables will be somehow correlated and, notably, there will be a linear transformation W defining the vector ξ ∈ R d , with d < D, which contains the still unknown latent variables such that y = Wξ.

The D × d transformation matrix W, which satisfies the orthogonality condition W T W = I d , is the main ingredient of the PCA and can be computed as detailed in Appendix A from the covariance matrix associated with a number (M) of snapshots y 1 , . . . , y M , which constitute the columns of matrix Y.

While PCA works with the covariance matrix (i.e., YY T ), multidimensional scaling (MDS) works with the Gram matrix containing scalar products (i.e., S = Y T Y) as described in Appendix A.

On the other hand, the k-PCA is based on the fact that data not linearly separable in D dimensions could be linearly separated if previously projected to a space in Q > D dimensions. However, the true advantage arises from the fact that it is not necessary to write down the analytical expression of that mapping as described in Appendix A.

An illustrative structural mechanics case study

Consider first a hypothetical mechanical system consisting of a prismatic beam whose three dimensions, height, width, and length, are denoted, respectively, by h, b, and L, all of them being measurable quantities. In what follows, we consider a particular output P that constitutes also a measurable quantity (buckling critical load, etc.) assumed related to those parameters from an existing but actually hidden model even if in what follows we will consider hypothetical, and most of the time, unphysical models. Thus, we consider a set of data composed of M measures y i = {h i , b i , L i , P i }, i = 1, . . . , M, with h i , b i , and L i being randomly chosen from a uniform probability distribution in their respective intervals of existence, I h , I b , and I L , respectively, defined from the following:

   I h = [h min , h max ] I b = [b min , b max ] I L = [L min , L max ]. (2)
Without any other pre-existing knowledge, one expects the output depending on the three geometrical parameters (i.e., P = P (h, b, L)). In what follows, we consider three different scenarios.

Output depending on a single parameter

In this section, we assume a quite simple model that relates the output to a single parameter,

P = αb, α ∈ R + , (3)
where α = 1000 in the numerical tests carried out. We perform M = 1000 measures, where The measures constitute a set of M = 1000 points in R 4 on which the k-PCA is applied by using the Gaussian kernel

I h ∈ [0.1, 0.
κ(y i , y j ) = exp - y i -y j 2 2β 2 , (4)
where β = 10.

Figure 2 depicts the highest eigenvalues among M resulting from the k-PCA, those lying between the highest value λ 1 and 10 -6 λ 1 .

The slow manifold associated with ξ is represented by selecting the first three reduced coordinates (ξ 1 , ξ 2 , ξ 3 ) as shown in Figure 3, where its one-dimensional (1D) intrinsic dimension is noted. This result was expected from the considered model expressed by [START_REF] Liu | Microstructural material database for self-consistent clustering analysis of elastoplastic strain softening materials[END_REF]. The points on the slow manifold are colored depending on the values of h, b, L, and P , evidencing that b constitutes the latent variable and that the output P scales (visually) linearly with it.

The process of coloring the data points in the embedded manifold deserves some additional comments due to the fact that this is used in all the analyses reported in the present paper. Manifold learning techniques look for a low-dimensional manifold defined by data points. As soon as the slow manifold is extracted, the different data points can be mapped on it. This visualization is only possible when the number of dimensions allows a simple graphical representation (as is the case for the problems addressed in the present paper). Then, these points can be colored depending on the value of the different initial coordinates, and one expects that if there is a correlation (direct or inverse and linear or nonlinear) between the initial and reduced coordinates, the colors must exhibit a certain grading.

Even if this analysis seems quite dependent on the low dimensionality of the embedding, in higher dimensional embeddings, the analysis can be performed by using local statistics. Thus, by considering a data point in the slow manifold and its closest neighbors, a local statistical analysis can be easily performed with the standard deviation indicating the dispersion of the data (equivalent to the local dispersion of colors).

One could be slightly surprised by the nonlinearity that the manifold exhibits despite the linearity of model ( 3). This nonlinearity is an artifact of the nonlinear kernel (4) used. As the model is linear, one could expect the ability of the PCA to address the problem at hand. For this purpose, it suffices transforming the kernel into its linear counterpart, giving rise to the PCA κ(y i , y j ) = y i • y j .

(5) In this case, as expected, a single nonzero eigenvalue results, and consequently the dimension of the reduced space becomes one (i.e., ξ = ξ). The associated manifold is depicted in Figure 4. Now, we consider a slightly different model, again depending on a single variable but in a nonlinear manner, according to where again α = 1000. Figure 5 depicts the 1D slow manifold, where points are colored according to the values of the variables h and P . Here, even if the direct relation can be noted, its nonlinearity is much less evident to visualize. Finally, we consider the model

P = αh 2 , (6)
P = α L 2 , ( 7
)
where α = 1000. Figure 6 depicts the 1D slow manifold, where points are colored according to the values of the variables L and P to emphasize the inverse relation between them.

Output depending on two parameters

In this section, we consider a model involving two of the three variables, in particular,

P = α h 3 L 2 , ( 8
)
where α = 1000. Figure 7 depicts the two-dimensional (2D) slow manifold, where points are colored according to the values of variables h, b, L, and P . Here, the direct and inverse effects of h and L with respect to P can be noted as well as the fact that parameter b seems, and in fact is, useless.

Identifying hidden variables

The previous case studies revealed the ability to extract the intrinsic dimensionality of the slow manifold as well as the possibility to identify useless parameters. The present case addresses a very different situation in which the model involves the three variables in the discussed case, h, b, and L. However, only two of them were measured, namely h and L with the output P , with b remaining inaccessible.

Thus, we have

P = α bh 3 L 2 , ( 9
)
where α = 1000. Therefore, the M = 1000 collected data y i , i = 1, . . . , M, reads as

y i = {h i , L i , P i } ∈ R 3 .
Figure 8 depicts the reduced points (ξ 1 , ξ 2 , ξ 3 ), which as can be seen are distributed in a domain ω ⊂ R 3 . However, no dimensionality reduction is noted, and the embedding remains 3D. A direct consequence is that many values of the output P exist for the same values on the measured inputs h and L, related to the different values of b, which affect the output P . However, as b is not measured, its value is not considered in the data points.

Such a multivalued output does not represent any conceptual difficulty. It indicates that even if both variables participate in the output, there may be others that were not considered or those considered may be useless for explaining the output. To conclude about the pertinence of these variables with respect to the considered output, we consider coloring the reduced data depending on the h and L values (refer to Figure 9). We compare them with the one where the color scales with the output P reported in Figure 8.

Thus, one could conclude that both variables h and L are relevant for explaining the output P . If we assume that the output P should be univocally explained from a small number of variables, clearly only two variables (here h and L) are not sufficient. One extra dimension suffices for recovering a single-valued output, that is, considering the reduced points in four dimensions R 4 .

As visualizing things in four dimensions is quite a difficult task, for the sake of clarity in the exposition, in what follows, we propose addressing a simpler model involving lower dimensional spaces.

We consider the simpler model

P = αbh 3 , (10)
where α = 1000. However, the M = 1000 collected data y i , i = 1, . . . , M, only deal with h and P , (i.e.,

y i = {h i , P i } ∈ R 2 ).
Figure 10 depicts the dataset y i = {h i , P i }, where it can be noted that many values of the output P are found for the same value of the variable h.

For univocally expressing the output P , we consider again the k-PCA. We compute the reduced dataset in a 3D space by considering the first three coordinates (ξ 1 , ξ 2 , ξ 3 ), while coloring these points with respect to h, to prove that h represents an explanatory variable, or with respect to P (refer to Figure 11). Figure 12 depicts the same manifold but now colored by using the hidden variable b, which proves that it contributes to explaining the output P and constitutes the hidden variable. Even if we just proved that a latent variable exists, and that it corresponds to b, as in practice we ignore this fact, we never measure the quantity b. Furthermore, it is even possible that we ignore its existence; the proposed procedure is only informative but not constructive.

Obviously, there is no unique choice. The same behavior is obtained by coloring the reduced dataset with respect to any function b p h q , where p, q ∈ R. Thus, any measurable variable ensuring such a uniform color grading could be used as the latent variable for constructing the model. There are an infinite number of possibilities but with certain constraints; in the present case, it must involve b.

To conclude this section, we address a similar but more complex situation. We consider now the richer model ( 9) with α = 1000. The M = 1000 collected data y i , i = 1, . . . , M, only deal with the input variable h and the output P (i.e., y i = {h i , P i } ∈ R 2 ), where two variables, b and L, involved in model ( 9) remain hidden. Figure 13 depicts the dataset y i = {h i , P i }, where again it can be noted that many values of the output P are found for the same value of the variable h.

For univocally expressing the output P , we consider again the k-PCA. We compute the reduced dataset in a 3D space by considering again the first three coordinates (ξ 1 , ξ 2 , ξ 3 ) while coloring these points with respect to h to prove that h represents an explanatory variable of the output P as Figure 14 proves. Figure 15 clearly reveals that by coloring the reduced points with respect to the two variables b and L taken solely, they do not represent the unique hidden variable able to explain using h the output P . The hidden latent variable according to model ( 9) should combine b and h. Figure 16 proves that the combined parameter b/L 2 perfectly works when compared with the manifold colored with respect to the output P .

It is important to note that the latent variable must involve the term b/L 2 up to any power and eventually a multiple of any power of h.

Discussion

The previous numerical experiments allow drawing a conclusion on the ability of unsupervised manifold learning techniques for identifying useless data and the existence of latent variables.

The procedure is based on the ability of embedding high-dimensional data in a lowdimensional space, where the dimension represents an approximation of the inherent dimensionality of the data. When the expected dimensionality reduction does not occur and the k-PCA reveals univocally that an extra dimension is required to accommodate the data, this indicates the existence of a hidden latent variable.

However, when the dimensionality reduction applies, by analyzing locally the manifold (any embedded data point with its closest neighbors) with respect to every original coordinate attached to each embedded data point, one can conclude that correlations exist between the initial and reduced coordinates (color grading or low statistical variance). When the dispersion increases to a large extent, one can infer that the initial coordinate analyzed is not directly correlated (even if it could be correlated in a more complex or combined manner) as discussed.

To prove the generality of these conclusions, Section 4 addresses a more complex scenario. It deals with the constitutive modeling of materials, particularly those whose behavior depends on the deformation history of the material.

Simple mechanical behaviors

In this section, we consider some simple 1D mechanical behaviors as that depicted in Figure 17.

Figure 18 shows a set of M = 1500 strain-stress pairs y i = {ε i , σ i }, i = 1, . . . , M, related to an hypothetical 1D elastic-plastic behavior with linear hardening. The depicted points represent the final mechanical state (strain-stress) of M loading-unloading random trajectories. It can be noted that with only the information of possible mechanical states, given a strain value, there are many probable values of permissible stresses and vice versa.

There is nothing physically inconsistent in having multivalued fields, but if one is looking for single-valued stress, then an extra latent variable must be introduced.

If, as in Section 3, we apply the k-PCA on the dataset y i , i = 1, . . . , M, as expected, the 2D manifold appears embedded in the R 3 space as depicted in Figure 19. This indicates that three mechanical variables, stress and strain completed by a latent extra variable p, are mutually related (e.g., σ = σ(ε, p)), giving rise to a 2D manifold embedded in R 3 .

As previously discussed, many latent variables could be certainly considered for obtaining for example a single-valued stress value. Historically, the plastic strain ε p , illustrated in Figure 20, was widely considered as a latent variable. Coloring the reduced points ξ i depicted in Figure 19 according to the plastic strain results in the colored manifold shown in Figure 21. This validates the choice of the plastic strain as a latent variable.

Toward alternative representations of the mechanical state

An alternative representation of the mechanical state, avoiding the choice of nonevident latent variables, consists in assuming that the stress at time t depends on the whole strain history, that is, on ε(τ), τ ≤ t .

When memory effects are neglected, the mechanical state can be described at time t from the strain and stress increments and the present state of stress-strain. That is, ∀t i = i∆t , y i = {ε i , σ i , ∆σ i /∆ε i }, where ∆ε i has a given magnitude and two possible signs (positive in loading and negative in unloading). Different scenarios can be analyzed:

• Linear elastic regime without plastic strain. In this case, as expected, the slow manifold constituted by the reduced states ξ i related to the mechanical states y i , depicted in Figure 22, in the 3D space (ξ 1 , ξ 2 , ξ 3 ) becomes 1D. • Nonlinear elastic behavior (with the tangent modulus E T = E ε 2 ) without plastic strain.

In this case, the results are similar to those just discussed. In Figure 23, the nonlinear behavior can be noted. • Linear elastic regime with nonzero plastic strain (perfect plasticity-no hardening). In the present case, in the mechanical state within the elastic domain, we recover from y i a reduced slow manifold of dimension two as expected, which is depicted in Figure 24. • Linear elastic regime with nonzero plastic strain (perfect plasticity-no hardening) with activated damage. In the present case, with respect to the previously discussed scenario, we presume that the material degrades with the magnitude of the plastic strain and the tangent modulus decreases accordingly. As Figure 25 depicts, this case is quite similar to the previous scenario. However, when coloring with respect to the tangent modulus, the expected uniform degraded map is noted. • The present case study considers, in addition to the points inside the elastic domain, another set of mechanical states on the elastic domain boundary with a different variable that is the product of the plastic strain with another extra variable. The latter, which could represent strain-rate sensitivity, is assumed to take arbitrary values here due to the fact that we are more interested in methodological aspects than in physical considerations. The considered mechanical states are depicted in Figure 28. When applying the k-PCA to the set of mechanical states y i = {ε i , σ i , ∆σ i /∆ε i }, the resulting manifold remains 3D. As expected, no dimensionality reduction is accomplished as Figure 29 reveals, where many values of the elastic tangent modulus can be found for the same values of the stress and strain. To investigate the nature of this behavior, we depict in Figure 30 the elastic tangent modulus versus the plastic strain, where as expected it can be noted that the former does not depend exclusively on the latter. By applying the k-PCA to the data shown in Figure 30, which are nonseparable in two dimensions, one expects to separate them by embedding in a 3D space as previously discussed and as Figure 31 proves. Finally, Figure 32 presents the slow manifold from Figure 31 but now colored with respect to the plastic strain or the extra latent variable.

Conclusion

The present work introduced, tested, and discussed issues related to manifold dimensionality with two major purposes: (i) first, when too many measurable variables are employed, manifold learning is able to discard the useless variables; (ii) second and more important, the same technique can be employed for discovering the necessity of employing and then measuring an extra latent variable that is able to recover and ensure single-valued outputs. Both issues were analyzed and discussed in two case studies, one with respect to structural mechanics and the other with respect to path-dependent material constitutive behaviors. Indeed, the physical interpretation of the discovered latent variable could imply the introduction of other measurable variables. This topic should be deeply analyzed.

The main interest in discovering these manifolds is that for a new accessible mechanical state, the output can be inferred by a simple interpolation from its neighbors on the manifold. The main aim of the present work is the construction and analysis of these manifolds. Future works, currently in progress, will focus on their use in performing data-driven simulations.

Figure 1 .

 1 Figure 1. Multidimensional data on one-(left), two-(center), and three-dimensional (right) manifolds embedded in R D .

  2], I b ∈ [0.15, 0.2], and I L ∈ [1, 1.5].

Figure 2 .

 2 Figure 2. Eigenvalues λ ∈ [λ 1 10 -6 , λ 1 ].

Figure 3 .

 3 Figure 3. Slow manifold in the 3D space defined by the first reduced coordinates (ξ 1 , ξ 2 , ξ 3 ). Each point is colored according to the value of the coordinate h (top left), b (top right), L (bottom left), or P (bottom right).

Figure 4 .

 4 Figure 4. Slow manifold ξ when considering the PCA linear dimensionality reduction. Each point is colored according to the value of the variable b.

Figure 5 .

 5 Figure 5. Slow manifold in the 3D space defined by the first three reduced coordinates (ξ 1 , ξ 2 , ξ 3 ). Each point is colored according to the value of the coordinate h (left) or P (right).

Figure 6 .

 6 Figure 6. Slow manifold in the 3D space defined by the first three reduced coordinates (ξ 1 , ξ 2 , ξ 3 ). Each point is colored according to the value of the coordinate L (left) or P (right).

Figure 7 .

 7 Figure 7. Slow manifold in the 3D space defined by the first three reduced coordinates (ξ 1 , ξ 2 , ξ 3 ). Each point is colored according to the value of the coordinate h (top left), b (top right), L (bottom left), or P (bottom right).

Figure 8 .

 8 Figure 8. Reduced representation ξ i ∈ R 3 of the data y i ∈ R 3 , where the reduced points are colored according to the value of the output P .

Figure 9 .

 9 Figure 9. Reduced points colored according to the values of the coordinates h (left) and L (right).

Figure 10 .

 10 Figure 10. Dataset y i = {h i , P i }, i = 1, . . . , M.

Figure 11 .

 11 Figure 11. Reduced points colored according to the values of the coordinates h (left) and P (right).

Figure 12 .

 12 Figure 12. Reduced points colored according to the value of the output b.

Figure 13 .

 13 Figure 13. Dataset y i = {h i , P i }.

Figure 14 .

 14 Figure 14. Reduced points colored according to the value of the coordinate h.

Figure 15 .

 15 Figure 15. Reduced points colored according to the values of the variables b (left) and L (right).

Figure 16 .

 16 Figure 16. Reduced points colored according to the values of the variables b/L 2 (left) and P (right).

Figure 17 .

 17 Figure 17. Elastic behavior (left), elastic-perfectly plastic behavior (center), and elastoplasticity with linear hardening (right).

Figure 18 .

 18 Figure 18. Strain-stress mechanical states related to elastic-plastic behavior with linear hardening.

Figure 19 .

 19 Figure 19. Elastic-plastic manifold embedded in the 3D space.

Figure 20 .

 20 Figure 20. Plastic strain definition (left) and mechanical states y i colored according to the plastic strain (right).

Figure 21 .

 21 Figure 21. Manifold colored according to the plastic strain.

Figure 22 .

 22 Figure 22. Slow manifold related to linear elastic behavior and colored according to the tangent modulus (i.e., ∆σ i /∆ε i ).

Figure 23 .Figure 24 .

 2324 Figure 23. Slow manifold related to nonlinear elastic behavior and colored according to the tangent modulus (i.e., ∆σ i /∆ε i ).

Figure 25 .

 25 Figure 25. Slow manifold related to damageable elastic-plastic behavior operating within the elastic domain colored according to the tangent modulus.

Figure 26 .

 26 Figure 26. Slow manifold related to elastic-plastic behavior operating within the elastic domain and on the elastic domain boundary colored according to the tangent modulus.

Figure 27 .

 27 Figure 27. Slow manifold related to elastic-plastic behavior operating within the elastic domain and on the elastic domain boundary, colored according to the tangent modulus, when damage is activated.

Figure 28 .

 28 Figure28. Mechanical states within the elastic domain, colored according to the tangent modulus, the last depending on the product of two latent variables, the plastic strain, and another arbitrarily chosen variable.

Figure 29 .

 29 Figure 29. Slow manifold related to elastic-plastic behavior operating within the elastic domain, colored according to the tangent modulus, with the last scaling with the product of the plastic strain and an extra latent variable.

Figure 30 .

 30 Figure 30. Elastic tangent modulus versus the plastic strain.

Figure 31 .

 31 Figure 31. Slow manifold related to the data consisting of the elastic tangent modulus and the plastic strain.

Figure 32 .

 32 Figure 32. Slow manifold of the elastic tangent modulus colored according to the plastic strain (left) and the extra latent variable (right).

Acknowledgments

The first, third, and fourth authors are supported by their respective ESI Group research chairs; their support is gratefully acknowledged. The first author is supported by CREATE-ID ESI-ENSAM research chair. The third author is supported by the ESI Group Chair at the University of Zaragoza. The fourth author is supported by CREATE-ID ESI-ENSAM research chair.

Dedication

This manuscript was written with the help of the contributions of all authors. All authors have given approval to the final version of the manuscript.

Conflicts of interest

The authors declare no competing financial interests.

Appendix A. From principal component analysis to its kernel-based counterpart

A.1. Principal component analysis

Let us consider a vector y ∈ R D containing experimental results or synthetic data from a numerical simulation. These results are often referred to as snapshots. If they are obtained by numerical simulation, they consist of nodal values of the essential variable. Therefore, these variables will be somehow correlated and, notably, there will be a linear transformation W defining the vector ξ ∈ R d , with d < D, which contains the still unknown latent variables, such that

The D × d transformation matrix W, which satisfies the orthogonality condition W T W = I d , is the main ingredient of the PCA [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF].

Assume that there exist M different snapshots y 1 , . . . , y M , which we store in the columns of a D×M matrix Y. The associated d × M reduced matrix Ξ contains the associated vectors ξ i , i = 1, . . . , M.

The PCA usually works with centered variables. In other words,

implying the necessity of centering data before applying the PCA.

The PCA proceeds by guaranteeing maximal preserved variance and minimal correlation in the latent variable set ξ. The latent variables in ξ are therefore uncorrelated, and consequently the covariance matrix of ξ,

should be diagonal.

To extract the d uncorrelated latent variables, we proceed from

Pre-and post-multiplying by W T and W, respectively, and making use of the fact that W T W = I, give us

The covariance matrix C y y can then be factorized by applying the singular value decomposition,

where V contains the orthonormal eigenvectors; Λ is a diagonal matrix containing the eigenvalues sorted in descending order. Substituting ( 16) into [START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF], we arrive at

This equality holds when the d columns of W are taken to be collinear with d columns of V. We then preserve the eigenvectors associated with the d nonzero eigenvalues,

which gives

We therefore conclude that the eigenvalues in Λ represent the variance of the latent variables (diagonal entries of C ξξ ).

A.2. Multidimensional scaling

The PCA works with the covariance matrix of the experimental results, YY T . However, the MDS works with the Gram matrix containing scalar products (i.e., S = Y T Y) [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF].

The MDS preserves pairwise scalar products:

Computing the eigenvalues of S, we arrive at

which in turn gives

A.

Kernel-based principal component analysis

The k-PCA is based on the fact that data not linearly separable in D dimensions could be linearly separated if they are previously projected to a space in Q > D dimensions. However, the true advantage arises from the fact that it is not necessary to write down the analytical expression of that mapping.

The symmetric matrix Φ = Z T Z, with Z containing the snapshots z i ∈ R Q , i = 1, . . . , M, associated with y i ∈ R D , has to be decomposed into eigenvalues and eigenvectors. The procedure for centering data z i is carried out in an implicit way. The eigenvector decomposition reads as

giving rise to

The difficulties of operating in a high-dimensional space of dimension, in general, Q D, and the mapping unavailability are circumvented by introducing the kernel functional κ (also known as the kernel trick). This allows computing scalar products in R Q while operating in R D by applying the Mercer theorem. This theorem establishes that if κ(u, v) (where u ∈ R D and v ∈ R D ) is continuous, symmetric, and positive definite, then it defines an inner product in the mapped space R Q . Many different kernels exist; some of them are reported in [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF].