François Boulier

François Lemaire

Alexandre Sedoglavic

Aslı Ürgüplü

Marc Lefranc

Towards an automated reduction method for polynomial ODE models in cellular biology

Keywords: Mathematics Subject Classification (2000). 12H05, 37N25, 62P10. Keywords. computer algebra, differential algebra, cellular biology, system modeling

This paper presents the first version of an algorithmic scheme dedicated to the model reduction problem, in the context of polynomial ODE models derived from generalized chemical reaction systems. This scheme, which relies on computer algebra, is implemented within a new MAPLE package. It is applied over an example. The qualitative analysis of the reduced model is afterwards completely carried out, proving the practical relevance of our methods.

Introduction

This paper is concerned by the model reduction problem in cellular biology. The concerned modeling approach is the well-established one, based on nonlinear differential equations [START_REF] Leloup | Modeling the molecular regulatory mechanism of circadian rhythms in Drosophila[END_REF][START_REF] Conrad | Modeling Molecular Interaction Networks with Nonlinear Differential Equations[END_REF]. However, this paper only considers parametric polynomial ordinary differential equation (ODE) systems, derived from generalized chemical reaction systems by means of the mass-action law [START_REF] Horn | General mass action kinetics[END_REF].

When modeling cellular processes by this approach, one rapidly gets very complicated overparameterized systems. Fitting methods to determine parameter values become difficult to carry out and unreliable. The parameters which reproduce some behaviour of interest are usually far from unique [START_REF] George Von Dassow | The segment polarity network is a robust developmental module[END_REF][START_REF] George Von | Exploring modularity with dynamical models of gene networks[END_REF]. Because of these reasons, the ODE systems need to be reduced for further analysis.

The model reduction problem is common in biological modeling. Though not formulated in these terms, most of the approaches listed in [START_REF] De | Qualitative Approaches to the Analysis of Genetic Regulatory Networks[END_REF] address this issue, in particular the one that considers the qualitative simulations of genetic networks [START_REF] Hidde De | Genetic Network Analyzer: qualitative simulation of genetic regulatory networks[END_REF]. Even in our particular setting, many methods exist for the model reduction, including lumping, sensitivity analysis and multiple time-scale analysis [START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF]. Among all these methods, this paper is concerned by the quasi-steady state approximation (QSSA) [START_REF] Kokotovic | Singular Perturbation Methods in Control: Analysis and Design[END_REF], combined with reparameterization techniques. In our setting, the QSSA relies on the assumption that some of the chemical reactions are much faster than the other ones. Its principle is simple: focusing on the dynamics of the slow reactions, assuming that the fast ones are at quasi-equilibrium. The QSSA has two advantages: it reduces the number of ODE occuring in the system under study and it transforms stiff ODE systems into non stiff ones. Even in the setting of chemical reaction systems, the QSSA has been extensively studied [START_REF] Van | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF].

The authors are developing software which aim is to make the model reduction process as automatic as possible. This paper presents a first version, still incomplete, of such a software. It carries out the following steps:

1. definition of the chemical reaction system to be studied ; 2. approximation by QSSA of the ODE system derived from the input system, leading to a raw reduced model ; 3. reduction of the parameter set of the raw reduced model, leading to a reduced model.

This software is implemented in the new MABSys package of the MAPLE computer algebra system. The QSSA step applies the method presented in [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF], which makes algorithmic the equivalent methods of [START_REF] Van | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF]. The reduction of the parameter set is an exact reduction step, relying on the computation of the system Lie symmetries. It relies on the ExpandedLiePointSymmetry MAPLE package [START_REF] Sedoglavic | Reduction of Algebraic Parametric Systems by Rectification of their Affine Expanded Lie Symmetries[END_REF][START_REF] Sedoglavic | Expanded Lie Point Symmetry (MAPLE package)[END_REF]. All these methods are presented in section 2.

In order to prove the relevance of our methods, an example, borrowed from [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF], is completely carried out in section 3. This example features a single gene regulated by an order n polymer of its own protein (the integer number n is a parameter of this model). The non reduced model involves n + 3 ODE depending on 2 n + 5 parameters. An interesting value is n = 8, which leads to a medium size system. The reduced model involves 3 ODE and n + 6 parameters only. Its qualitative analysis is carried out and one proves that the model exhibits Poincaré-Andronov-Hopf bifurcations if and only if n ≥ 9. Observe that this problem is already solved in [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF] but with only a sketched proof in the case n ≥ 9, a different reduction in the case n ≤ 8, and much less automatic computations.

The methods

Performing QSSA by means of differential elimination

For differential systems arising from generalized chemical reactions systems, there exists a standard way to perform the QSSA, provided that the set of chemical reactions is divided in two parts: the fast ones and the slow ones.

As far as we know, the first clear relationship between this method and the Tikhonov theorem [START_REF] Kokotovic | Singular Perturbation Methods in Control: Analysis and Design[END_REF]Theorem 3.1] was established in [START_REF] Van | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF]. Afterwards, close variants of the same method were rediscovered more or less independently [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF]. Though all these papers present methods, none of them is fully presented in an algorithmic manner. This may at least partly be due to the fact that some steps of the methods require the inversion of a matrix over a residue class ring, a non obvious task which may imply splitting cases. Indeed, it turns out that the whole method is equivalent to a differential elimination process, as shown for the first time in [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. See section 4 for an introduction to differential elimination.

Algorithm 1 DifferentialModelReduction( Ẋ = N V )

Input: The initial parametric ODE system Ẋ = N V derived from a generalized chemical reaction system involving n chemical species and p reactions. The stoichiometry matrix N has dimension n × p and integer entries. The vector X of dependent variables has dimension n. The vector of reaction rates V has dimension n. Its entries are power products of system parameters and dependent variables. Each reaction rate of V (hence each column of N ) is assumed to be tagged "fast" or "slow". Output: a list of dynamical systems in the dependent variables X obtained by quasi-steady state approximation from the initial system or Fail.

1: Split N into two matrices N f (columns of the fast reactions) and N s (columns of the slow reactions). Split V into V f and V s so that the initial ODE model writes Ẋ = N s V s + N f V f .

2: By (say) Gaussian elimination, determine a maximal linearly independent set of columns of N f and remove the other ones, giving N f . Update the vector of reaction rates V f , giving a new vector V f , such that the initial ODE model writes as follows. The entries of V f are linear combinations of elements of V f with rational number coefficients.

Ẋ = N s V s + N f V f .
3: Build a vector F of new dependent variables, having the same dimension as V f . 4: Build the following DAE of the differential polynomial ring K{X ∪ F } where K is the field of the rational fractions in the system parameters:

Ẋ = N s V s + N f F , V f = 0 .
(2.1)

5: R := a ranking F ≫ X eliminating F w.r.t. X. 6: [C 1 , . . . , C t ] := Rosenfeld-Gröbner((2.1), R) 7:
If there exists some dependent variable X i and some regular differential chain

C k such that NF(X i , C k ) is not a rational fraction in the variables X then Fail. 8: Return the list Ẋ = NF( Ẋ, C k ) for 1 ≤ k ≤ t.
The DifferentialModelReduction algorithm is summarized in Algorithm 1. One illustrates it over a famous example: the beginning of the Henri-Michaelis-Menten reduction of the basic enzymatic reaction system:

E + S k1 ----→ ← ---- k2 C k3 ----→ E + P .
(2.

2)

The initial system of ODE writes:

Ẋ = N V i.e.     Ė Ċ Ṡ Ṗ     =     -1 1 1 1 -1 -1 -1 1 0 0 0 1     •   k 1 E S k 2 C k 3 C   (2.3)
where X is the vector of the chemical species, N is the system stoichiometry matrix and V is the vector of the reaction rates. The stoichiometry matrix is built as follows: it involves one row per species and one column per reaction. The entry at row r, column c is the number of molecules of species r produced by the reaction c (i.e. the number of times species r occurs on the reaction right-hand side minus the number of times it occurs on the reaction left-hand side). The rate of a reaction is the product of the left-hand side species (with multiplicities) times the reaction rate constant (the parameter over the arrow). Split the stoichiometry matrix N into two matrices N f and N s putting the columns which correspond to fast reactions in N f and the ones which correspond to slow reactions in N s . Split accordingly the rows of the vector V into two vectors V f and V s . One gets a formula Ẋ = N s V s + N f V f . Over system (2.3), one gets:

    Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 1 -1 -1 1 0 0     • k 1 E S k 2 C . (2.4)
Determine a maximal linearly independent set of columns of N f (i.e. a basis of that matrix) and remove the other ones, giving a new matrix N f . Update the vector of reaction rates V f , giving a new vector V f such that N f V f = N f V f . Over the example, removing the second column, one gets a new formula Ẋ =

N s V s + N f V f which is equivalent to formula (2.3):     Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 -1 0     • k 1 E S -k 2 C . (2.5)
Replace the vector V f by a vector F of new dependent variables F i . The slow variety 1 is defined by letting the entries of V f all equal to zero. The DAE to be considered for quasi-steady state approximation is

Ẋ = N s V s + N f F , V f = 0 . (2.6)
Over the example, one gets:

    Ė Ċ Ṡ Ṗ     =     1 -1 0 1     • k 3 C +     -1 1 -1 0     • F 1 , k 1 E S -k 2 C = 0 .
(2.7) Expand this system in order to get a set of differential polynomials in the differential polynomial ring K{X ∪ F } where K = Q(k 1 , k 2 , k 3 ) and, using a ranking F ≫ X eliminating the new unknowns F , simplify it by means of a differential elimination process. One gets the following system:

F 1 = k 3 k 1 E S (k 1 S + k 2 ) k 2 (k 1 S + k 1 E + k 2 ) , Ė = k 2 1 E 2 k 3 S k 2 (k 1 S + k 1 E + k 2 ) , Ṗ = k 3 k 1 E S k 2 , Ṡ = - k 3 k 1 E S (k 1 S + k 2 ) k 2 (k 1 S + k 1 E + k 2 ) , C = k 1 E S k 2
Readers used to the Michaelis-Menten formula do probably not recognize it in the above result. Indeed, some further simplifications need to be done but these simplifications are actually not related to the QSSA. They involve considerations on initial conditions and parameter renaming. See [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF] for a complete study.

Reduction of the parameter set and reparameterization

Two packages are available for reducing the parameter set and reparameterizing differential systems. The first one is the ExpandedLiePointSymmetry package which is not specifically designed for the manipulation of biochemical systems. The second one is the MABSys package which is mostly oriented for the manipulation of models coming from biochemical reactions. They perform automatically two kinds of changes of coordinates:

1. changes of coordinates which reduce the number of parameters (by computing Lie symmetries of the differential system) ; 2. changes of coordinates which make some parameters appear as factors in the right-hand sides of differential equations (by computing Lie symmetries of the nondifferential system which defines the system steady points).

To illustrate these features, consider the following ODE, which is borrowed from system (3.4) and slightly simplified:

Ġ = θ (1 -G) -α G.
Greek letters denote parameters. The solutions of the steady point equation 0 = θ (1 -G) -α G are not changed when α and θ are both multiplied by any nonzero constant λ. This suggests to replace α by α θ, which leads to the new ODE:

Ġ = θ (1 -G -α G).
The parameter θ appears as a factor of the right-hand side of the ODE. The steady point does not depend on θ anymore. Now, the ODE is left unchanged when the time t is multiplied by λ and θ is divided by λ. This suggests to replace t by t/θ. This leads to the new ODE:

Ġ = 1 -G -α G.
An important feature of the MABSys package is that it restricts itself to symmetries of type scaling because they preserve the positivity of the system variables and parameters: a crucial property for studying the system qualitative behaviour.

An important issue is not yet satisfactorily solved: the parameters that can be removed or rewritten by means of symmetries is not uniquely defined. Designing a package interface which precisely defines the package output is thus far from easy. Different choices may lead to different systems which lead to qualitative analyses of various difficulties.

The example

One considers the genetic circuit depicted in Figure 1. The single gene is regulated by an order n polymer of its own protein. The integer number n is a parameter of the system. This study was motivated by the activity of a working group aiming at modeling the circadian clock of the green alga ostreococcus tauri. See [START_REF] Mcclung | Plant Circadian Rhythms[END_REF] for a survey on circadian rhythms and [START_REF] Fall | [END_REF]Chapter 9] or [START_REF] Goldbeter | Biochemical Oscillations and Cellular Rhythms: The Molecular Bases of Periodic and Chaotic Behaviour[END_REF][START_REF] Françoise | Oscillations en biologie[END_REF] for more general texts about oscillations in biology.

