Qiang Chen

George Chatzigeorgiou

Gilles Robert

Fodil Meraghni
email: fodil.meraghni@ensam.eu

viscoelastic-viscoplastic and damage mechanisms

Simulating the nonlinear response of structures under a large number of cycles has been a bottleneck problem in the literature, especially when tens of thousands to millions of fatigue cycles are involved [START_REF] Russo | A robust cumulative damage approach for the simulation of delamination under cyclic loading conditions[END_REF][START_REF] Bian | Elasto-plastic analysis of critical fracture stress and fatigue fracture prediction[END_REF][START_REF] Qin | Inhomogeneous deformation growth of a metal under cyclic loading and its influence on fatigue[END_REF]. A simple load cycle is normally divided into many increments and iterations. This applies in particular to largescale heterogeneous systems that contain many dissipative mechanisms, where the cycle-by-cycle simulations are barely affordable. An example in point is the injection-molded short glass fiber-reinforced polyamide 66 composites (PA66/GF), frequently used for automotive components subjected to cyclic loading. These composites have a specific skin-core-shell microstructure with short fibers distributed randomly in the laminate plane during the injection molding manufacturing process [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF][START_REF] Hessman | Microstructural analysis of short glass fiber reinforced thermoplastics based on x-ray micro-computed tomography[END_REF]. Moreover, the PA66/GF composites exhibit remarkable time-dependent inelastic deformation due to the viscoelastic and viscoplastic mechanisms in the polyamide phase. Another important consideration that needs to be taken into account is the progressive matrix microcracking and interface debonding formed in this class of materials under a variety of loading and environmental conditions when the internal stress reaches the strength of the polyamide matrix or the interface [START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF][START_REF] Arif | In situ damage mechanisms investigation of PA66/GF30 composite: effect of relative humidity[END_REF].

Homogenization theories are good predictive tools for the composite's macroscopic response using as information the response of the individual phases and the microstructures [START_REF] Pindera | Micromechanics of spatially uniform heterogeneous media: a critical review and emerging approaches[END_REF][START_REF] Saeb | Aspects of computational homogenization at finite deformations: a unifying review from Reuss' to Voigt's bound[END_REF][START_REF] Chen | Homogenization and localization of nanoporous composites-a critical review and new developments[END_REF][START_REF] He | Finite volume based asymptotic homogenization of viscoelastic unidirectional composites[END_REF][START_REF] Schneider | A review of nonlinear FFT-based computational homogenization methods[END_REF]. Homogenization techniques, in general, can be categorized into two broad categories, namely the full-field (i.e., periodic) and mean-field homogenization techniques. The periodic homogenization techniques, such as finite element [START_REF] Saeb | Aspects of computational homogenization at finite deformations: a unifying review from Reuss' to Voigt's bound[END_REF][START_REF] Zhi | Direct FE2 modeling of heterogeneous materials with a micromorphic computational homogenization framework[END_REF][START_REF] Chen | Homogenization of size-dependent multiphysics behavior of nanostructured piezoelectric composites with energetic surfaces[END_REF] and finite volume-based approaches [START_REF] Chen | Finite-volume homogenization of elastic/viscoelastic periodic materials[END_REF][START_REF] He | Finite volume based asymptotic homogenization theory for periodic materials under anti-plane shear[END_REF][START_REF] Chen | Homogenization and localization of elastic-plastic nanoporous materials with Gurtin-Murdoch interfaces: an assessment of computational approaches[END_REF], are considered to be computation standards as they can take into account the detailed interaction of the microstructures. A major limitation of the periodic homogenization technique is that extensive mesh refinement is required in order to model the three-dimensional complex microstructures. For describing high-cycle fatigue behavior, the computation time of each loading cycle becomes essential for developing any micromechanics model. Therefore, the full-field homogenization methods are not feasible for high-cycle calculations of complex PA66/GF composites.

The mean-field homogenization techniques are analytical. Thus, they are extremely advantageous to periodic homogenization in terms of computation efficiency. These methods are accurate enough for the elastic moduli estimations, but appear to overestimate the stress-strain response when significant nonlinearities occur [START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: unified formulation, application and benchmark[END_REF][START_REF] Barral | Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: theory and numerical validation[END_REF][START_REF] Haddad | Viscoelastic-viscoplastic polymer composites: development and evaluation of two very dissimilar mean-field homogenization models[END_REF]. Our early attempts aimed at capturing the viscoelastic-viscoplastic deformation and the ductile damage of PA66/GF employ an extended micromechanics approach combining Mori-Tanaka with the transformation field analysis (TFA) [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF][START_REF] Chen | Extended mean-field homogenization of viscoelastic-viscoplastic polymer composites undergoing hybrid progressive degradation induced by interface debonding and matrix ductile damage[END_REF]. In this model, an interphase layer was introduced to correct the overestimation of stress-strain behavior in the nonlinear regime. Additionally, a phenomenological model consisting of four Kelvin-Voigt branches and a viscoplastic branch, formulated under the thermodynamics framework by Praud et al. [START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], was utilized to describe the local deformation mechanisms of the polyamide phase. The interface debonding was described by the Weibull probabilistic density function to simulate the initiation and coalescence of the void-type discrete damage in the vicinity of the fiber/matrix interphase. It has been clearly revealed [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF] that the extended Mori-Tanaka TFA approach for short fiber-reinforced composites with random fiber orientations, accounting for damage, viscoelastic and viscoplastic coupling mechanisms, is capable of reproducing the oligocyclic experimental response of PA66/GF composites with good accuracy. However, since a considerable amount of phases (for instance 36 orientations) needs to be considered in the extended Mori-Tanaka approach in order to properly capture the fiber orientation distribution, conducting the high-cycle analysis passing through every load cycle is still prohibited due to the high computation costs.

The cycle jump approach is a viable technique to circumvent the computational challenges posed by the cycle-by-cycle simulation. The cycle jump technique was originally proposed by Troschenko et al. [START_REF] Troshchenko | Calculation of fatigue and life of crack-bearing structural elements under cyclic load[END_REF] for the calculation of fatigue and life of crack-bearing structural elements, and further extended by Lesne and Savalle [START_REF] Lesne | An efficient cycles jump technique for viscoplastic structure calculations involving large number of cycles[END_REF] for viscoplastic structure calculations. It has been adopted by numerous researchers for accelerated highcycle computations [START_REF] Abdul-Latif | New hybrid cycle jump approach for predicting low-cycle fatigue behavior by a micromechanical model with the damage induced anisotropy concept[END_REF][START_REF] Loew | Accelerating fatigue simulations of a phase-field damage model for rubber[END_REF][START_REF] Titscher | Efficient higher-order cycle jump integration of a continuum fatigue damage model[END_REF][START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF]. For describing the high-cycle response of composite materials, Van Panpegem and coworkers [START_REF] Van Paepegem | Finite element approach for modelling fatigue damage in fibre-reinforced composite materials[END_REF][START_REF] Van Paepegem | Fatigue degradation modelling of plain woven glass/epoxy composites[END_REF] have incorporated the cycle jump scheme into the continuum damage models to extrapolate the stiffness degradation of plain woven composites. Fish et al. [START_REF] Fish | Multiscale fatigue life prediction model for heterogeneous materials[END_REF] have introduced the cycle jump technique into a two-scale asymptotic homogenization approach to study the multiscale fatigue life of high-temperature ceramic matrix composites. Luders et al. [START_REF] Lüders | Adaptive cycle jump and limits of degradation in micromechanical fatigue simulations of fibre-reinforced plastics[END_REF] have developed an adaptive cycle jump method for micromechanical fatigue simulations of fiber-reinforced plastics. Sally et al. [START_REF] Sally | An efficient computational strategy of cycle-jumps dedicated to fatigue of composite structures[END_REF] have proposed an efficient computational strategy of cycle jumps dedicated to the fatigue of composite structures. Whatever the developed algorithm is, a common assumption adopted in various cycle jump techniques is that a single cycle is very short compared to the full cycle lifetime; therefore, the change of material states in one cycle is negligible compared to the change over the entire cycle lifetime. The general idea of the cycle jump method is then to avoid conducting cycle-by-cycle analysis with only a set of loading cycles simulated. At the interval where a large number of loading cycles can be skipped, the state-dependent variables (SDVs) necessary for obtaining an appropriate homogenized behavior for any given load increments at each material phase are extrapolated to a certain jump length using the Euler forward approximation. The latter is based on the information of the simulated cycles before the attempted jumps. As a result, the computation time can be significantly reduced. Those extrapolated SDVs are subsequently utilized as initial conditions to perform a new set of control cycles after the cycle jumps.

In this study, a major interest is the implementation of the cycle jump procedure into the extended Mori-Tanaka TFA micromechanics method for the first time and employing the extended theory to predict the high-cycle response of PA66/GF composites under various loading conditions. This extension is based on performing Mori-Tanaka TFA simulation considering the viscoelastic-viscoplastic and damage mechanisms for a number of cycles to establish the global evolution functions of all the SDVs, extrapolating the SDVs, such as strain, stress, damage variables, viscoelastic and viscoplastic strains, which span over significant cycle lengths, and conducting a set of control cycles using the Mori-Tanaka TFA after the cycle jumps to obtain the updated global evolution functions. It should be noted that the determination of the proper cycle jump length is a key issue for the performance of the cycle jump technique since a high jump length reduces the computational cost, but simultaneously can increase the extrapolation error if not the numerical instability. Therefore, to enable efficient implementation of the cycle jump procedure, the extrapolation scheme with the control function proposed by Cojocaru and Karlsson [START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF] is adopted, which permits an automatic jump length and extrapolating error control over the cycle jump period. To the best knowledge of the authors, thus far, there has been no reported work of micromechanical simulations of the inelastic behavior of short glass fiber-reinforced polyamide composites with such complex behavior under high-cycle loading in the open literature. Such work is important to reduce the cost-prohibitive fatigue test of PA66/GF composites in automotive applications.

The remaining work is organized as follows: Sect. 2 recalls the theoretical framework of the extended Mori-Tanaka TFA approach with random fiber orientations and describes the implementation of the cycle jump algorithm into the extended Mori-Tanaka TFA micromechanics. Section 3 investigates the effect of cycle jump control parameters on the accuracy of the simulated response. Numerical verifications of the predictive capabilities of the cycle jump approach for simulating the inelastic response of short fiber-reinforced polyamide composites under high-cycle fatigue loading are also performed in this section. Comparisons of the extended cycle jump micromechanics simulations with the experimental results of actual composites under high-cycle loading are presented in Sect. 4. Limitations of the present work are discussed in Sect. 5. Section 6 draws the pertinent conclusions.

Theoretical development

2.1 Extended Mori-Tanaka TFA approach for randomly oriented short-fiber composites In this section, the extended Mori-Tanaka TFA approach, developed by Chen et al. [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF], accounting for the viscoelastic-viscoplastic and damage mechanisms is briefly recalled. For an N-orientation coated fiberreinforced composite, as shown in Fig. 1, the averaged phase strain can be obtained using the transformation field analysis technique proposed by Dvorak [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF] and Dvorak and Benveniste [START_REF] Dvorak | On transformation strains and uniform fields in multiphase elastic media[END_REF] as follows:

ε i r T i r : ε 0 + T i r 0 : ε in 0 + T i,i r 1 : ε in,i 1 + T i,i r 2 : ε in,i 2 , ( 1 )
where the superscripts i 1, ..., N are associated with variables in the i th fiber orientation while the subscripts r 1, 2 represent the fiber and interphase phases, respectively. ε 0 and ε in 0 denote the average matrix total strain and the matrix inelastic strain. ε in r (r 1, 2) represents the inelastic strain of the r th phase. T i r (r 1, 2) are the elastic interaction tensors of the i th fiber and i th coating, respectively. T i r p ( p 0, 1, 2) are the inelastic interaction tensors. For the expressions of these elastic and inelastic interaction tensors, the readers are referred to Chen et al. [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF].

The homogenized strain for the overall composites is obtained as the weighted sum of the averaged strains over all the phases: where c 0 is the volume content of the matrix phase and c i r (r 1, 2) are the volume contents of the ith fiber and coating. Substituting Eq. (1) into Eq. ( 2) yields the homogenized strain in the following form:

ε c 0 ε 0 + N i 1 c i 1 ε i 1 + c i 2 ε i 2 , ( 2 )
ε c 0 I + N i 1 c i 1 T i 1 + c i 2 T i 2 : ε 0 + N i 1 c i 1 T i 10 + c i 2 T i 20 : ε in 0 + N i 1 c i 1 T i,i 11 + c i 2 T i,i 21 : ε in,i 1 + N i 1 c i 1 T i,i 12 + c i 2 T i,i 22 : ε in,i 2 . ( 3
)
To obtain the homogenized constitutive relations of the overall composites, the averaged strains in the matrix, fiber and coating are expressed in terms of the homogenized strain for the overall composites and inelastic strains of each phase as follows:

ε 0 A 0 : ε + A 00 : ε in 0 + N i 1 A i 01 : ε in,i 1 + N i 1 A i 02 : ε in,i 2 , ε i 1 A i 1 : ε + A 10 : ε in 0 + N j 1 A j,i 11 : ε in, j 1 + N j 1 A j,i 12 : ε in, j 2 , ε i 2 A i 2 : ε + A 20 : ε in 0 + N j 1 A j,i 21 : ε in, j 1 + N j 1 A j,i 22 : ε in, j 2 . ( 4
)
In the above equation, the various A r and A r p matrices denote the elastic and inelastic strain concentration tensors, respectively, which can be considered as the generalization of Eshelby's problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] for single inclusion to multi-inclusions interacting together in an average sense. The expressions for those concentration tensors revolve around the Eshelby tensors [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] that are expressed in terms of the matrix and interphase secant moduli, and the geometrical parameters of the inclusions. They can be evaluated analytically in the case of infinitely long cylindrical fibers or prolate spheroids embedded in an isotropic matrix [START_REF] Mura | Mechanics of Elastic and Inelastic Solids[END_REF]; otherwise, numerical technique developed by [START_REF] Gavazzi | On the numerical evaluation of Eshelby's tensor and its application to elastoplastic fibrous composites[END_REF] should be employed.

It should be noted that just like all other mean-field micromechanics methods, the extended Mori-Tanaka TFA is based on the averaged field per phase. Therefore, the predictions of the overall inelastic behavior of a composite based on the extended Mori-Tanaka TFA method are considerably stiffer than the experimentally measured response or the response of full-field homogenization methods. The latter are considered to be computation standards in the literature. This was expected because the non-uniform stress/strain field concentration, predicted by the full-field homogenization, produces localized inelastic strains in the vicinity of the fibers. To account for the differences between the ε in 0 and ε in 2 , a correction tensor Y is introduced as follows [START_REF] Barral | Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: theory and numerical validation[END_REF]:

ε in 2 Y : ε in 0 . ( 5
)
For a short glass fiber composite with fibers aligned in x 3 direction, the correction tensor Y takes the following form:

Y ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ γ N 0 0 0 0 0 0 γ N 0 0 0 0 0 0 γ L 0 0 0 0 0 0 γ ST 0 0 0 0 0 0 γ SL 0 0 0 0 0 0 γ SL ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 6
)
where γ N , γ L , γ ST and γ SL represent the differences between the inelastic strains in the matrix and the inelastic strains in the virtual interphase in transverse normal, longitudinal, transverse shear and axial shear directions, respectively, which can be calibrated from the macroscopic experimental response or the full-field simulations. The averaged stress in each phase can be expressed as