The addressed question is: does there exist biologically meaningful (i.e. positive) parameter values which make this circuit oscillate ? A related but easier problem consists of searching for the existence of parameter and variable positive values which give rise to Poincaré-Andronov-Hopf bifurcations. See [START_REF] Hale | Dynamics and Bifurcations[END_REF]Chapter 11], or [START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF]Section I.16]. In the neighborhood of a Poincaré-Andronov-Hopf bifurcation indeed, a stable steady point of the model under study gives birth to a small stable limit cycle under some general hypotheses. Note that searching for Poincaré-Andronov-Hopf bifurcations is not as general as searching for limit cycles: first, some Poincaré-Andronov-Hopf bifurcations (the subcritical ones) do not imply the existence of stable limit cycles; second, there may exist limit cycles not related to Poincaré-Andronov-Hopf bifurcations.

Our approach is decomposed in three parts. First one reduces the initial model to a three-variable model by means of quasi-steady state approximation. Second, one reduces the parameter set of the three-variable model. Third one proves by computer algebra methods that Poincaré-Andronov-Hopf bifurcation occur for positive values of the parameters if and only if n ≥ 9. As pointed out previously, this is not sufficient to prove the presence of limit cycles in the case n ≥ 9. However, this result is confirmed by extensive numerical simulations. Other approaches could have been applied. There exist software packages such as AUTO or XPPAUT [START_REF] Doedel | AUTO software for continuation and bifurcation problems in ODEs[END_REF][START_REF] Ermentrout | Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students[END_REF] which locate Poincaré-Andronov-Hopf bifurcations by means of numerical calculations. They allow one to evidence the existence of Poincaré-Andronov-Hopf bifurcations but not to prove their absence, and thus cannot be used to discard a model. Moreover, they can only tackle particular values of n. Theoretically, the existence or the absence of Poincaré-Andronov-Hopf bifurcations can be decided algebraically [START_REF] Kahoui | Deciding Hopf bifurcations by quantifier elimination in a software-component architecture[END_REF][START_REF] Wang | Stability Analysis of Biological Systems with Real Solution Classification[END_REF][START_REF] Gatermann | Computational algebra for bifurcation theory[END_REF][START_REF] Gatermann | Toric Ideals and graph theory to analyze Hopf bifurcations in mass action systems[END_REF][START_REF] Niu | Algebraic Approaches to Stability Analysis of Biological Systems[END_REF]. In particular, it can be decided by means of computer algebra methods which rely on Sturm sequences computations and algebraic elimination. Practitioners usually seem to avoid these methods because of their huge complexity in the worst case. See however [START_REF] Anai | Algebraic Methods for Solving Real Polynomial Constraints and their Applications in Biology[END_REF][START_REF] Yoshida | Inquiry into Conditions for Cell-Type Diversity of Multicellular Organisms by Quantifier Elimination[END_REF] for applications in biology. In particular, the QEPCAD [START_REF] Brown | QEPCAD B: a program for computing with semi-algebraic sets using CADs[END_REF] package, which is based on quantifier elimination methods, could not solve the problem addressed in this section. The REDLOG package [START_REF] Dolzmann | Redlog: computer algebra meets computer logic[END_REF] and the software described in [START_REF] Kahoui | Deciding Hopf bifurcations by quantifier elimination in a software-component architecture[END_REF] rely on QEPCAD for the quantifier elimination process. An attempt to solve the addressed problem using the RAGLib library [START_REF] Safey | RAGLib (Real Algebraic Library Maple package[END_REF] did not succeed. According to specialists, these methods are optimized to tackle the worst case while our problems exhibit more generic features.

The abstract model depicted in Figure 1 is closely related to models studied by Goodwin and Griffith in the 60's [START_REF] Goodwin | Temporal Organization in Cells[END_REF][START_REF] Goodwin | Advances in Enzyme Regulation[END_REF][START_REF] Griffith | Mathematics of Cellular Control Processes. I. Negative Feedback to One Gene[END_REF]. It features a negative feedback loop, one of the core ingredients for generating oscillations [START_REF] Fall | [END_REF]. Griffith considered a model of a gene regulated by a polymer formed of n copies of its own protein.

The same problem is studied here, but in a slightly more general case, where gene activation is not assumed to be fast. Although we do focus here on biology, it should be stressed that a cooperativity of order 9 is not as unrealistic as it may seem. In particular, gene regulation by an octamer has been reported [START_REF] Selleck | A histone fold TAF octamer within the yeast TFIID transcriptional coactivator[END_REF]. Moreover, an effective cooperativity of order 9 may also be obtained as a consequence of reducing a higher-dimensional, more realistic, model to that of Figure 1. Finally, our conclusions are consistent with those of Griffith [START_REF] Fall | [END_REF] and of other works devoted to more sophisticated variants of the Goodwin model [START_REF] Ruoff | The Temperature-Compensated Goodwin Model Simulates Many Circadian Clock Properties[END_REF][START_REF] Ruoff | The Goodwin Model: Simulating the Effet of Cycloheximide and Heat Shock on the Sporulation Rhythm of Neurospora crassa[END_REF][START_REF] Kurosawa | Comparative Study of Circadian Clock Models, in Search of Processes Promoting Oscillation[END_REF].

In the next sections, all the computations are sketched. A detailed trace of these computations is available at [20].

The initial model

The model of Figure 1 is translated as a system of generalized chemical reactions (observe that transcription and translation are not balanced reactions). The variables G and H represent the state of the gene. The mRNA concentration and the concentration of the protein translated from the mRNA are represented respectively by M and P . The n types of polymers of P are denoted by P = P 1 , P 2 , . . . , P n . Greek letters and k - i , k + i (1 ≤ i ≤ n-1) represent parameters:

G + P n α ----→ ← ---- θ H , G ρ f ----→ G + M , H ρ b ----→ H + M , M β ----→ M + P , M δM ----→ ∅ , P δP ----→ ∅ , P i + P k + i ----→ ← ---- k - i P i+1 (1 ≤ i ≤ n -1). (3.1)
This generalized chemical reaction system can now be canonically translated as a system of parametric ordinary differential equations, denoting

A i = (k - i P i+1 - k + i P i P ).
Variables G, H, M, P = P 1 , . . . , P n are dependent variables. They all represent species concentrations except G and H, which should rather be viewed as "random variables".

Ġ = θ H -α G P n , Ḣ = -θ H + α G P n , Ṁ = ρ f G + ρ b H -δ M M , Ṗ = β M -δ P P + 2 A 1 + A 2 + • • • + A n-1 , Ṗi = -A i-1 + A i (2 ≤ i ≤ n -1) , Ṗn = -A n-1 + θ H -α G P n . (3.2)
This system involves n + 3 differential equations depending on 2 n + 5 parameters.

The raw reduced model

In order to apply a quasi-steady state approximation, it is assumed that the n -1 chemical reactions describing the polymerization of the protein are fast compared to the other ones. Then, according to the technique sketched in section 2.1, one gets an approximation of system (3.2) by replacing each expression A i by a new dependent variable F i (1 ≤ i ≤ n -1) and by augmenting this system by the n -1 following algebraic equations:

0 = k + i P P i -k - i P i+1 , (1 ≤ i ≤ n -1). (3.3)
It is now sufficient to eliminate the F i from the so obtained differential-algebraic system.

The MABSys package allows us to define the initial model and to compute the raw reduced model. Unfortunately, computations cannot be performed by keeping a symbolic value for n (this is a classical restriction of symbolic computation methods). However, the MABSys package can be applied for many different values of n and the general formula can be inferred.

Computations are as follows. The chemical reaction system (3.1) is entered by using the NewReaction function. Then the ModelReduce function (which is an enhanced implementation of DifferentialModelReduction, based on the RegularChains package) is applied. Then, H is replaced by γ 0 -G, introducing a new positive parameter γ 0 . Last n new K i parameters are introduced for legibility with the convention K 0 = 1.

K i = k + 1 • • • k + i k - 1 • • • k - i •
The obtained raw reduced model writes:

Ġ = θ (γ 0 -G) -α K n-1 P n G, Ṁ = ρ b (γ 0 -G) + ρ f G -δ M M, Ṗ = n θ (γ 0 -G) -n α K n-1 P n G -δ P P + β M n-1 i=0 (i + 1) 2 K i P i • (3.4)
Observe that, in principle, there is no need to introduce the K i parameters since this simplification should be performed in the next section. However, the MABSys package does not find this particular change of coordinates but a slightly more complicated one. We chose to perform this natural simplification interactively for the sake of the legibility of our paper.

The reduced model

The raw reduced model (3.4) can now be simplified by rescaling all parameters and variables. Observe that computations are not completely automatic: the practitionner needs to choose the parameters to keep and the ones to eliminate.

Computations are made in two main steps as follows. The first step uses the function RemoveParameterByScalings which is an interface to some functionnalities of the ExpandedLiePointSymmetry package. It removes the two parameters β and δ P . The second step uses MABSys. The CylindrifySteadyPoints function is called in order to make some parameters appear as factors in the right-hand sides of the differential equations. Thus, these parameters do not appear anymore in the algebraic steady points system. For n ≥ 2, two further changes of coordinates are applied interactively for the sake of the legibility of the result (the case n = 1 is slightly different, see [20]). In all cases, the reduced model writes:

Ġ = θ (γ 0 -G -G P n ), Ṁ = ((ρ f -1) G + γ 0 -M ) δ M , Ṗ = M -P + n ρ b (γ 0 -G -G P n ) n 2 α n-1 P n-1 + n-2 i=0 (i + 1) 2 K i α i P i • (3.5)
The changes of coordinates, summarized below, are obtained by composing the ones automatically computed by the MABSys package with the extra ones, performed interactively. The old variables are expressed as functions of the new (overlined) ones. The bars are removed afterwards for legibility. For n ≥ 2, the change of variables are as follows (they are slightly different in the case n = 1, see [20]):

M = M α δ P K n-1 β , t = t δ P , P = P α K n-1 , G = G ρ b α K n-1 θ , α = K n-1 θ δ P α n , θ = θ δ P , ρ f = ρ f θ δ M δ 2 P ρ b β , ρ b = θ δ M δ 2 P ρ b β , δ M = δ M δ P , γ 0 = γ 0 ρ b α K n-1 θ , K n-1 = K n-1 n-1 , K i = K i K i n-1 (1 ≤ i < n) . (3.6)

Qualitative analysis of the reduced model

In this section and the next ones, one proves the following proposition. Though the MABSys package provides functions related to the qualitative analysis of differential systems (in particular, functions for studying the presence of Poincaré-Andronov-Hopf bifurcations), most of the study performed in this section needed to be performed interactively. Observe that one does not need to distinguish the roles of the variables from the ones of the parameters at this step. The zeroes of this system provide the steady points of system (3.5). In general, one cannot compute a Gröbner basis if a symbolic n is left as an exponent, but in our case, a generic Gröbner basis exists. The ordering was chosen carefully in order to achieve two important properties:

1. the leading monomials are plain variables ; 2. the right-hand sides of the Gröbner basis equations are positive. The first property implies that the quotient ring defined by the Gröbner basis is a free algebra: a polynomial ring. The second property implies that there are no constraints on the values that can be assigned to the variables and parameters occuring in the right-hand sides of the Gröbner basis equations.

Computations were performed using the Basis function of the MAPLE Groebner package. The leading monomials appear on the left-hand sides of the equations:

γ 0 = P (1 + P n ) P n + ρ f , M = P, G = P P n + ρ f •

The Jacobian matrix

In order to study Poincaré-Andrononv-Hopf bifurcations of system (3.5), one needs to consider the Jacobian matrix of that system, evaluated over the system steady points. Thanks to the striking properties of the Gröbner basis computed in the above paragraph, one just needs to replace each element of the generic Jacobian matrix by its normal form w.r.t. the Gröbner basis and to forget the steady point equations. The normal form of the Jacobian matrix writes as follows:

J =         -θ (1 + P n ) 0 - n θ P n P n + ρ f (-1 + ρ f ) δ M -δ M 0 - n ρ b (1 + P n ) B 1 B - n 2 ρ b P n + P n + ρ f B (P n + ρ f )         where B = n 2 α n-1 P n-1 + n-2 i=0 (i + 1) 2 K i α i P i .
The parameters α and K i only occur in B. It is thus possible to assign arbitrary positive values to B without perturbating the values of the other expressions involved in the matrix elements. One can thus consider B as a new parameter.

If n ≤ 8 then no Poincaré-Andronov-Hopf bifurcation arises

In order to prove that no Poincaré-Andronov-Hopf bifurcation arises for positive values of variables and parameters whenever n ≤ 8, it is sufficient to prove that the three Hurwitz determinants c 0,0 , c 1,0 and c 2,0 are positive, thanks to Proposition 4.4 and Definition 4.6. These determinants are defined as follows:

c 0,0 = 1, c 1,0 = a 1 , c 2,0 = a 1 a 2 -a 3 .
where x 3 + a 1 x 2 + a 2 x + a 3 denotes the characteristic polynomial of the matrix J. The determinant c 0,0 is positive. The determinant c 1,0 is positive also since its numerator and denominator are linear combinations of power products of positive variables and parameters, with positive coefficients:

c 1,0 = θ B P 2 n + (1 + n 2 ρ b + δ M B + θ B + θ ρ f B) P n + (1 + δ M B + θ B) ρ f B (P n + ρ f ) •
In the sequel, one proves that c 2,0 is positive if n ≤ 8. The denominator B (P n +ρ f ) 2 of this rational fraction is positive. The numerator is a sum of 59 monomials, only two of which have negative coefficients:

-n δ M θ ρ f B P n (P n + ρ f ).

One thus studies the positivity of the numerator, which is a polynomial in P n . One performs a change of coordinates, renaming P n as P . This polynomial has the form

d 0 ρ 2 b + d 1 ρ b + d 2
where d 0 and d 1 are linear combinations of power products of positive variables and parameters, with positive coefficients. Thus c 2,0 is positive for each ρ b > 0 if and only if d 2 is nonnegative. One thus studies the nonnegativity of d 2 , factoring out its positive coefficient P + ρ f . This polynomial has the form:

d 2 P + ρ f = e 0 ρ f + e 1
where e 1 is a linear combination of power products of positive variables and parameters, with positive coefficients. Thus d 2 is nonnegative for each ρ f > 0 if and only if e 0 is nonnegative. One thus studies the nonnegativity of e 0 , which has the form:

e 0 = f 0 P 2 + f 1 P + f 2
where f 0 and f 2 are linear combinations of power products of positive variables and parameters, with positive coefficients. Thus e 0 is nonnegative if and only if e 0 (viewed as a univariate polynomial in P ) has no nonnegative root. Since f 2 /f 0 is the product of the roots, the two roots have the same sign and are nonzero. For both roots to be positive, it is necessary and sufficient to have f 1 negative (since f 1 is the opposite of the sum of the roots) and a positive discriminant of e 0 w.r.t. P (in order to have real roots). These polynomials write:

f 1 θ = 2 B (B δ M + 1) θ + (δ M B) 2 + (2 -n) δ M B + 1
and disc (e 0 , P )

θ 2 = -4 n B (δ M B) (δ M B + 1) θ +(δ M B) 4 -2 n (δ M B) 3 + (n 2 -4 n -2)(δ M B) -2 n (δ M B) + 1.
Towards an automated reduction method for polynomial ODE models 13

Denoting δ M B by δ, one has e 0 nonnegative if and only if the conditions 0 < θ < θ 0 and 0 < θ < θ 1 are not satisfied simultaneously, i.e. if and only if θ 0 and θ 1 are not positive simultaneously, where:

θ 0 = δ 4 -2 n δ 3 + (n 2 -4 n -2) δ 2 -2 n δ + 1, θ 1 = -δ 2 + (n -2) δ -1.
These two polynomials are reciprocal polynomials. Thus, denoting λ = δ + 1/δ, one has:

θ 0 /δ 2 = λ 2 -2 n λ + n 2 -4 n -4, θ 1 /δ = -λ + n -2.
For θ 1 to be positive, it is necessary that n > λ + 2. For θ 0 to be positive, it is necessary that

n < λ + 2 -2 √ λ + 2 or n > λ + 2 + 2 √ λ + 2.
Thus, for both θ 0 and θ 1 to be positive, it is necessary that n > λ + 2 + 2 √ λ + 2. Since λ = δ + 1/δ with δ > 0, one has λ ≥ 2 whence n > 8. Thus c 2,0 is positive if and only if n ≤ 8. This concludes the proof of the left to right implication of Proposition 3.1.

K i = 1 for each 1 ≤ i < n.
The initial values are: P (0) = 1.17, M (0) = 0.8 and G(0) = 0.002.

Let n be an integer number greater than or equal to 9. One exhibits a Poincaré-Andronov-Hopf bifurcation by applying Proposition 4.5. Take δ = 1. Then, since δ = δ M B one can take δ M = 1/2 and B = 2 (observe that B denotes an expression that is greater than 1). The conditions θ < θ 0 and θ < θ 1 then permit to take θ = 1/20. Then the polynomial e 0 has two positive roots (in the P variable) and one can take the value P = 10, which is enclosed between the two roots of e 0 , in order to ensure e 0 < 0. Now, the curve d 2 = 0 is a decreasing function of n, bounded by (say), 600. Taking ρ f = 600 thus ensures that d 2 < 0. The positive root of d 0 ρ2 b + d 1 ρ b + d 2 = 0 provides a value of ρ b which cancels c 2,0 . Its analytic formula is:

ρ b = -3843 + √ 1461560 n + 1640961 20 n 2 •
These values ensure that c 0,0 > 0, c 1,0 > 0, c 2,0 = 0 and c 2,1 = -a 3 < 0 i.e. that a Poincaré-Andronov-Hopf bifurcation occurs. This concludes the proof of Proposition 3.1. Figure 2 shows an oscillation in the neighborhood of a Poincaré-Andronov-Hopf bifurcation.

Appendix

On differential elimination

The differential elimination theory is a subtheory of the differential algebra [START_REF] Fels | Differential Algebra[END_REF][START_REF] Robert | Differential Algebra and Algebraic Groups[END_REF].

See also [START_REF] Wang | Elimination Practice: Software Tools and Applications[END_REF]. The differential elimination processes that are presented in this paper take as input two parameters: a system of polynomial (thus nonlinear) differential equations, ordinary or with partial derivatives 2 and a ranking. They produce on the output an equivalent finite set of polynomial differential systems, which are simpler, in the sense that they involve some differential equations which are consequences of the input system but were somehow hidden. The output may consist of more than one differential system because the differential elimination process may need to split cases. The set of the differential equations which are consequences of the input system forms a so-called differential ideal of some polynomial differential ring. Since this ideal is an infinite set, a natural question arises: how does the process select the finitely many differential equations which appear in the output system ? This is indeed the role of the rankings.

A differential ring (resp. field) is a ring (resp. field) R endowed with a derivation (this paper is restricted to the case of a single derivation but the theory is more general) i.e. a unitary mapping R → R such that (denoting ȧ the derivative of a):

˙ (a + b) = ȧ + ḃ , ˙ (a b) = ȧ b + a ḃ . (4.1)
Observe that, theoretically, the derivation is an abstract operation. For legibility, one views it as the derivation w.r.t. the time t. Algorithmically, one is led to manipulate finite subsets of some differential polynomial ring R = K{U } where K is the differential field of coefficients (in practice, K = Q, Q(t) or Q(k 1 , . . . , k r ) where the k i denote parameters that would be assumed to be algebraically independent)

and U is a finite set of dependent variables3 . The elements of R, the differential polynomials are just polynomials in the usual sense, built over the infinite set, denoted ΘU , of all the derivatives of the dependent variables.

Definition 4.1. A differential ideal of a differential ring R is an ideal of R, stable under the action of the derivation.

Let F be a finite subset of a differential ring R. The set of all the finite linear combinations of various orders derivatives of elements of F , with elements of R for coefficients, is a differential ideal. It is called the differential ideal generated by F . An ideal A is said to be radical if a ∈ A whenever there exists some nonnegative integer p such that a p ∈ A. The radical of an ideal A is the set of all the ring elements a power of which belongs to A. The radical of a (differential) ideal is a radical (differential) ideal.

Theorem 4.2. Let R be a differential polynomial ring and F be a finite subset of R. A differential polynomial p of R lies in the radical of the differential ideal generated by F if and only if it vanishes over every analytic solution of F .

Proof. [START_REF] Fels | Differential Algebra[END_REF]chap. II,[START_REF] Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Cox | Ideals, Varieties and Algorithms. An introduction to computational algebraic geometry and commutative algebra[END_REF] or [START_REF] Boulier | A computer scientist point of view on Hilbert's differential theorem of zeros[END_REF].

The Rosenfeld-Gröbner algorithm [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF] solves the membership problem to radical differential ideals. To present it, one needs to define the concept of ranking. Definition 4.3. If U is a finite set of dependent variables, a ranking over U is a total ordering over the set ΘU of all the derivatives of the elements of U which satisfies: a < ȧ and a < b ⇒ ȧ < ḃ for all a, b ∈ ΘU .

Let U be a finite set of dependent variables. A ranking such that, for every u, v ∈ U , the ith derivative of u is greater than the jth derivative of v whenever i > j is said to be orderly. If U and V are two finite sets of dependent variables, one denotes U ≫ V every ranking such that any derivative of any element of U is greater than any derivative of any element of V . Such rankings are said to eliminate U w.r.t. V .

Assume that some ranking is fixed. Then one may associate with any differential polynomial f ∈ K{U } \ K the greatest (w.r.t. the given ranking) derivative v ∈ ΘU such that deg(f, v) > 0. This derivative is called the leading derivative or the leader of f . Rankings permit to define leaders. Leaders permit to use differential polynomial as rewrite (substitution) rules. Assume that f = a d v d + • • • + a 1 v + a 0 is a differential polynomial with leader v (the coefficients a i are themselves differential polynomials). Then the equation f = 0 can be written:

v d -→ - a d-1 v d-1 + • • • + a 1 v + a 0 a d • (4.2)
It can be used afterwards as a rule to simplify any differential polynomial g such that deg(g, v) ≥ d or deg(g, v (k) ) > 0 where v (k) denotes any proper derivative of v. There are precise algorithms for performing these sorts of substitution by finite sets of rewrite rules: Ritt's reduction algorithm or the normal form algorithm [4, algorithm NF].

The Rosenfeld-Gröbner algorithm gathers as input a finite system F of differential polynomials and a ranking. It returns a finite family (possibly empty) C 1 , . . . , C r of finite subsets of K{U } \ K, called regular differential chains. Each system C i defines a differential ideal C i (it is a characteristic set of C i ) in the sense that, for any f ∈ K{U }, we have

f ∈ C i iff NF(f, C i ) = 0 . (4.
3)

The relationship with the radical A of the differential ideal generated by F is the following:

A = C 1 ∩ • • • ∩ C r . ( 4
.4) When r = 0 we have A = K{U }. Combining both relations, one gets an algorithm to decide membership in A. Indeed, given any f ∈ K{U } we have:

f ∈ A iff NF(f, C i ) = 0, 1 ≤ i ≤ r . (4.5)
The differential ideals C i do not need to be prime. They are however necessarily radical. The NF( • , C i ) function permits to compute canonical representatives of the residue classes of the differential ring R/C i .

4.2. On Poincaré-Andronov-Hopf bifurcations 4.2.1. Non parametric systems. Let ẋ = F (x) be a differential system in m dependent variables. The steady points of the differential system are the zeros of the system (that we assume to be polynomial or rational) F (x) = 0. To each steady point, one may associate a linear system ẋ = J x where J is the m × m jacobian matrix of the differential system, evaluated over the steady point. The stability of the steady state is determined by the eigenvalues of J. It is stable if and only if all eigenvalues have negative real parts. Thus to each steady point, one may associate the characteristic polynomial

C(σ) = σ m + a 1 σ m-1 + • • • + a m (a 0 = 1) of J.
Thanks to the Routh-Hurwitz criterion, the stability of the steady points can be studied by analyzing the sign of the Hurwitz determinants c k,0 . These ones can be directly computed from the coefficients of the characteristic polynomial, as shown below. Following [28, Section I.13], compute the Sturm sequence:

p 0 (ω) = ℜ C(i ω) i m , p 1 (ω) = -ℑ C(i ω) i m (4.6) p k+2 (ω) = -rem(p k , p k+1 , ω) (k ≥ 0). Denote p k (ω) = c k,0 ω m-k + c k,1 ω m-k-2 + c k,2 ω m-k-4 + • • •
Observe that the computation of p k must be performed carefully (e.g. using subresultant sequences) to ensure that c k,0 actually is a Hurwitz determinant. See [START_REF] Kahoui | Deciding Hopf bifurcations by quantifier elimination in a software-component architecture[END_REF]. Indeed,

c 0,0 = 1, c 1,0 = a 1 , c 2,0 = a 1 a 2 -a 3 , . . . , c m,0 = a m c m-1,0 .
Towards an automated reduction method for polynomial ODE models 17

The two following propositions are well known. The first one is nearly a corollary to the Routh Theorem [START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF]Theorem 13.4].