σ i 1 C i 1 : ε i 1 -ε in,i 1 , σ i 2 C i 2 : ε i 2 -ε in,i 2 C i 2 : ε i 2 -Y : ε in 0 , σ 0 C 0 : ε 0 -ε in 0 , ( 7 )
where C 0 denotes the secant modulus of the matrix. C i 1 and C i 2 denote the secant moduli of the fiber and interphase in the i th direction, respectively. The computation of the inelastic strains and secant moduli depends on the chosen phase constitutive law, which will be given in Sect. 2.2.

Ultimately, the homogenized constitutive relation of the composites in the presence of damage and inelastic strains can be obtained as

σ c 0 σ 0 + N i 1 c i 1 σ i 1 + c i 2 σ i 2 C * ε -σ in , ( 8 )
where σ in denotes the homogenized inelastic stress. C * is the effective secant modulus of the N-orientation coated short fiber-reinforced composite given by the following relation:

C * c 0 C 0 A 0 + N i 1 c i 1 C i 1 : A i 1 + c i 2 C i 2 : A i 2 . ( 9
)

Constitutive laws for the matrix and interphase

The extended Mori-Tanaka TFA approach employs a phenomenological model consisting of four Kelvin-Voigt branches and a viscoplastic branch connected in series, formulated under the thermodynamics framework by Praud et al. [START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], to characterize the local deformation mechanisms of the polyamide matrix phase. The damage mechanism in the matrix was described by the continuum damage mechanics theory based on the wellestablished principle of effective stress proposed by Lemaitre and Chaboche [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. Under the thermodynamics framework, the Helmholtz free energy for the matrix medium can be expressed as

ρψ ε, ε vi , ε vp , p, d 1 2 ε - N i 1 ε vi -ε vp : (1 -d)C e : ε - N i 1 ε vi -ε vp + 1 2 N i 1 ε vi :(1 -d)C vi : ε vi + p 0 R(ξ ) dξ , ( 1 0 )
where C e (representing the initial matrix stiffness tensor C 0 ) and C vi denote the elastic stiffness tensor and the viscoelastic tensor of the i th branch, respectively. ε, ε vi and ε vp are total, viscoelastic and viscoplastic strains, respectively. p is the viscoplastic scalar variable, and d denotes the ductile matrix damage variable represented as a scalar quantity. R is the hardening function. The evolution laws for the interval variables (ε vi , ε vp , p and d) and their activation criteria are presented in Praud et al. [START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF]. Under the continuum damage mechanics theory, the secant modulus of the matrix in Eq. ( 7) can be expressed as:

C 0 (1 -d)C e .
Additionally, virtual interphases are assumed to exist between the fibers and the matrix. They are occupied by the same polyamide material initially as in the matrix. The average viscoelastic-viscoplastic strains within the virtual interphases are directly linked to the matrix inelastic strains through Eq. ( 5). The damage mechanism in the interphase layers, however, is remarkably different from the matrix phase. While the damage in the matrix is characterized by gradual stiffness degradation that can be described by the continuum damage mechanics theory, the interphase debonding is catastrophic with the rapid occurrence and the coalescence of voids or defects [START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF][START_REF] Arif | In situ damage mechanisms investigation of PA66/GF30 composite: effect of relative humidity[END_REF]. Therefore, the stiffness reduction D(γ c ) [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF], induced by the void creation, is introduced to mimic the interphase debonding in the virtual coating layers.

C 2 (γ c ) C 0 -D(γ c ) (1 -γ c )C 0 : A m (γ c ) (11)
In the above equation, A m (γ c ) is the strain concentration tensor of the net matrix embedding voids and γ c is the crack density [START_REF] Chen | Extended mean-field homogenization of viscoelastic-viscoplastic polymer composites undergoing hybrid progressive degradation induced by interface debonding and matrix ductile damage[END_REF][START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF]. To naturally track the crack initiation and growth, the cumulative probability of the interphase crack density is described by the following Weibull-type probabilistic density function:

γ c γ max 1 -exp - σ eff 2 σ c κ , ( 1 2 )
where γ max indicates the saturation limit of the microcrack density. σ e f f 2 denotes the homogenized interphase effective stress. σ c and κ are the Weibull parameters controlling how fast the interphase crack can grow.

Extended Mori-Tanaka TFA approach with controlled cycle jump: new developments

The extended Mori-Tanaka TFA approach with randomly oriented fibers can be implemented using the radial return mapping algorithm for the nonlinear micromechanics analysis [START_REF] Chatzigeorgiou | Computational Methods[END_REF][START_REF] Chatzigeorgiou | Chapter 12-Nonlinear composites[END_REF]. At each loading step, many loading increments and iterations are utilized to ensure the numerical convergence and stability of the nonlinear analysis. It should be noted that for an N-orientation coated fiber composite, the total number of interaction tensors is 8N , while the total number of the concentration tensors is 4N 2 + 6N + 2. These tensors need to be recalculated at each iteration of every loading increment. For describing high-cycle fatigue behavior, which is of particular interest in this work, the computation time becomes crucial for every single cycle simulation. The determination of these tensors in the Mori-Tanaka TFA approach with consideration of damage and various inelastic mechanisms relies on the numerical integration technique, hence consuming a considerable amount of computation resources and time.

To alleviate this issue, an accelerated Mori-Tanaka TFA scheme is introduced herein by assuming that there are no drastic variations in the state-dependent variables between the consecutive loading cycles under high-cycle loading. The general goal of the accelerated numerical scheme is to avoid simulating all individual cycles in a high-cycle analysis through the use of the "cycle jump" concept, therefore reducing significantly the need for repetitive evaluation of the interaction and concentration tensors. As illustrated in Fig. 2, the cycle jump procedure consists of the following steps:

Training cycle: a number of training cycles are conducted using the Mori-Tanaka TFA approach to obtain the global evolution functions of all the SDV s, which include the macroscopic stress, strain, inelastic strain, as well as viscoelastic strain, viscoplastic strain, damage variable in each phase; Cycle jump: establishing the maximum common jump length of all the SDVs and extrapolating the SDV s based on the global evolution function to a certain jump length; Control cycle: the extrapolated SDVs are used as the initial state for a new set of the Mori-Tanaka TFA micromechanics simulations after the cycle jump. The global evolution functions of all the SDVs are updated according to the control cycle analysis; Repeat steps 2 and 3 to obtain the response till the desirable number of loading cycles or time.

Hereafter, the allowable jumped cycle in the extended Mori-Tanaka TFA analysis is determined based on the extrapolation scheme proposed by Cojocaru and Karlsson [START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF]. The latter permits automatic determination of the jump length during the cyclic loading. When a fast variation in the SDVs is observed, the extrapolation scheme automatically performs shorter or no jumps. Otherwise, a longer cycle jump will be conducted while the accelerated solutions remain close to the cycle-by-cycle simulation within a given tolerance. Let us assume that the solutions for at least three successive loading cycles have been obtained using the Mori-Tanaka TFA computation, as shown in Fig. 3. To establish the global evolution functions of state variables, the values for each SDV are extracted from the same relative position in time within each loading cycle, which are denoted by S DV (t 1 ), S DV (t 2 ) and S DV (t 3 ). The differences in the state variables between the last two cycles are:

S DV (t 1 ) S DV (t 1 ) -S DV (t 2 ), (13)
S DV (t 2 ) S DV (t 2 ) -S DV (t 3 ).

The corresponding slopes are expressed as:

s 12 (t 1 ) S DV (t 1 ) t cycle , ( 15
)
s 23 (t 2 ) S DV (t 2 ) t cycle , ( 16
)
where t cycle t 1t 2 t 2t 3 represents the period of one cycle. The allowable jump length t jump for the designated SDV is given by the following criterion:

s p t 1 + t jump -s 12 (t 1 ) |s 12 (t 1 )| ≤ q, ( 17
)
where q is a user-specified relative error greater than zero. t jump is the time spanned by the cycle jump for the given SDV. s p t 1 + t jump indicates the predicted slope at the time after the cycle jump, which can be obtained by linear extrapolation as:

s p t 1 + t jump s 12 (t 1 ) + s 12 (t 1 ) -s 23 (t 2 ) t cycle t jump . ( 18
)
Hence, Eq. ( 17) ensures that the predicted slope of each SDV at the time after the cycle jump is not too far beyond the known slopes before the attempted jump if a proper relative error q is prescribed.

Substituting Eq. ( 18) into Eq. ( 17) and considering that the time increment is non-negative yield the following expression for the determination of allowable jump length per SDV :

t jump q t cycle |s 12 (t 1 )| |s 12 (t 1 ) -s 23 (t 1 )| . ( 19
)
It is worth noting that Eq. ( 19) only provides the allowed cycle jump length for one SDV at each phase. Nonetheless, the cycle jump lengths will not be necessarily the same for all the SDVs considered since the rate of changes of SDVs in short fiber-reinforced composites with viscoelastic-viscoplastic and damage mechanisms in different orientations can be drastically different. Thus, the common (or global) jump length t common jump is obtained as the minimum of the computed allowed jump lengths for all the SDVs:

t common jump t cycle ⎢ ⎢ ⎢ ⎣ min t jump t cycle ⎥ ⎥ ⎥ ⎦ , ( 2 0 )
where denotes the floor function. Finally, the Heun integrator [START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF] is employed to perform the extrapolation of all the variables necessary for the next Mori-Tanaka TFA micromechanics analysis (control cycle):

S DV (t 1 + t common jump ) S DV (t 1 ) + 1 2 s 12 (t 1 ) + s p (t 1 + t common jump ) t common jump . ( 21
)
Substituting Eq. ( 18) into Eq. ( 21) yields the following formula for the computation of SDVs after the cycle jump:

S DV

Numerical study

Having described the accelerated Mori-Tanaka TFA procedure, it is now of interest to test the accuracy of the cycle jump technique combined with the Mori-Tanaka TFA micromechanics vis-à-vis the cycle-by-cycle reference analysis and understand how various factors might influence the computation accuracy and efficiency.

The material system considered in this work is a thermoplastic short glass fiber-reinforced PA66/GF composite containing 30% fiber mass content, which corresponds to 15.5% fiber volume fraction. The average fiber aspect ratio is 22. Figure 4 shows the skin-shell-core microstructure of the PA66/GF composites, which has been frequently employed in a thin-plate structure of thermoplastic composites obtained by the injection molding manufacturing process. The skin-shell-core microstructure of the injection-molded PA66/GF composites was taken into account in an average sense by considering specific fiber orientation distribution following a density function. The latter is extracted from the micro-CT statistical measurement of an actual skin-shell-core microstructure in the entire through-thickness direction, as shown in Fig. 5. In the mean-field homogenization techniques, the reinforcement is taken into account through its volume fraction, orientation and aspect ratio. The actual position of fibers in the RVE is not considered in Eshelby-based micromechanics as in the Mori-Tanaka scheme.

The material parameters at room temperature and relative humidity of 50%, taken from Chen et al. [START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF], have been directly utilized in generating the results that follow (see Table 1). It should be noted that the damagerelated parameter β is typically positive for metallic materials [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF] for which these types of evolution laws were originally developed. The damage propagation in the case of positive values β occurs in a rather slow manner at the initial loading stage and propagates very fast at the latter loading stage. For thermoplastic polyamide, however, an inverse trend is usually experimentally observed [START_REF] Detrez | Plasticity/damage coupling in semi-crystalline polymers prior to yielding: micromechanisms and damage law identification[END_REF], i.e., fast damage initiation followed by slow growth. This tendency can be well captured by setting negative values for the damage-related parameter β [START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF]. The allowable jump length is determined by a user-defined parameter q which controls the difference of the predicted slope for each SDV the moment after the cycle jump with respect to the slope of the SDV before the attempted jump. The appropriate value for q needs to be determined by comparing the cycle-by-cycle simulation to the cycle jump simulation for a reasonable number of cycles. It should be also noted that the value of q for each SDV can be time-dependent, depending on the evolution mechanisms of each SDV during the cyclic loading. For the sake of simplicity, q is considered a constant value in this work. In the present study, the PA66/GF composites are subjected to piecewise linear stress-controlled cyclic loading with a frequency of f , as shown in Fig. 6. Before passing to the numerical study, let us first define the stress amplitude ratio in cyclic loading as

R σ min σ max , ( 23
)
where σ max and σ min represent the maximum and minimum stresses during the cyclic loading, respectively. Therefore, the cyclic loading history can be described by three parameters: R, σ max and f . First of all, the cycle jumps performed without the control parameter q are compared with cycle-by-cycle reference calculations. The accumulated effective plastic strain p and ductile damage d in the polyamide matrix are obtained for a PA66/GF composite, which is cut from a plate in 45°orientation with respect to the mold flow direction. Stress-controlled cyclic loading in the composite by the normal in-plane stress σ x x 0 has been applied with R 0.1, σ max 50 MPa and frequency of 3 Hz. The global evolution function for each SDV is established based on 30 complete training cycles. The extrapolations are conducted repeatedly using Eq. ( 22) with two fixed jump lengths, namely 50 cycles and 100 cycles, respectively, which are followed by 6 complete control cycle simulations. For all the simulations, 3000 cycles were calculated to reach the target time of 1000 s. Comparing the differences between the cycle jump and cycle-by-cycle simulations in Fig. 7 clearly indicates that the effective plastic strains and ductile damage diverge immediately after the training cycles in the 100-cycle jump length case. The 50-cycle jump length produces stable response, but the solutions are visibly off from the reference solution. This suggests that the cycle jump without error control cannot guarantee that the correct values of the SDVs will be captured and hence may lead to unreliable solutions.

Next, the cycle jump simulations are obtained for four control parameters q 0.2, 0.5, 1.0, 1.5. Cursory examination of Fig. 8 reveals that q 0.2 and q 0.5 predict indistinguishable results relative to the cycleby-cycle reference calculations. Increasing q to 1.0 and 2.0 tends to slightly increase the differences between the cycle jump and the reference solutions, but does not produce a diverged response.

Comparisons of the computation time and the total number of jumped cycles for PA66/GF composites with 0°, 45°and 90°orientations between the cycle jump and cycle-by-cycle simulations are summarized in Table 2. The maximum stresses are 80 MPa, 50 MPa and 40 MPa, respectively, for 0°, 45°and 90°specimens. The stress amplitude ratio R and frequency f for the cyclic loading are 0.1 and 3 Hz for all the cases, respectively. It is observed the cycle jump procedure substantially reduces the number of cycles that need to be conducted relative to the cycle-by-cycle calculations. Increasing q can significantly increase the number of jumped cycles and hence an increase in computational efficiency. Hereafter, to balance the computation efficiency and accuracy, q 0.5 will be employed unless otherwise stated. 8 Comparison of accumulated effective plastic strain and ductile damage in the matrix phase between cycle jump simulation and cycle-by-cycle reference solution for control parameters q 0.2,0.5,1.0,1.5 Table 2 Comparison of the computation time and total number of jumped cycles to reach 3000 cycles between the cycle jump simulation and cycle-by-cycle reference solution for control parameters q 0.2,0.5,1.0,1.5

Loading configuration

Control parameter Number of jumped cycles Computation time 0°No jump 0 23 h 11 min q 0.2 1629 10 h 53 min q 0.5 2140 7 h 3 min q 1.0 2393 5 h 10 min q 1.5 2517 3 h 34 min 45°No jump 0 20 h 48 min q 0.2 1869 7 h 53 min q 0.5 2259 5 h 31 min q 1.0 2453 4 h 1 min q 1.5 2543 3 h 20 min 90°No jump 0 13 h 21 min q 0.2 1714 5 h 54 min q 0.5 2231 3 h 20 min q 1.0 2464 2 h 16 min q 1.5 2567 1 h 51 min

Comparisons of cycle-by-cycle and cycle jump simulations for different loading orientations

Figure 9 presents the comparison of the strain time response and strain envelope between cycle jump simulation and cycle-by-cycle reference solution for three PA66/GF specimens with different orientations, namely θ 0°, 45°and 90°. As observed, the cycle jump solutions show good accordance with the cycle-by-cycle reference solutions for all the loading cases. The employed extrapolation technique enables the automatic determination of a suitable jump length during the cyclic loading. As such, when there are fast variations in the SDVs, such as in the initial portion of loading cycles where significant nonlinearity occurs, the extended Mori-Tanaka TFA approach conducts shorter or no jumps. When the SDVs are stabilized, the extended Mori-Tanaka TFA approach performs longer jumps, thereby eliminating the need for extensive numerical computation during the period of the cycle jump.

Figure 10 shows the stress-strain response for the selected loading cycles for the above three orientations generated using the cycle jump procedure and the cycle-by-cycle reference solution. It is observed nearly identical results between the two solutions, providing good support for the cycle jump technique. Figure 11 presents the comparison of the effective plastic strain and the ductile damage in the matrix phase as a function of time generated by the cycle jump procedure and the cycle-by-cycle analysis. It is concluded that the control cycles between the jumps tend to bring back the solution to the reference results, indicating conducting more intermediate simulations will automatically decrease errors of the cycle jump simulations.

It should be mentioned that, thus far, the extended Mori-Tanaka TFA approach with cycle jump has been verified numerically against the cycle-by-cycle calculation. From the physical problem point of view, the cycle jump micromechanics approach developed in this work is a reliable technique as it gives satisfactory and consistent results compared to the reference solutions obtained without cycle jumps.

Experimental verification

In this section, the full potential of the extended Mori-Tanaka TFA approach with the cycle jump capability is further demonstrated by comparison with the experimental data under high-cycle loading conditions and extensively verified. These experiments were performed for PA66/GF specimens under stress-controlled loading with a frequency of 3 Hz. As mentioned by Arif et al. [START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF], the maximum strain rate is ranging between 0.07 and 0.2 s -1 . The specimens were placed in a climatic chamber with air containing 50% of relative humidity at 65°C. The conditioning process was carried on till all the specimens have a uniform water concentration (∼2.8%). The specimens were sealed prior to testing to maintain a controlled humidity. Figures 12,13 and 14 show the comparison of the maximum and minimum strains as a function of the loading cycle between the cycle jump predictions and experimental response for the PA66/GF composites loaded, respectively, in 0°, 45°and 90°orientations from the mold flow direction. The strain envelopes generated by cycle-by-cycle simulation up to 10 k loading cycles have also been enclosed in each figure for comparison. In the case of 0°specimens, two different stress magnitude ratios and maximum stresses, namely R 0.1, σ max 80 MPa and R 0.3, σ max 90 MPa, are applied. In the case of 45-and 90°specimens, several stress magnitude ratios and maximum stresses, namely R 0.1, σ max 50 MPa; R 0.3, σ max 60 MPa;R 0.5, σ max 70 MPa and R 0.1, σ max 40 MPa; R 0.3, σ max 54 MPa, are performed, respectively. In general, the extended Mori-Tanaka TFA approach appears to reproduce the experimental response for the first 1000 cycles with sufficient accuracy. Particularly worthy of notice is the fact that the experimental response for the 45°specimens under R 0.1, σ max 50 MPa and R 0.3, σ max 60 MPa loading conditions are captured remarkably well. It should be also emphasized that the cycle-by-cycle and cycle jump simulations predict almost identical strain envelopes for every loading configuration, indicating the differences between the experimental and simulated results are not attributed to the cycle jump technique.

Remarks and limitations

Despite their successful applications in automotive industries, the characterization of nonlinear cyclic response of injection-molded short glass fiber polymer composites relies extensively on experiments. In addition to computational efficiency and implementability, the challenges in the modeling of PA66/GF composites include the need to take into account: (1) the skin-shell-core microstructures induced by the injection molding process;

(2) the viscoelastic-viscoplastic response of the polymer matrix; and (3) damage mechanisms in constituent phases. The present development of the extended Mori-Tanaka TFA framework with cycle jump technique attempts to shed light on efficiently simulating the high-cycle response of PA66/GF composites, accounting for the above-mentioned factors that may affect significantly the composite response for the first time. From the physical problem point of view, the cycle jump micromechanics approach developed in this work can be considered a reliable technique as it gives satisfactory and consistent results, but with reduced computational cost compared to the cycle-by-cycle reference solutions. An added advantage is that the proposed micromechanics model, contrary to phenomenological approaches, takes into account explicitly the microstructural characteristics and local constitutive laws.

The successful simulation of the experimental response of actual PA66/GF composites under different loading configurations demonstrates the potential of the proposed technique to predict the response of such composites under high-cycle loading. In general, the extended Mori-Tanaka TFA approach appears to reproduce the experimental response at the initial loading stages with sufficient accuracy for the 0°and 90°specimens and the entire loading stages for the 45°specimens. However, a close examination of Figs. 12, 13 and 14 also reveals that greater differences may occur at high loading cycles. This is particularly true for the specimen loaded with high mean stress, such as 45°specimen with R 0.5, σ max 70 MPa and 90°specimens with R 0.3, σ max 54 MPa. The reason for the observed discrepancies at the high-cycle range and high-meanstress loading condition is several-fold. First of all, the parameters for the viscoelastic polyamide phase used in the extended Mori-Tanaka TFA simulations have been identified from experiments at short characteristic times and thus can only represent the short-or medium-term response. To capture long-term fatigue creep, more viscoelastic branches with higher characteristic times should be included in the polyamide constitutive model. Secondly, the experiment response manifests itself important "banana effect" characterized by larger hysteresis loops with increasing loading cycles, induced by the friction of the cracked phases combined with the viscoelasticity nature of the polyamide, which the extended Mori-Tanaka TFA approach does not capture. It should be emphasized that besides the microstructure variability induced by the injection process, the fatigue experiments exhibit significant discrepancies, especially for composites with random microstructures, at high loading cycles since they were not conducted in a constant environment. The relative humidity and temperature affecting the material response can evolve remarkably during the full cycling test due to several factors, among them energy dissipation and self-heating under the high-cycle range. As far, the extended micromechanics framework does not take into account the effect of temperature and humidity variations, and hence may further add to the differences between the experimental response and numerical simulation at the high-cycle range. It should be mentioned also that at the last stage of fatigue loading, fiber breakage may occur leading to the final macroscopic failure of the composite specimen [START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF][START_REF] Arif | In situ damage mechanisms investigation of PA66/GF30 composite: effect of relative humidity[END_REF]. This mechanism is not taken into account in the present multiscale damage model either. A cycle jump-aided micromechanics methodology has been developed for simulating the high-cycle viscoelastic-viscoplastic response of short glass fiber-reinforced composites with ductile damage. The developed approach employs the previously developed extended Mori-Tanaka method and the transformation field analysis technique to obtain the homogenized behavior of composites with randomly oriented fibers. The response of the polyamide matrix is described by the viscoelastic-viscoplastic-ductile damage constitutive law. The extended theory circumvents the computation challenges imposed by the cycle-by-cycle simulation during the high-cycle loading. It simulates only a minimum number of cycles at the chosen time intervals while most of the cycles between these intervals are skipped using linear extrapolation. The implementation of the cycle jump technique is further facilitated by introducing an adaptive extrapolation scheme, which permits not only the automatic control of cycle jump length, but also good accuracy of the accelerated simulations.

Comparison of the cyclic response of randomly oriented short fiber-reinforced PA66/GF composites with 3 different off-axis orientations from the mold flow direction reveals excellent accordance of the cycle jump and the cycle-by-cycle simulations. The computational effort based on the former approach is drastically reduced relative to the latter. The subsequent successful simulation of the experimental response of the PA66/GF composites with reasonable accuracy further demonstrates the potential of the developed cycle jump-aided micromechanics methodology in predicting the high-cycle (> 100,000 cycles) response of the actual composites with complex microstructures and constitutive behaviors.

Fig. 1

 1 Fig. 1 Modified Mori-Tanaka scheme with multiorientation coated inclusions

Fig. 2 Fig. 3

 23 Fig. 2 Illustration of the cycle jump algorithm

Fig. 4 Fig. 5

 45 Fig. 4 Skin-shell-core microstructure in the through-thickness of PA66/GF composites

Fig. 6

 6 Fig. 6 Piecewise linear stress-controlled cyclic loading

Fig. 7

 7 Fig. 7 Comparison of accumulated effective plastic strain and ductile damage in the matrix phase between cycle-by-cycle solution and cycle jump simulation without cycle control

Fig. 9

 9 Fig. 9 Comparison of the differences in the strain time response (left) and strain envelope (right) between cycle jump simulation and cycle-by-cycle reference solution at three different loading configurations

Fig. 10

 10 Fig. 10 Comparison of stress-strain response between cycle jump simulation and cycle-by-cycle reference solution at the selected loading cycles for three different orientations

Fig. 11

 11 Fig. 11 Comparison of accumulated effective plastic strain and ductile damage in the matrix phase between cycle jump simulation and cycle-by-cycle reference solution for three different orientations

Fig. 13

 13 Fig. 13 Comparison of the maximum and minimum strains at each loading cycle for 45°PA66/GF30 composites under a R 0.1, σ max 50 MPa, b R 0.3, σ max 60 MPa and c R 0.5, σ max 70 MPa

Fig. 14

 14 Fig. 14 Comparison of the maximum and minimum strains at each loading cycle for 90°PA66/GF30 composites under a R 0.1, σ max 40 MPa and b R 0.3,σ max 54 MPa

Table 1

 1 Material parameters[START_REF] Chen | Viscoelastic-viscoplastic homogenization of short glass-fiber reinforced polyamide composites (PA66/GF) with progressive interphase and matrix damage: New developments and experimental validation[END_REF] used in the random Mori-Tanaka TFA simulation (T 23 o C, RH 50%)

		Glass fiber	Polyamide interphase layer	Polyamide matrix
	E	72.4 GPa	2731 MPa	2731 MPa
	υ	0.22	0.3	0.3
	E v1 η v1		8766 MPa 1395 MPa s	8766 MPa 1395 MPa s
	E v2 η v2		13,754 MPa 165,601 MPa s	13,754 MPa 165,601 MPa s
	E v3 η v3		15,010 MPa 457,955 MPa s	15,010 MPa 457,955 MPa s
	E v4 η v4		11,634 MPa 1,307,516 MPa s	11,634 MPa 1,307,516 MPa s
	R 0		4.86 MPa	4.86 MPa
	K		1304.33 MPa	1304.33 MPa
	n H		0.674 47.35 MPa s	0.674 47.35 MPa s m
	m		0.068	0.068
	S			21.607
	β			-1.105
	γ L		1
	γ N		6
	γ ST		2.64
	γ SL		3.7
	γ max κ		0.99 0.62
	σ c		55.20 MPa

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.


domain: spi chim phys sdv

text:
S Fehaili

B Rega

M Courel

P Giampaoli

C Brandam

X Meyer

C Bonazzi

Reaction Engineering for Sponge Cake Baking: Development of a Methodology to extract an Apparent Identifiable Reaction Scheme

The food properties often results from a set of complex chemical reactions initiated by technological treatments. The ultimate objective of this work is to develop rational operating strategies for foods in order to preserve the integrity of the nutritional compounds or to enhance the production of positive newly formed compounds. A "reaction engineering" approach will allow building robust and reliable knowledge on the mechanisms and on the reaction kinetics and then formalizing them into models to optimize the product/process system. The present paper describes the first stage of this approach. It consists in extracting an apparent identifiable reaction scheme on the basis of theoretical knowledge of the system and observable experimental data and on the use of statistical methods.

Introduction

The role of the Maillard reaction in building food quality has been a matter of interest for years. Despite abundant literature, few studies are dealing with the simultaneous effect of physical and chemical variables on the development of Maillard reaction products; in papers dealing with thermal treatment, the latter is merely considered as a mass and energy transfer operation and the link with the chemical variables is not approached [START_REF] Broyart | Modelling heat and mass transfer during the continuous baking of biscuits[END_REF]. One of the main objectives of this project is to understand how the processes act on reactions and can favour the appearance, disappearance or preservation of newly formed compounds during the baking of sponge cake. The reaction system can be described as a complex network of many chemical reactions. But only a small part of all the molecules synthesized by these reactions are observable. A major stake of this work will be to extract a simplified reaction pathway which can give a reliable representation of the complex reaction system based on observable compounds. It will be necessary to determine the reactions which really take place and to specify the occurrence of these reactions during the process. The theoretical knowledge of the system and the use of statistical methods should make it possible to extract an apparent reaction pathway from available experimental information. The methodological approach is as follows: i) Selection of the key-reaction pathways related to previously defined organoleptic and nutritional qualities of the sponge cake and the links occurring between the different reaction pathways [START_REF] Rega | Searching the missed flavour: chemical and sensorial characterisation of flavour compounds released during baking[END_REF]. Three types of reactions will be focussed on: Maillard, lipid oxidation and caramelization reactions. ii) Definition and selection of relevant chemical markers to follow the reactions pathways. iii) Extraction of an apparent but identifiable reaction scheme regarding the observable components from the complete set of theoretical reactions by making as little assumptions as possible. iv) Acquisition of experimental data taking into account the key-variables likely to influence the reaction. v) Validation of the reaction scheme from experimental data. vi) Development of a knowledge based model and selection of the right assumptions to be made: which physical phenomena have to be integrated in the unit operation model and which reactions are to be taken into account and with which kinetic model?