Proposition 4.4. With the same notations, if all the Hurwitz determinants c k,0 are positive, apart perhaps c m,0 , then J has no pure imaginary eigenvalue.

Proof. If all the Hurwitz determinants c k,0 are positive (0 ≤ k < m) then they are a fortiori nonzero. Assume J has pure imaginary eigenvalues ±i ω (they are necessarily conjugate). These values ±ω are then common zeros of p 0 and p 1 . The gcd of p 0 and p 1 has thus degree greater than or equal to 2. This gcd is the last nonzero polynomial in the sequence p 0 , . . . , p m-1 . Thus one polynomial p k with 0 ≤ k < m must vanish identically. Therefore the corresponding Hurwitz determinant c k,0 must vanish also.

Parametric systems.

The differential systems encountered in biological modelling involve parameters. Let ẋ = F (x, θ) be a differential system in m variables and p parameters θ. If some real values are assigned to the parameters then one gets a system such as the one described in the previous section. If these real values continuously vary then the steady points and their associated eigenvalues continuously vary also. Definition 4.6. With notations as above, a Poincaré-Andronov-Hopf bifurcation arises for a steady point when all the eigenvalues associated to the steady point have negative real parts except one complex conjugate pair, which crosses the imaginary axis because of a variation in the system parameters.

In computer algebra, an important point is to avoid to compute the steady points, i.e. not to solve the system F (x, θ) = 0. The Hurwitz determinants can be computed generically. They depend on the system parameters. Their sign is studied modulo the ideal I generated by the polynomial system F (x, θ) = 0. The absence of Poincaré-Andronov-Hopf bifurcation is established, thanks to Proposition 4.4 and Definition 4.6, by proving that the Hurwitz determinants c 0,0 , . . . , c m-1,0 are positive for all x and θ, considering that x and θ satisfy F (x, θ) = 0 plus, usually, some extra (positivity) conditions such as x, θ > 0.

The Hurwitz determinants c k,0 get reformulated by computing their normal forms ck,0 w.r.t. any Gröbner basis of the ideal I. Reference books for the Gröbner basis theory are [START_REF] Cox | Ideals, Varieties and Algorithms. An introduction to computational algebraic geometry and commutative algebra[END_REF][START_REF] Becker | Gröbner Bases: a computational approach to commutative algebra[END_REF]. Indeed, the difference c k,0 -ck,0 belongs to I. Over any steady point of the differential system, it is thus zero, thus the two polynomials c k,0 and ck,0 have the same value hence the same sign.

In practice moreover, Gröbner bases can be computed in dimension zero. Computing in dimension zero corresponds to some generic computation, which may be false for particular values of the system variables and parameters. However, in biological models, parameters (and thus variables) have no accurate values and zero dimensional computing makes sense.