Only the first three points will be developed in this paper.

State of the art on the kinetics of Maillard reaction

The majority of the published work dealing with the kinetics of the Maillard reaction rely on simple kinetic models to describe either loss of reactants (sugars, amino acids, amino acid residues in proteins) or formation of products (Amadori, HMF...). On the use of this approach, van Boekel concluded that the results are very much dependent on experimental conditions and on the extent of the reaction [START_REF] Van Boekel | Kinetic aspects of the Maillard reaction : a critical review[END_REF]. To better understand reaction mechanisms, an advanced approach is necessary, considering reaction pathways in more detail. Kinetic modelling using multiresponse models has proven to be a useful tool to enable further unravelling of reaction mechanisms in the Maillard reaction [4, A first proposed reaction mechanism is tested by the multiresponse analysis where the degradation of the reactants is analysed simultaneously with the formation of the intermediates [START_REF] Martins | Kinetic modelling of Amadori N-(1-deoxy-Dfrucyos-1-yl)-glycine degradation pathways. Part I-Reaction mechanism[END_REF]. It is a very interesting approach which makes it possible to identify the mechanisms of formation/disappearance of the compounds. The multiresponse approach has already been validated for acrylamide formation in model solutions of glucose and asparagines [START_REF] Knol | Toward a kinetic model for acrylamide formation in a glucose-asparagine reaction system[END_REF] and on Maillard reaction in milk [8]. However, this method is not easily applicable in a real food system like sponge cake: first, it is a complex formulated solid where reaction compounds do not react in the same way as in a simple model solution; second, only a small part of all the molecules that can be synthesized by thermal reactions are observable and it would be almost impossible to identify the whole reaction scheme.

Selection of key reaction pathways

The early Maillard reaction consists of the non-enzymatic formation of a covalent bond between the carbonyl group (sugar) and the free amine group (R-NH 2 ). The occurring chemical reactions can be broadly divided into three main stages: (i) the early stage, consisting of the formation of the Amadori product; (ii) the advanced stage, comprising degradation of the Amadori product, and (iii) the final stage, typified by the production of brown polymers and co-polymers, the melanoidins. The degradation of the Amadori product depends on the pH of the system. At pH 7 or below, it mainly undergoes 1,2-enolisation with the formation of 3-deoxyosones. At pH above 7 the degradation of the Amadori compound is thought to involve mainly 2,3 enolisation, where 1-deoxyosones are formed. The key intermediate compounds of this reaction are the αdicarbonyl (3-deoxyosones and 1-deoxyosones) which are also formed by the caramelization reaction [START_REF] Kroh | Caramelisation in food and beverages[END_REF]. Thus, the compounds derived from the degradation of the 1-and 3-deoxyosones are indistinctly generated either by the caramelization or by the Maillard reaction (Fig. 1).

Fig. 1: Interaction between Maillard and caramelization reactions

Compounds derived from the degradation of the deoxyosones

The 3-deoxyosone will undergo successive dehydrations and a displacement of the double bond to form hydroxymethylfurfural (HMF) [START_REF] Kroh | Caramelisation in food and beverages[END_REF], while the 1-deoxyosone generates the 2,3-dihydro-3,5-dihydroxy-6-methyl-4H-pyran-4-one (pyranone) as a typical marker of this reactive pathway [START_REF] Davidek | Simultaneous quantitative analysis of Maillard reaction precursors and products by high-performance anion exchange chromatography[END_REF]. Pure 2,3-dihydro-3,5-dihydroxy-6-methyl-4H-pyran-4-one was found to be converted into maltol, hydroxymaltol, and isomaltol [START_REF] Yaylayan | Stereochemical control of maltol formation in Maillard reaction[END_REF]. Maltol can also be formed directly from the Amadori compound and from disaccharides [START_REF] Yaylayan | Stereochemical control of maltol formation in Maillard reaction[END_REF]. The deoxyosones are also degraded by the Strecker reaction. This reaction puts in presence the deoxyosones on the one hand and the amino acids on the other hand. This is an autocatalytic reaction, which occurs in food with high concentrations in amino acids and strong temperatures or under pressure [START_REF] Cremer | Investigation of the formation of Strecker aldehydes from the reaction of amadori rearrangement products with alpha-amino acids in low moisture model systems[END_REF]. The formed α-aminoketone can induce the genesis of pyrazines.

Compounds derived from the interaction between Maillard reaction and lipid oxidation

Lipid autoxidation reactions generate a multitude of new compounds, in particular, a great number of aliphatic aldehydes. Nevertheless these aldehydes can further react with heterocyclic compounds formed by Maillard reaction. For instance aliphatic aldehydes may react with pyrazines and result in the formation of compounds with long alkyl chains, the alkylpyrazines [START_REF] Shi | The flavour of poultry meat[END_REF].

Schematic overall reaction

A first scheme of the main reaction pathways involved in the three studied reactions and their interactions is proposed in Fig 2 . The concentration of the products as well as various processing variables are likely to influence the advancement of the Maillard, lipid oxidation and caramelization reactions. Manipulating these variables can affect the balance of the various chemical pathways making up the three reactions of interest.

Fig. 2: Proposal of a first scheme of the main pathways involved in the three studied reactions and their interactions, i.e. Maillard, lipid oxidation and caramelization reactions

Selection of relevant markers to follow the reaction pathways

The reaction markers are selected according to their ability to point out different advancement steps of the chemical reactions, taking into account the availability of reliable methods of analysis. A reaction marker can be either a substrate or a product of a reaction pathway, possibly with an organoleptic interest. The Amadori compound and aliphatic aldehydes were chosen as Maillard reaction and lipid oxidation markers respectively. No marker has been selected for the caramelization reaction since the deoxyosone formation intermediates are unstable and therefore difficult to analyse. HMF et 2,3-dihydro-3,5-dihydroxy-6-methyl-4H-pyran-4-one were selected as typical markers of 3-deoxyosone and 1-deoxyosone degradation [START_REF] Davidek | Simultaneous quantitative analysis of Maillard reaction precursors and products by high-performance anion exchange chromatography[END_REF] and maltol, because of its organoleptic interest. To track the Strecker degradation and alkylpyrazine formation, two Strecker aldehydes and pyrazines were chosen.

Extraction of observable reaction pathways

The available reaction markers do not permit observing all the reaction pathways presented in the initial scheme (Fig. 2). The reaction scheme the current state is not identifiable. To reach observable reaction pathways, a stoichiometric matrix is built. The matrix rank indicates the number of independent reactions. An identifiable reaction system is a system which matrix rank equals the number of total reactions [START_REF] Fillon | Data reconciliation based on elemental balances applied to batch experiments[END_REF][START_REF] Moreau | Metabolic engineering of Brevibacterium linens at a macroscopic level[END_REF].

Construction of the stoichiometric matrix

The construction of the stoichiometric matrix requires that both the reaction scheme and the values of the stoichiometric coefficients are known. For example, if we assume that A, B, C, D and E components obey to the following reaction scheme:

A + B → C + D (I) D → B + E (II) A → E (III)
Then the matrix of stoichiometric coefficients is:

The matrix constructed from the initial reaction scheme consists of 21 lines corresponding to the chemical reactions pathways and 13 columns corresponding to selected reaction markers.

Extraction of an observable reaction pathways

The obtained stoichiometric matrix has a rank of 8, which means that only 8 out of the 21 reactions are linearly independent, and therefore observable. Consequently, the building of an observable reaction scheme requires reducing the initial one by: i) eliminating the unlikely reactions based on literature and ii) searching linear dependencies. For instance, in the previous example, reaction III is the sum of reactions I and II; as these reactions are linearly dependent and cannot be distinguished, they must be condensed into reaction III. This implies that some assumptions must be drawn a priori and validated a posteriori.