Conclusion

We have presented the first version of an algorithmic scheme dedicated to the model reduction problem. By carrying out a complete medium size example, including its qualitative analysis, we have proven that computer algebra tools may be most useful in biological modeling. Embedding these methods in an easy to use package should help practitioners to model cellular processes by means of a much wider variety of functions than the classical Michaelis-Menten or Hill functions, and to make the hypotheses leading to model reductions more explicit than they do. However, the model reduction part of our package still needs many improvements. It also misses tools which make automatic the small set of heuristics which permitted us to carry out the qualitative analysis of our example.

2 Figure 1 .

 21 Figure 1. A single gene regulated by a polymer of its own protein.

Proposition 3 . 1 .

 31 For positive values of the variables and parameters, the reduced model (3.5) exhibits a Poincaré-Andronov-Hopf bifurcation if and only if n ≥ 9.

3. 4 . 1 .

 41 The steady point equations. One computes a Gröbner basis of the ideal generated by the right-hand sides of the differential equations (3.5) w.r.t. the lexicographical ordering G > M > γ 0 . The other variables and parameters are considered as algebraically independent elements of the base field of the equations.

3. 7 .Figure 2 .

 72 Figure 2. A limit cycle between P (t) (horizontally) and G(t) (vertically) obtained by simulating system (3.5) numerically. The parameter values are: n = 15, δ M = 1/2, θ = 1/20, ρ f = 600, ρ b = 1/40, γ 0 = 1/50, α = 6/25 andK i = 1 for each 1 ≤ i < n.The initial values are: P (0) = 1.17, M (0) = 0.8 and G(0) = 0.002.

Proposition 4 . 5 .

 45 With the same notations, if all the Hurwitz determinants c k,0 are positive (0 ≤ k ≤ m -2) and c m-1,0 = 0 and c m-2,1 < 0 then all the eigenvalues of J have negative real parts except a purely imaginary conjugate pair.Proof. The polynomial p m-1 has the special form p m-1 = c m-1,0 ω. We have c m-1,0 = 0. Then p 0 and p 1 have a degree two gcd, p m-2 , which has the special form p m-2 = c m-2,0 ω 2 + c m-2,1 . We have c m-2,1 < 0 and c m-2,0 > 0 thus, the common roots ±ω of p 0 and p 1 are real. Therefore J has one pair of purely imaginary conjugate eigenvalues ±i ω. Now, compute the Sturm sequence (4.6) over the polynomial C(σ) = C(σ)/(σ 2 + ω2 ). This Sturm sequence p0 , p1 , . . . , p m can actually be derived from that of C:p0 (ω) = p 0 σ 2 + ω2 , p1 (ω) = p 1 σ 2 + ω2 , . . . , p m(ω) = c m-2,0 .All the corresponding Hurwitz determinants are positive. According to the Routh Theorem[START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF] Theorem 13.4], all the roots of C have negative real parts. This concludes the proof of the proposition.For m = 3 we have c m-2,1 = -a 3 . For m = 4 we have c m-2,1 = -a 1 a 4 .

More precisely, its approximation M 0 , following the notations of[START_REF] Kokotovic | Singular Perturbation Methods in Control: Analysis and Design[END_REF] Sect. 1.4].

This paper is only concerned by the ordinary case.

In the differential algebra theory, the terminology differential indeterminates is preferred to dependent variables for derivations are abstract and differential indeterminates are not even assumed to correspond to functions. In order not to mix different expressions in this paper, the second expression, which seems to be more widely known, was chosen.

and Pierre-Emmanuel Morant. This work was partially supported by the PPF Bioinfo of the University Lille I.


A Rogani

P Navarro

S Marguet

J-F Ferrero

C Lanouette

Study of post-impact behaviour of thin hybrid carbon/epoxy and glass/epoxy woven composite laminates under fatigue tensile loading -Part II : numerical study

Keywords: Fatigue, Woven composites, thin laminates, hybrid laminates, post-impact behaviour, explicit F.E modelling, semi-continuous approach, damage

This article concerns the modelling of post-impact damage propagation in thin carbon/epoxy and glass/epoxy hybrid woven composite laminates loaded in fatigue tension. The modelling is based on the semi-continuous approach implemented into the explicit finite element code RADIOSS. The bundles are modelled with rod elements and the resin is modeled with a specific damageable shell element, used to stabilize the bundles. This approach is extended to post-impact fatigue tensile loadings. Resin damaging is implemented for the shell elements through two fatigue laws Y xy -N resin a and Y xy -N resin b based on thermodynamical forces and identified through the experimental results obtained in Part I. Fiber breakage is introduced for rod elements through Basquin curves based on tensile and compressive strains.

The damage accumulation is taken into account by a Miner's rule. For several impact energies and amplitude of fatigue loading, the modelling well correlates the experimental results in terms of damage propagation scenario, path and speed.

Introduction

This article deals with the modelling of post-impact damage propagation in thin hybrid carbon/epoxy and glass/epoxy woven composite laminates under fatigue tensile loading. In this paper, the semi-continuous approach, already validated for quasi-static tensile loading on mono-material and hybrid woven composite laminates in [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF], is extended to fatigue loading. It is the second second part of a two-parts article : the phenomenons which occur during the fatigue tensile test have been investigated through several experimental tests. They are described in Part I.

The damage generated by an impact on a woven composite structure is typically composed of matrix cracking, delamination and fibre breakage.

When the structure is subject to fatigue tensile loading, these types of damage propagate, what leads to the final failure of the structure. In order to optimize the post-impact behaviour of their woven composite structure without expensive testings campaigns, the aeronautical industries try to develop numerical tools able to predict the damage propagation.

However, the literature on the subject is barely existing. Most of the studies have been conducted on laminates with an initial defect such as notches, holes or delamination. In these cases, the resin damaging is modeled through a continuum damage modelling (CDM) [START_REF] Miot | Rupture de structures composites stratifiées sous chargements statique et de fatigue[END_REF][START_REF] Hochard | A generalized damage model for woven ply laminates under static and fatigue loading conditions[END_REF]. This approach has been introduced in quasi-static for non-damaged unidirectional composite laminates by Talreja, Ladevèze and LeDantec [START_REF] Talreja | Transverse cracking and stiffness reduction in composite laminates[END_REF][START_REF] Ladevèze | Damage modelling of the elementary ply for laminated composites[END_REF]. The principle is to come degrade the stiffness matrix by one or some damage variables, calculated through thermodynamical forces. Hochard et al. [START_REF] Payan | Damage modelling of laminated carbon/epoxy composites under static and fatigue loadings[END_REF][START_REF] Hochard | A progressive first ply failure model for woven ply cfrp laminates under static and fatigue loads[END_REF][START_REF] Miot | Rupture de structures composites stratifiées sous chargements statique et de fatigue[END_REF][START_REF] Hochard | A generalized damage model for woven ply laminates under static and fatigue loading conditions[END_REF] have adapted this approach to carbon woven composites in quasi-static [START_REF] Payan | Damage modelling of laminated carbon/epoxy composites under static and fatigue loadings[END_REF][START_REF] Hochard | A progressive first ply failure model for woven ply cfrp laminates under static and fatigue loads[END_REF] and in fatigue [START_REF] Miot | Rupture de structures composites stratifiées sous chargements statique et de fatigue[END_REF][START_REF] Hochard | A generalized damage model for woven ply laminates under static and fatigue loading conditions[END_REF].

They have introduced a fatigue damage variable depending on the variation of thermodynamical forces on one fatigue cycle and the number of cycles. Moreover, they have introduced a threshold value corresponding to the minimal thermodynamic force necessary to initiate matrix crackings. A similar modelling (Onera Damage modelling) has been developed by the ONERA, especially for organic matrix and 3D woven [START_REF] Maire | A new formulation of continuum damage me-chanics (cdm) for composite materials[END_REF][START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF]. Their modelling takes into account thermal strains and viscoelasticity. In fatigue, the damage variable includes a static part and a damaging accumulation. For that, they have two modulus : first a static modulus for the static rise which creates a resin damaging and then a fatigue modulus which increases this damaging depending on the number of cycles. Then the variables are updated and the calculation comes back to the static modulus to evaluate the residual performances of the laminate.

To propagate fibers breakages in a glass/woven composite laminate, Rouault [START_REF] Rouault | Propagation de coupure en fatigue sur composites tissés : Etude expérimentale et modélisation[END_REF] based on ε -N curves of its woven fabric identified on fatigue testings on non-impacted samples. He calculates the strains in each element of its modelling. With the curve, he knows the number of cycles necessary to break each element. Then, the principle is delete the element which requires the lower number of cycle to reach the failure. In this case, the damage accumulation is taken into account with a Miner's rule [START_REF] Miner | Cumulative damage in fatigue[END_REF]. Some authors use empirical criterions as the point stress criterion (PSC) and the average stress criterion (ASC). They have been introduced by Whitney and Nuismer [START_REF] Whitney | Stress fracture criteria for laminated composites containing stress concentrations[END_REF]. They observed that, for holed laminate, the maximal stress at the edge of the hole in fatigue can overtake the quasi-static maximal stress without damage propagation. For the PSC, the stress to evaluate the damage propagation is calculated at a distance d 0 f the hole. For the ASC, the stress is a mean stress calculated along a distance d 0 . Tan et al. [START_REF] Tan | Fracture of composite laminates containing cracks due to shear loading[END_REF] used the PSC method to evaluate the stress at failure of samples with an elliptical hole. Miot [START_REF] Miot | Rupture de structures composites stratifiées sous chargements statique et de fatigue[END_REF] used the ASC method for holed and notched carbon laminates. She uses a characteristic volume in order to calculate a mean stress and evaluate the damage propagation. But, even if the implementation of these criterions is simple, they are very dependant on the testings parameters.