The observable system obtained (Table 1) is less exhaustive than the overall reaction scheme but it should give a reliable representation of the real and complex reaction system.

Conclusion

The successful implementation of such an approach is highly dependent on the quality of the experimental data that will be obtained. Appropriate experimental tools have been designed to guaranty the reliability of the kinetic models that will be derived from experiments: i) reliable analytical methods giving access to reaction products ii) A discontinuous lab-scale baking oven was specially designed to achieve a good reproducibility of the baking conditions. Hence, the model of heat and mass transfer has to be associated to the kinetic model in order to relate the control variables of the process to the building of the quality through chemical reactions.

  5 th International Conference on Simulation and Modelling in the food and bio-industry, 26-28 June 2008, Dublin, Irland Extended abstract 3/5

  5 th International Conference on Simulation and Modelling in the food and bio-industry, 26-28 June 2008, Dublin, Irland Extended abstract 4/5

Table 1 :

 1 Observable reactions issued from the stoichiometric matrix of rank 8 Strecker aldehyde + α-aminoketone + CO 2 7 2 α-aminoketone → Pyrazines + 2 H 2 O 8 Pyrazines + Aliphatic aldehyde → alkylpyrazine + H 2 O FOODSIM'2008 5 th International Conference on Simulation and Modelling in the food and bio-industry, 26-28 June 2008, Dublin, Irland

	1 Sugar + RNH 2 ↔ Amadori compound + H 2 O
	2 2 Amadori compounds → 4 H 2 O + Maltol + HMF
	3 Sugar → 2 H 2 O + pyranone
	4 Sugar → 3 H 2 O + HMF
	5 Disacharide → 2 H 2 O + Maltol + monosacharide
	6 RNH 2 + deoxyosone →

Acknowledgment

This work is carried out in the framework of the project ANR-06-PNRA-023 with the financial support of the French National Research Agency within the « Programme National de Recherche en Alimentation et nutrition humaine ».