Concerning the delamination modelling, cohesive elements are a lot used for impact loading. Their concept have been introduced by Hillerborg et al. [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. A cohesive law links the stress state in the element with the relative displacement of the nodes and a critical energy release rate G c manages the dissipated energy in the cracking opening. Robinson et al. [START_REF] Robinson | Numerical simulation of fatigue-driven delamination using interface elements[END_REF] used them to model the propagation of delamination in fatigue. They represent the delamination propagation in fatigue as a combination of the propagation due to the amplitude of the static loading and due to the cyclic repetition of the load. To do that, they introduce a fatigue damage variable in addition to the damage variable of the classic law. It depends on the relative displacement of the nodes ans the number of cycles. Some authors use the cohesive elements to model the propagation of matrix crackings. For example, Hallett et al. [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF] used them for a tensile simulation on a holed sample. For this purpose, they have first noted the position of the crackings after an experimental test. Then, they meshed the sample placing the cohesive elements on potential paths of the crack propagation.

The objective of the study is to model fatigue tensile tests on impacted woven laminates with the semi-continuous approach [START_REF] Navarro | Semi-continuous approach for the study of impacts on woven composite laminates: Modeling interlaminar behavior with a specific interface element[END_REF][START_REF] Pascal | On the modelling of low to medium velocity impact onto woven composite materials with a 2d semi-continuous approach[END_REF][START_REF] Pascal | Impact damage prediction in thin woven composite laminates -part i: Modeling strategy and validation[END_REF][START_REF] Pascal | Impact damage prediction in thin woven composite laminates -part ii: Application to normal and oblique impacts on sandwich structure[END_REF][START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part i: Experimental study[END_REF][START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF].

For this purpose, the existing semi-continuous strategy has been improved.

In fact, add-ons, based on experimental observations made in Part I, are introduced through propagation laws in fatigue of resin damaging (tows/resin splittings and intra-tows crackings) and fibres breakages. A calculation strategy has been implemented in order to perform their propagation in the explicit finite element code Radioss. The model is calibrated and validated, for several impact energies and levels of solicitation, by simulations on the same woven laminates used in Part I (Tables 12).

Semi-continuous modelling in quasi-static tension

Presentation

In this study, the semi-continuous modelling is extended to fatigue tensile simulation. It has been developed in the research team for impact loadings [START_REF] Navarro | Semi-continuous approach for the study of impacts on woven composite laminates: Modeling interlaminar behavior with a specific interface element[END_REF][START_REF] Pascal | On the modelling of low to medium velocity impact onto woven composite materials with a 2d semi-continuous approach[END_REF][START_REF] Pascal | Impact damage prediction in thin woven composite laminates -part i: Modeling strategy and validation[END_REF][START_REF] Pascal | Impact damage prediction in thin woven composite laminates -part ii: Application to normal and oblique impacts on sandwich structure[END_REF] and extended to post-impact quasi-static tensile loadings on woven laminates [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF] into the explicit finite element code Radioss.

The modelling strategy is to decouple the resin and tows behaviour in order to represent simultaneously the damaged resin and the undamaged fibers tows (Figure1). For this purpose, a specific finite element is developed.

It is made up a truss structure of rod elements, modeling the woven fabric bundles and following the woven pattern geometry, and a shell element, representing the epoxy matrix. An interface element, with an elasting damaging behaviour based on a bilinear cohesive law, connects two plies and represent the delamination.

The rod elements have an elastic and brittle behaviour. Their breakage are calculated through tensile and compressive strain criterions :

If ε(x) > ε tension max or |ε(x)| > |ε comp max | then break of the bar element (1)
Furthermore, the cross-point of the woven ply are taken into account for these criterions. Indeed, at a cross-point, the fibres bundles are in local bending.

So, they need a more (resp. less) important tensile (resp. compressive) strain to spread out and reach the tensile (resp. compressive) strain criterion. So, at a cross-point, ε tension max is increased and |ε comp max | is reduced. For this purpose, ε comp max cross and ε tension max cross are introduced in the model such as :

|ε comp max cross | |ε comp max | = ε tension max ε tension max cross (2)
The resin has an elastic and isotropic behaviour with an anisotropic damaging in each direction of the woven through two damage variables d 1 and d 2 , based on thermodynamical forces (or energy release rates). It has also a pseudo-plastic behaviour in plane shear with a shear damaging law. It is calculated through two damage variables d a xy and d b xy which represent respectively tows/resin splittings and cracks into bundles [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF]. This original shear damaging law is based on experimental observations made by [START_REF] Karahan | Investigation of damage initiation and propagation in 2x2 twill woven carbon/epoxy multi-layer composites[END_REF][START_REF] Bogdanovich | Quasi-static tensile behavior and damage of carbon/epoxy composite reinforced with 3d non-crimp orthogonal woven fabric[END_REF] during quasi-static tensile tests on carbon/epoxy woven laminate oriented at ±45 • : first there is an initiation of tows/resin splitting and then, intra-tows crackings initiate and propagate in meta-delamination (delamination between warp and weft bundles). The splittings quickly evolve at the beginning of the test then they stabilize. The crackings and the meta-delaminations slowly grow at the beginning and suddenly accelerate once the splittings no longer evolve. However, the modelling scale is not enough small to represent meta-delamination. The two damage variables depend on thermodynamical function, as it is commonly used [START_REF] Hochard | A progressive first ply failure model for woven ply cfrp laminates under static and fatigue loads[END_REF][START_REF] Hochard | A generalized damage model for woven ply laminates under static and fatigue loading conditions[END_REF] :

             d a xy = d a max xy × 1 -exp - Y xy - √ Y 0 ≥0 √ Y c d b xy = max exp Y xy - √ Y 0 bis ≥0 √ Y c bis , 1 (3)
where Y 0 (respectively Y 0 bis ) and Y c (respectively Y c bis ) are parameters of damage threshold and damage speed of d xy a (respectively d xy b ). d max xy a is the maximal limit of d xy a . This limit has been introduced for tows/resin splittings because they are not responsible for the final break of the laminate.

The coupling of the shear stress with the stresses in the other directions is neglected. So, shear strain energy is only due to shear stress and the thermodynamical function Y xy is defined by :

Y xy = ∂W xy ∂d xy (4)
where d xy = d a xy + d b xy and W xy is the shear strain energy defined by :

W xy = τ 2 xy 2G 0 xy (1 -d xy ) (5)
Hence :

Y xy = 1 2 G 0 xy γ e xy 2 (6)
Finally, the two shear damage variables are used to degrade the shell and the bending matrix :

               (1 -d 1 ) 1 -d a xy -d b xy E 1 1 -(1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy νxy ν 21 (1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy ν 21 E 1 1 -(1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy νxy ν 21 0 (1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy νxy E 2 1 -(1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy νxy ν 21 (1 -d 2 ) 1 -d a xy -d b xy E 2 1 -(1 -d 1 )(1 -d 2 ) 1 -d a xy -d b xy νxy ν 21 0 0 0 1 -d a xy -d b xy G 0 xy               
The irreversible strains are taken into account by the introduction of a pseudo-plastic criterion on the in-plane shear stress based on an isotropic hardening using a classic elastic prediction and a plastic correction. An elastic field is defined by :

f = |τ xy | -K plas p β -τ 0 (7)
where τ xy is the in-plane shear stress, τ 0 is the plastic strength, (K plas , β) are parameters of the plastic hardening law and p the cumulative plastic strain.

Thus, if f > 0, the plastic correction is performed through the update of the plastic variables γ p xy and p, based on a Newton-Raphson iterative scheme, with γ p xy the plastic shear strain defined by :

γ xy = γ e xy + γ p xy (8)
where γ xy is the total shear strain and γ e xy is the elastic shear strain.

Validation

The semi-continuous modelling has been validated in quasi-static for all the mono-material and hybrid woven laminates presented in the Tables 12[1]. The model well represents the propagation of the post-impact damage.

The fracture surfaces and the load -crosshead displacement curves well correlate the experimental results. They are expressed as Basquin curves :

Propagation laws in fatigue

N resin a = Y a Y xy -Y 0 a 1 α and N resin b = Y b Y xy -Y 0 b 1 β (9)
where Y a , α, Y b and β are parameters of the fatigue laws. Threshold values Y 0 a and Y 0 b have also been added for each damage type [START_REF] Miot | Rupture de structures composites stratifiées sous chargements statique et de fatigue[END_REF][START_REF] Hochard | A generalized damage model for woven ply laminates under static and fatigue loading conditions[END_REF]. These parameters will be identified through severeal post-impact fatigue tensile simulations on laminates C45C45 and C0C0.

Fibres bundles breakages

The fatigue tensile loadings on impacted samples of configuration G0G0, G0C45C45 and G0C45C0, performed in Part I, have shown fibres bundles breakages in carbon and glass plies. They are taken into account in the semicontinuous modelling through the implementation of two Basquin curves (for each material) for the rod elements based on their tensile strain :

N breakage bundle = ε tension 6 ε bundle + 1 γ (10)
where ε bundle + is the tensile strain of the rod element during the simulation. ε 6 tension , the tensile strain at 10 6 cycles, and γ are the law parameters. They have been identified by Airbus Helicopters by tensile fatigue testings piloted in strain on non-impacted samples of configuration C0C0 and G0G0. In fact, in tension, the sample behaviour is governed by the fibres.

The tensile loadings on impacted samples of configuration G0C45C45, performed in [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part i: Experimental study[END_REF] and Part II, indicate a similar behaviour in quasi-static and fatigue. In the glass ply, the damage propagation is due to tensile and compressive breakage of warp and weft bundles [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part i: Experimental study[END_REF][START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF]. In fact, the glass ply elongation is piloted by the two carbon plies elongation. So, the weft bundles of the glass ply are compressed due to the Poisson's effect of the carbon plies oriented at ±45 • . In the quasi-static modelling, this phenomenon is represented by a compressive strain criterion of the rod element (see 2.1).

In fatigue, this phenomenon is represented by the introduction of a Basquin curve for the rod elements based on their compressive strain :

N breakage bundle = |ε comp 6 | | ε bundle -| 1 γ (11)
where ε bundle -is the compressive strain of the rod element during the simulation. The value of γ is the same as in tension.

ε comp 6 is chosen such as the ratio |ε comp 6 |/ε tension 6
is the same as in quasi-static :

|ε comp 6 | ε tension 6 = |ε comp max | ε tension max (12)
The cross-points are also taken into account in fatigue. In fact, the experiments of Part I have indicated than the propagation of glass fibres breakages in the impacted laminates of configuration G0G0 and G0C45C45 was guided by the cross-points of the woven plies. In tension and compression, two different curves are used whether the fibre is located at a cross-point or not such as :

ε tension 6 ε tension 6 cross = |ε comp 6 cross | |ε comp 6 | = ε tension max ε tension max cross (13)

Damage accumulation

The damage accumulation in the shell and rod elements during the fatigue tensile simulation is take into account by a simple Miner's rule. For each rod element and shell element of the specific woven finite element (Figure 1), a variable D m is introduced such as :

D m = i n i N i (14)
For a rod element, n i is the number of cycles realized at a strain ε i for which the number of cycles at failure is N i . For a shell element, n i is the number of cycles realized at an energy release rate Y i for which the number of cycles at failure is N i . Thus, if D m = 1, the rod or the shell element is deleted.

Strategy of post-impact fatigue tensile simulation

Principe

The principle of the simulation strategy is to find in the model the rod or the shell element which requires the lowest number of cycles to break. That amounts to find the number of cycles N elem which verifies :

N elem = min j∈mesh (1 -D m (j)) * ε 6 ε j 1 γ , (1 -D m (j)) * Y i Y j -Y 0 i 1 k (15)
with i = a or b and i = α or β. In this part, the propagation laws presented above are identified and validated on the mono-material laminates. For the tows/resin splittings and the intra-tows crackings, they are identified and validated through simulations on samples of configuration C45C45 and C0C0. For glass fibres breakages, they are validated through simulations on samples of configuration G0G0.

Strategy

The objective is to obtain numerical results which well correlate the experimental results in terms of damage propagation scenario, damage propagation velocity and number of cycles at failure. In order to test the robustness and the validity of the modelling, the same impact energies and displacement levels as in the experimental study are used for the simulations. They are presented on the Table 3.

Configurations C0C0 and C45C45

The propagation laws in fatigue Y xy -N cycles , identified for the tows/resin splittings and the intra-tows crackings, are presented on the Figure 2.

Configuration C45C45.

For the configuration C45C45, the post-impact damage scenario is presented on the Figure 3. As for the experimental testings, the damage propagates at ±45 • . Tows/resin splittings initiate from each tip of the initial damage transversely to the initial fibres breakages. They join into a unique point and form a rhombus-shape damage. Then, they propagate at ±45 • from this point. Once this propagation is started, intra-tows crackings initiate also and propagate at ±45 • following the tows/resin splittings. Their propagation continue with a slower progression of the intra-tows crackings.

Finally, the splittings reach the edges of the samples following some cycles later by the intra-tows crackings. The samples is entirely cracked on a ±45 • band. The numbers of cycles required to reach the edges of samples obtained with the modelling are very close to the values experimentally obtained (Tables 4-5).

For the couples (∆L 2 ; 2.25 J) and (-20% ∆L 2 ; 4 J), the post-impact damage barely propagates (Figure 4). In fact, the values of energy release rates calculated in the sample are lower than the value required to initiate intra-tows crackings (Figure 5). So, they need a high number of cycles to emerge.

Configuration C0C0.

For the laminate C0C0, the post-impact damage propagation scenario obtained with the modelling is in accordance with the experimental results obtained in Part I. Tows/resin splittings and then intra-tows crackings propagate vertically during the fatigue simulation (Figure 6). Furthermore, the velocity of damage propagation numerically obtained well correlates the velocity experimentally obtained. The evolution curves depending on the number of cycles are presented on the Figure 7 for the tows/resin splittings and then intra-tows crackings. They are close to the experimental curves for all the energies and displacement levels tested.

They propagate out of step very quickly in the 100 000 -200 000 first cycles before stabilizing progressively. This stabilization is due to the decrease of the shear strains in the resin at the tip of the damage propagation as the simulation progresses. That leads to a decrease of the energy release rates in the resin and an increase of the number of cycles required to initiate the splittings. Finally, for 10 7 cycles, the energy release rates calculated at the tip of the splittings are lower to the damages threshold values Y a 0 and Y b 0 . So, the damage no longer evolves For the highest displacement level (+30% ∆L 3 ), the fibres reach their quasi-static tensile strain criterion during the quasi-static rise. That leads to an horizontal propagation of fibres breakages as in quasi-static (Figure8) [START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF].

Configuration G0G0

For the configuration G0G0, the fracture surface numerically obtained for each impact energy and displacement level is in accordance with the experimental results (Figure 9). Fibres breakages propagate horizontally following the cross-points of the woven fabric (Figure 10).

The damage propagation initiation and speed are close to those experimentally obtained. Indeed, the experimental and numerical curves of damage length depending on the number of cycles have the same appearance (Figure 11) and the number of cycles at the sample failure are close.

Conclusion

The semi-continuous modelling is validated for the mono-material laminates. For C45C45, C0C0 and G0G0, the model well correlates the experimental results of Part I in terms of damage propagation scenario, damage propagation velocity and number of cycles at failure for several impact energies and displacement levels.

In the next part, the semi-continuous modelling is tested on the hybrid configuration G0C45C45.

Validation on hybrid laminate G0C45C45

As for the mono-material configurations, several impact energies and displacement levels are tested for G0C45C45. The values are presented in the Table 6.

For each impact energy and displacement level, the fracture surfaces are in accordance with the experimental results of Part I (Figure 12). Fibres breakages propagate at ±45 • following the cross-points of the upper glass ply (Figure 13) and and resin damaging propagates at ±45 • in the lower carbon ply.

The post-impact damage propagation scenario is precisely obtained (Figure 14) . The simulation brings some elements not visible during the experimental tensile tests in Part I. Fibres breakages initiate in the upper glass ply while no resin damaging has emerged in the two carbon plies. Warp and weft bundles are highly loaded respectively in tension and compression at the tip of the damage propagation. That leads to a propagation a fibres breakages in tension for warp bundles and in compression for weft bundles. They propagate at ±45 • following the cross-points of the glass ply along an area highly loaded in plane-shear in the carbon plies. In fact, the carbon plies elongation pilots the glass ply elongation. Fibres breakages propagate gradually in the glass ply and, out of step, resin damaging propagates also in the carbon plies. From a certain number of cycles, the splittings in carbon plies overtakes the fibers breakages in the glass ply and lead their propagation.

When the propagation comes closer to the edge of the sample, intra-tows crackings emerge in the carbon plies and each damage quickly propagate in the three plies until the edges of the sample and the break of the glass ply. However, once the glass ply is broken, the carbon plies are cracked but not broken. Indeed, the carbon fibres are not enough loaded. The resin keeps damaging along the band already cracked.

The fibres breakages propagate at a velocity close to the experimental velocity noted in Part I. The curves of fibres breakages length depending on the number of cycles are similar whatever the initial damage size (impact energy) or the imposed displacement (Figure 15). They progress into three steps : an initiation of the propagation, a propagation quite slow and a fast propagation when the damage comes closer to the edges of the sample, due to a structural effect.

The number of cycles required to initiate the damage propagation in the glass ply and the break of the glass ply with the modelling are also close to values experimentally noted (Tables 78).

The post-impact damage propagation scenario explains the influence of the impact energy and the displacement level on the damage initiation and propagation. In fact, the propagation starts with fibres breakages in tension and compression in the glass ply. With a greater initial damage or higher imposed displacement, the glass fibres at the edges of the impacted area deform more easily in tension and compression (Figure 16).

As for the matrix damaging of the carbon plies overtakes the fibers breakages of the glass ply from a certain number of cycles, the final fast evolution of the damage, when it comes closer to the edge of the sample, is lead by the matrix damaging.

Limitation of the modelling

The semi-continuous modelling is validated for fatigue tensile simulations on carbon/epoxy woven and glass/epoxy woven laminates made up of two identical plies and on the hybrid woven laminate G0C45C45 for which the two carbon plies have the same orientation and no delamination is observed during the experimental fatigue tests. However, a limitation is reached when the two plies have a different orientation and delamination is observed during the experimental tests. For instance, the Figure 17 shows the result of simulations performed on a sample of configuration G0C45C0. First, the modelling manages to obtain two different scenario according the impact energy and the displacement level as for the experimental observation. Even the laminate does not break (after 7.10 6 cycles) : there are only fibers breakages in the glass ply and vertical matrix damaging the lower carbon ply. Even, with the increase of impact energy and displacement level, propagation of fibres occurs simultaneously in three plies and leads to a rapid failure of the laminate. But, in the first case, the matrix damaging in the carbon ply is not enough developed in comparison to the experimental damaging. This difference is due to the lack of delamination, which has not been implemented yet in the modelling, contrary to what has been done for example in [START_REF] Robinson | Numerical simulation of fatigue-driven delamination using interface elements[END_REF].

Nevertheless, the modelling brings an explanation concerning the change of scenario. When the impact energy and the displacement level increases, the fibers of the lower carbon plies are more loaded and fibres breakages can emerge (Figure 18). So, at the experimental level, when the impact energy and the displacement level increase, the fibers breakages propagation can take over the delamination, what leads the failure of the laminate. In the other case, the propagation of delamination between the two carbon plies unload the carbon fibers. So, the resin damaging takes over the fibers breakages.

Conclusion

In this article, post-impact fatigue tensile tests have been modeled with the semi-continuous strategy. The modelling manages to represent the propagation of resin damage (tows/resin splittings and intra-tows crackings) in the carbon resin and the propagation of bundles breakages in the glass ply during the fatigue tensile simulation.

The main contribution of this article is to provide a description of a modelling strategy that is able to represent the damage induced by an impact and also the propagation of this damage in fatigue loading with the same model.

The original scale and the choice of the damage laws based on experimental observation and on the litterature allows to well represent the physics, the specific role of the resin and that of the fibres in the behaviour.

In order to do that, the semi-continuous modelling has been extended to fatigue tensile simulations on impacted thin woven carbon/epoxy laminates.

3. 1 .

 1 Resin damagingThe fatigue tensile loadings on impacted samples of configuration C45C45 and C0C0, performed in Part I, have shown that the post-impact damage propagation was governed by shear damaging : tows/resin splittings and intra-tows crackings propagate in the laminate at ±45 • for the configuration C45C45 and vertically for the configuration C0C0. In both cases, the two resin damaging types follow the same direction and have a similar behaviour.But, the splittings emerge first and propagate faster than the crackings.As for the quasi-static modelling, the propagation in fatigue of each type of resin damaging is based on thermodynamical forces. Two fatigue laws Y xy -N resin a and Y xy -N resin b are implemented for the shell element, only for carbon material. N resin a and N resin b are the number of cycles necessary to initiate respectively the tows/resin splittings and the intra-tows crackings.

Finally, the

  post-impact fatigue tensile simulation strategy is implemented in FORTRAN through a user law in the explicit finite element code Radioss. The user law is made up of a quasi-static and a fatigue modules. The steps of the simulation are : 1. A quasi-static rise is numerically realized until the maximum displacement of the experimental fatigue cycle : ∆L = ∆L sta + ∆L dyn . A first calculation of the strains and the thermodynamical forces are performed for each element of the model. 2. Then, the calculation goes into the fatigue module : (a) First, the breakage number of cycles à rupture is calculated for each rod element and shell element (only for carbon) of the model. (b) If the lowest number of cycles is reached for :a rod element : then it is deleted of the model as well as the shell element. Indeed, it assumed that the bundle breakage lead to the breakage of the resin in the bundle and at the interface between the bundle and the resin. a shell element : either the damage corresponds to tows/resin splitting so the shell element does not have to be entirely deleted because the intra-tows crackings did not appear yet.The damage variable of the shell element is set to 0.5. Or the damage corresponds to intra-tows crackings and the shell element is entirely deleted.(c) Once the rod or the shell element is damaged, the number of cycles and all the variables D m (in order to simulate the damage accumulation) are uploaded.

3 .

 3 Before moving to the next iteration, the calculation goes again into the quasi-static module in order to get the new strains and thermodynamical forces of each element of the model.5. Fatigue tensile simulations on impacted samples5.1. Identification and validation on mono-material laminates 5.1.1. Introduction

  The resin damaging is implemented for the shell elements through two fatigue laws Y xy -N resin a and Y xy -N resin b based on thermodynamical forces. Fiber breakage is introduced through Basquin curves based on tensile and compressive strains for rod elements. Damage accumulation is taken into account by a simple Miner's rule. The strategy is to find in the sample the shell or rod element which requires the lowest number of cycles to reach the breakage and delete it. The resin fatigue laws have been identified through several simulations on impacted samples of configuration C45C45 and C0C0. The tensile and compressive fiber breakage Basquin curves have been identified through simulations respectively on G0G0 and G0C45C45. The modelling is validated in terms of damage propagation scenario, number of cycles at initiation and speed for mono-material carbon/epoxy and glass/epoxy laminates made up of plies with the same orientation. It also has been validated for the hybrid configuration G0C45C45. In this case, the modelling can highlight the precise scenario of the post-impact damage propagation. It shows especially that, from a certain number of cycles, the matrix damaging of the carbon plies overtakes the fibers breakages of the glass ply and leads the global propagation of the damage. However, as the delamination has not been implemented yet, the semicontinuous modelling has not been validated for the configuration G0C45C0, what shows the importance of the delamination consideration in the mod-elling to correctly represent the post-impact damage propagation. Nevertheless, two different propagation scenario are obtained according the impact energy or the level displacement as for the experimental tests in Part I. When the impact energy and the displacement level increase, fibers breakages can appear in the lower carbon ply oriented at (0/90) • . Their propagation can take over the delamination, what leads the failure of the laminate. In the other case, the propagation of delamination between the two carbon plies unload the carbon fibers and the resin damaging is predominant.

Figure 1 :

 1 Figure1: Modelling strategy of the woven ply[START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF]

Figure 2 :

 2 Figure 2: Resin damaging laws identified with fatigue simulations on impacted samples of configuration C45C45 and C0C0

Figure 3 : 29 Figure 4 :

 3294 Figure 3: Post-impact damage propagation scenario numerically obtained for the configuration C45C45

Figure 5 :

 5 Figure 5: Intra-tows crackings propagation depending on the energy release rate Y xy calculated in the resin for the configuration C45C45

Figure 6 :

 6 Figure 6: Post-impact damage propagation scenario for the configuration C0C0

Figure 7 :Figure 8 :

 78 Figure 7: Length evolution of the tows/resin splittings (Damage 1) and the intra-tows crackings (Damage 2) depending on the number of cycles for the configuration C0C0

Figure 9 :Figure 10 :

 910 Figure 9: Comparison numerical/experimental of the fracture surfaces obtained after fatigue tension on impacted samples of configuration G0G0

Figure 11 :

 11 Figure 11: Length evolution of fibres breakages depending on the number of cycles for the configuration G0G0

Figure 12 :

 12 Figure 12: Numerical/experimental comparison of the fracture surfaces obtained after fatigue tension on impacted samples of configuration G0C45C45

Figure 13 :Figure 14 :

 1314 Figure 13: Propagation in fatigue of fibres breakages following the cross-points in the upper glass ply of the configuration G0C45C45

Figure 15 :Figure 16 :Figure 17 :Figure 18 :

 15161718 Figure 15: Length evolution of fibres breakages depending on the number of cycles in the upper glass ply of the configuration G0C45C45

Table 1 :

 1 Configurations for single-material laminates[START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF]

	Configuration	Representation	UpperLply	MiddleLply	LowerLply
	G0C45C45		lassL(0/90)°CarbonLLL45°+	arbonLLL45°+
	G0C45C0		GlassL(0/90)°G	CarbonLLL45°+	CarbonL(0/90)°C

Table 2 :

 2 Configurations for hybrid laminates[START_REF] Rogani | Tensile post-impact behaviour of thin carbon/epoxy and glass/epoxy hybrid woven laminates -part ii: Numerical study[END_REF]

Table 3 :

 3 Impact energies and displacement levels used for the fatigue simulations on C45C45, C0C0 and G0G0 E impact N edges EXP N edges NUM

	3 J	120 000	100 000
	4 J	60 000	50 000
	6.25 J	15 000	8 000
	9 J	20 000	10 000

Table 4 :

 4 Number of cycles necessary to reach the edges of the sample depending on the impact energy for the configuration C45C45

Acknowledgment

This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2019-P09105.

Displacement N edges EXP N edges NUM ∆L 2 60 000 50 000 +10% ∆L 2 15 000 12 000