Xiao Ma

Yang Chen

Modesar Shakoor
email: modesar.shakoor@imt-nord-europe.fr

Dmytro Vasiukov

Stepan Lomov

Chung Hae Park

Simplified and complete phase-field fracture formulations for heterogeneous materials and their solution using a Fast Fourier Transform based numerical method

Keywords: Fast Fourier Transform (FFT), Phase-field, Damage modeling, Heterogeneous material

This paper focuses on the numerical implementation of phase-field models of fracture using the Fast Fourier Transform based numerical method. Recent developments in that field rely on the separate solution of a coupled problem where the mechanical equilibrium problem is solved first, and then the phase-field evolution equation. The latter involves a diffusion term which has been simplified in previous

Introduction

A good comprehension of the failure mechanics of materials requires a deep study on the microscale, especially for heterogeneous materials. One of the study directions of cracks in heterogeneous materials is fracture mechanics, which principally focuses on the crack initiation and propagation. This subject has been extensively studied since the milestone works of Griffith [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids[END_REF] and Irwin [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF]. In the work of Griffith, a global energy approach is provided to find quantitative relations between the crack length, the resistance to the crack growth of the materials, and the criterion of crack propagation that regards fracture as a competition between the surface energy of propagation and the elastic energy stored in the bulk material.

If we look at the microscale, right before the material attains its strength, the macro-crack does not propagate directly. First, there is an accumulation and evolution of local microdefects, such as micro-cracks and micro-cavities. Due to these micro-defects, the surface for bearing the load reduces. Hence, the material cannot carry the load as before even if the effective macroscopic load increases. From a global point of view, this phenomenon can be regarded as the degradation of local material properties. The degree of material degradation can be described by a continuous parameter d called the damage variable. This variable d varies in the range [0, 1], where 0 means intact state and 1 means fully degraded state.

Phase-field models, one of the branches of CDMs, also utilize a damage variable d to describe the failure state. In the phase-field approach, sharp cracks are regularized as diffusive crack bands by a function of d [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]. On the one hand, the evolution of the phase-field variable itself completely describes the initiation and propagation of cracks in the material. On the other hand, phase-field models are closely related to the variational approach to brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. This approach enables a simultaneous approximation of the displacement field and the cracks by minimizing the total potential energy of the solid.

Comparisons between gradient-based non-local damage models and the phase-field approach can be found in the literature [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF][START_REF] Mandal | Phase field and gradient enhanced damage models for quasi-brittle failure: A numerical comparative study[END_REF]. These studies show that the two approaches are similiar in terms of equations and solutions. Even some gradient-based damage models have been proposed that tend to be closer to the phase-field model [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. That is why some researchers consider the phase-field model as a gradient-based non-local damage model. However, as discussed in Refs. [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF][START_REF] Wu | Phase-field modeling of fracture[END_REF], their fundamental ideas are different. The gradient-based non-local damage model was proposed to average the local damage field to solve the mesh dependency, while the phase-field model departs from the discontinuous description of a crack, where a distribution function is introduced to smear this discontinuity. For example, the regularization from sharp crack to diffusive crack is not introduced in Ref. [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF] even if the final equations are almost the same as in the phase-field approach. More details can be found in the review paper [START_REF] Wu | Phase-field modeling of fracture[END_REF].

Compact formulations of the phase-field model can be found in Refs. [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. These works and others have promoted the variational phase-field models in the engineering community. According to the literature [START_REF] Wu | Phase-field modeling of fracture[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF], the major characteristics of the phasefield model can be summarized as follows:

• The phase-field model is purely based on energy minimization, and there is no need for a pre-defined crack.

• Multi-cracks merging and branching problems can be solved without additional effort.

• The model can easily be applied in 3D cases, and numerical implementation is straightforward.

• The model is non-local and intrinsically avoids mesh dependency.

The phase-field approach, hence, has gained increased attention among the damage modeling techniques. Many applications of this approach can be found in the recent literature.

Espadas-Escalante analyzed the fracture behavior of fiber-reinforced composites on the microscopic scale [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF]. In Ref. [START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF], a phase-field model has been applied to study cracks in concrete materials. For a multi-phase system, phase-field models can be combined with other phase-field [START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF] or cohesive zone models [START_REF] Chen | Fft phase-field model combined with cohesive composite voxels for fracture of composite materials with interfaces[END_REF] to study the interface debonding problem.

Although numerical methods based on the phase-field approach is generally considered mesh independent, a condition on element size must be ensured. It has been proven in Ref. [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF] that the characteristic length of the phase-field model needs to be twice larger than the element size of the fractured zone. In this context, especially for 3D simulations, conventional finite element solvers can be cumbersome and hence require an efficient parallel implementation to extend the computational limits [START_REF] Nguyen | Largescale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. There is a demand for massive parallelization in complex microstructure cases to get a reasonable computational time.

The Fast Fourier Transform (FFT)-based method was initially proposed by Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] in 1994 as a voxel-based methodology that does not need stiffness matrix assembling, unlike conventional finite element solvers. In the FFT-based method, the local strain tensor is calculated by a convolution product with a fourth-rank Continuous Green Operator (CGO) and a polarization term. The convolution product is transformed into a simple multiplication in Fourier space. The multiplication of the polarization term and the CGO are local operations that can be easily parallelized [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF]. The FFT operation itself is not a local operation but efficient parallel implementations are available (such as FFTW [START_REF]FFTW[END_REF]). Consequently, large-scale simulations based on full-resolution images can be performed with the FFT-based method. Since its proposition, FFT-based methods have been improved and applied to investigate a wide range of physical phenomena in heterogeneous materials [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF][START_REF] Ma | Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods[END_REF].

As it was shown in Refs. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF][START_REF] Ernesti | Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures[END_REF], the combination of the FFT method with the phasefield model for fracture problems can be attractive. However, a fundamental issue must be taken into account for heterogeneous materials, which is the main subject discussed in this paper.

When studying fracture problems, a model can often be categorized by the degree of isotropy and homogeneity. Isotropic models, in our studies, have two main families of properties that can help to describe the degree of homogeneity: elastic properties, such as Young's modulus (E), and damage properties, such as critical energy release rate (G c ). In our recent studies, if the damage properties are homogeneous, the solver in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] works well no matter the homogeneity degree on elastic properties. However, if the damage properties become heterogeneous, the results often show non-physical damage diffusion between phases. This is the main topic that will be discussed in the present paper, and the word heterogeneous herein is mainly linked to the damage properties.

Based on the paper of Jeulin [START_REF] Jeulin | Towards crack paths simulations in media with a random fracture energy[END_REF], due to the fact that the G c values are not constant, one term in the phase-field equation has been omitted in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] for heterogeneous materials. Although various works have adopted a heterogeneous setting for material fracture properties, it is still not fully clear whether such heterogeneity can be justified, for which rigorous discussions are seldom in the literature. A recent work [START_REF] Vicentini | Phase-field modeling of brittle fracture in heterogeneous bars[END_REF] discussed this aspect using 1D heterogeneous bars.

The present work focuses on Miehe's initial formulation of the phase-field approach [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], although the readers are to be aware that numerous improvements and extensions have been proposed in the later literature. For the sake of brevity, in the subsequent text, the algorithm used in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] will be called the simplified formulation, and the modified algorithm after adding the heterogeneity correction will be called the complete formulation. Note that a recent paper [START_REF] Ernesti | Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures[END_REF] utilizes a multi-phase-field system for heterogeneous mediums, which may also be able to solve the mentioned diffusion issue, but it is not the scope of this paper.

The paper is organized as follows: first, a brief review of Miehe's phase-field model, on which this paper will focus, is presented in Sec. 2. Then, the implementation of the complete formulation in the FFT-based solver is presented in Sec. 3. Third, several numerical experiments are carried out in Sec. 4 to compare the difference between simplified and complete formulations. Finally, Sec. 5 presents general conclusions.

Review of Miehe's phase-field model

The starting point of Miehe's phase-field model [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF] is smearing the sharp crack topology by a diffusive crack band, as illustrated in Fig. 1. The crack is smeared using the exponential function

d (x) = exp - |x| l c . ( 1
)
The solid Ω with a crack set Γ (or Γ l ) is presented in Fig. 1. The conventional sharp crack is shown in Fig. 1(a) and the diffusive crack used in phase-field is shown in Fig. 1(b).

The damage variable d represents the local damage state with d = 1 for fully damaged material and d = 0 for the intact state.

Since the model in Fig. 1(a) possesses a sharp crack, we will get the delta function shown in the yellow window if we draw the evolution of d along the AB line. This contrasts with the diffusive crack presented in Fig. 1(b), where the crack is approximated as a crack band characterized by a length parameter l c .

The influence of l c values on the damage band diffusion is presented in Fig. 1(c). The smaller l c is, the narrower the crack band is. The curves in Fig. 1(c) all satisfy that d = 1 Looking at Miehe's formulation, in a 1D case, it can be easily verified that the diffusive crack function in Eq. ( 1) is a solution to the following ordinary differential equation

(ODE) 1 l c d (x) -l c d (x) = 0, (2)
subjected to the Dirichlet-type boundary conditions: d (x = 0) = 1 and d (x = ±∞) = 0.

On the one side, in consistence with the variational approach to fracture, the crack surface function can be obtained from the minimization of the crack surface that yields the ODE shown in Eq. (2):

d (x) = arg {inf [Γ l ]} , with Γ l expressed as Γ l (d) = +∞ -∞ 1 2 1 l c d 2 + l c (d ) 2 dx = +∞ -∞ γdx,
where γ is the crack surface density. The extension to a higher dimension is straightforward:

γ (d, ∇d) = 1 2 1 l c d 2 + l c (∇d) 2 , ( 3
)
where d has been replaced by ∇d.

On the other side, the phase-field is also in line with Griffith's theory [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids[END_REF]. In this energetic approach, the damage evolution and crack propagation result from the competition between the bulk energy stored in the body and the energy dissipation from the opening of the crack surface. From this point of view, in quasi-static loading, the total energy Π can be expressed as:

Π := Φ s + Φ d -P, ( 4
)
where Φ s is the strain energy stored in the cracked body, Φ d is the energy dissipated for opening the crack surface, and P is the external loading. As illustrated in Fig. 1(c), we consider a domain Ω ⊂ R n (n = 1, 2, 3) which contains a diffusive crack set Γ l with Γ l ⊂ R n-1 and is under a volumetric loading by a body force f * , a displacement u * on the boundary ∂Ω u , and a surface force t * on the complementary boundary ∂Ω t (∂Ω = ∂Ω u ∂Ω t ), where the symbol () * means a prescribed term. A detailed formulation of the components in Eq. ( 4) is presented as follows:

Φ s = Ω ϕ (ε, Γ l ) dΩ = Ω g (d) ϕ 0 (ε) dΩ, (5a)
Φ d = Γ l G c d∂Ω ≈ Ω G c γdΩ, ( 5b
) P = Ω f * • udΩ + ∂Ωt t * • ud∂Ω + ∂Ωu (σ • n) • u * d∂Ω, ( 5c
)
where ϕ is the elastic strain energy density stored in the cracked body, g is the energetic degradation function, ϕ 0 is the initial strain energy, ε is the strain tensor, G c is the critical energy release rate, σ is the stress tensor and n is the outward normal vector.

The displacement field and the crack set can be solved by minimizing the total energy with the variational approach [START_REF] Bourdin | The variational approach to fracture[END_REF]:

(u (t) , Γ l (t)) = arg {min [Π (u, Γ l )]} , ( 6
)
under the constraint of irreversibility, which is stated as Γl ≥ 0, and under the boundary conditions: u(x) = u * at ∂Ω u , and σ • n = t * at ∂Ω t . Assuming a zero body force, the detailed expression of Eq. ( 4) is

Π (u, d) = Ω ϕ (ε (u) , d) dΩ+ Ω G c γ (d, ∇d) dΩ- ∂Ωt t * •ud∂Ω- ∂Ωu (σ •n)•u * d∂Ω. ( 7
)
Eq. ( 7) is then rewritten in the variational form for determining the displacement and phase-field (u, d):

δΠ (u, d) = Ω σ : δεdΩ + Ω g (d)δdϕ 0 (ε) dΩ + Ω G c ∂γ ∂d δd + ∂γ ∂∇d .δ∇d dΩ -∂Ωt t * • δud∂Ω, (8)
where σ = ∂ϕ ∂ε denotes the Cauchy stress, and δu = 0 for ∀x ∈ ∂Ω u . In general, Eq. ( 8) can be divided into two parts δΠ (u, d) P 1 and δΠ (u, d) P 2 . The first one is the mechanical part without damage that is written as follows:

δΠ (u, d) P 1 = Ω σ : δεdΩ - ∂Ωt t * • δud∂Ω.
Regarding δΠ (u, d) P 2 , it can be written as:

δΠ (u, d) P 2 = Ω g (d) δdϕ 0 (ε) dΩ + Ω G c l c dδddΩ + Ω G c l c ∇d.∇δddΩ. (9)
Applying the divergence theorem, Eq. ( 9) would yield

δΠ (u, d) P 2 = Ω g (d) δdϕ 0 (ε) dΩ + Ω G c l c dδddΩ - Ω ∇ • (G c l c ∇d) δddΩ+ ∂Ω (G c l c ∇dδd) • nd∂Ω. ( 10
)
The surface integral vanishes by considering the crack boundary condition: ∇d (x)•n (x) = 0. Thus, we have:

δΠ (u, d) P 2 = Ω g (d) δdϕ 0 (ε) dΩ + Ω G c l c dδddΩ - Ω ∇ • (G c l c ∇d) δddΩ. ( 11
)
Since the solution of the crack set, the damage variable herein, is a stationary point of the energy functional, the local damage evolution equation of the phase-field can be written in strong form as:

g (d) ϕ 0 (ε) + G c l c d -∇ • (G c l c ∇d) = 0, (12)
when δd > 0. In the case of homogeneous medium, the G c and l c are constant. Thus,

∇ • (G c l c ∇d) = G c l c
∆d, and the damage evolution function can be rewritten as:

g (d) ϕ 0 (ε) + G c l c d -G c l c ∆d = 0, (13)
In Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], Eq. ( 13) was applied. However, in the case of heterogeneous medium, G c and l c are functions of x for x ∈ Ω. Thus, these two parameters cannot be directly put outside the divergence operation, and Eq. ( 12) should be adopted. Comparing Eq. ( 13)

with Eq. ( 12),

∇ • (G c l c ∇d) = ∇ (G c l c ) • ∇d + G c l c ∆d. ( 14
)
That means the term (∇ (G c l c ) • ∇d) was omitted in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]. This is the reason why the phase-field in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] is called the simplified formulation, and the phase-field presented in this work, using Eq. ( 13), is called the complete formulation. Two particular situations should be considered for the gradient of the fracture parameters ∇ (G c l c ): (i) a multiphase material with homogeneity inside each phase, and (ii) a single-phase material with continuously varying property (for instance, a functionally-graded material). In situation (i), the gradient ∇ (G c l c ) is non-zero only at the interfaces between different phases; whereas it is non-zero everywhere in situation (ii). The analyses presented in this work focus on situation (i), which is more often encountered in engineering problems.

For an isotropic elastic body, the initial free energy density in Eqs. ( 12) and ( 13) can be expressed as:

ϕ 0 (ε) = 1 2 λtr 2 (ε) + µε : ε ( 15
)
where λ and µ represent the Lamé coefficients of the material. In order to activate the damage evolution only under tension, a split of the initial strain energy into tensile and compressive parts is necessary under complex loading conditions.

In the work of Miehe [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], a spectral decomposition of the strain tensor has been applied in order to fully distinguish the tensile and compressive strains:

ε = 3 n=1 ε n p n ⊗ p n = ε + + ε - (16)
where ε n and p n denote the eigenvalues and eigenvectors of the strain tensor, respectively.

The operator ⊗ denotes the outer product of two vectors. Note that ε + and ε -represent the tension and compression parts of the strain tensor, respectively:

ε + = 3 n=1 ε n + p n ⊗ p n , ε -= 3 n=1 ε n -p n ⊗ p n (17)
with the bracket operators x + = x+|x| 2 , and

x -= x-|x| 2 .
Miehe's spectral decomposition can be performed with the help of fourth-order projection tensors, P + and P -, which are shown as follows:

ε + = P + : ε, ε -= P -: ε, and P + + P -= I, (18)
where I denotes the fourth-order identity tensor. As detailed in Ref. [START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF] and Appendix A of Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], the calculation of the projection tensors is mainly based on the eigenvectors (w 1 , w 2 , w 3 ) and the eigenvalues (v 1 , v 2 , v 3 ) of ε:

P + = W T SW, W =      w 11 w 12 w 13 w 21 w 22 w 23 w 31 w 32 w 33      , S =      s 1 0 0 0 s 2 0 0 0 s 3      , s i =    1, v i > 0, 0, v i ≤ 0. (19)
The tensile and compressive part of the initial strain energy can then be expressed as:

ϕ ± 0 = 1 2 λ ( tr (ε) ± ) 2 + µε ± : ε ± . ( 20
)
As damage only affects the tension part, the stress tensor can be defined as:

σ = g (d) [λ tr (ε) + I + 2µε + ] + [λ tr (ε) -I + 2µε -] , (21)
with I the identity matrix.

Damage models should reflect the physics, where a fully damaged material cannot be damaged further, and the damage state cannot be recovered. Consequently, the boundedness d ∈ [0, 1] and irreversibility condition ḋ ≥ 0 should be met. For Miehe's model, the boundedness d ∈ [0, 1] is dealt with using a so-called history field (H) of maximum positive initial strain energy that replaces ϕ 0 in Eqs. ( 12) and ( 13). This history field is given as:

H (x, t n ) := max 0≤t≤tn ϕ + 0 (x, t) . ( 22
)
Thus, Eq. ( 12) should rewritten as follows:

g (d) H + G c l c d -∇ • (G c l c ∇d) = 0, (23)
as well as Eq. ( 13), which can be reexpressed as:

g (d) H + G c l c d -G c l c ∆d = 0. ( 24
)
Note that most studies using this model rely on Eq. ( 24) and a finite element solver. The weak form on which the finite element method is based implies an integration by parts which eliminates the divergence operator from the phase-field evolution equation. It does not make sense, therefore, to distinguish between a simplified and a complete formulation for finite element solvers. The following discussions are specific to FFT-based solvers. It is possible that they could be generalized to finite difference solvers as well but it is outside of the scope of our work.

Implementation of the complete formulation into an FFT-based solver

In the work of Chen et al. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], a fixed-point algorithm was adopted to solve the phasefield governing equation, and the mechanical part was solved separately. A weak coupling scheme was applied to bridge the phase-field and mechanical parts. This method is relatively straightforward to implement. Thus, the same algorithm is adopted in this work and is not detailed. It is reminded that periodic boundary conditions are intrinsic to FFT-based methods, and that those methods can only rely on regular (Cartesian) grids composed of hexahedral elements that we coin voxels. All fields (including the strain field, the stress field and the phase-field) are discretized as voxel-wise constant variables. More details can be found in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF].

Due to the complexity of the complete formulation, the solving procedure is slightly more intricate and is hence detailed in the following. The main ideas are inspired from the work of Jeulin [START_REF] Jeulin | Towards crack paths simulations in media with a random fracture energy[END_REF].

First, a term Q(x), Q(x) = G c (x)l c (x), (25)
is introduced to describe the multiplication of G c and l c . A fluctuation term Q(x) is then introduced, and Eq. ( 25) can be reexpressed as:

Q(x) = Q 0 + Q (x) , ( 26
)
where

Q 0 = max x∈Ω (Q (x)) + min x∈Ω (Q (x)) 2 . ( 27
)
Therefore, Eq. ( 23) can be rewritten as:

-2 (1 -d (x)) H (x) + G c (x) l c (x) d (x) -∇ • (Q 0 ∇d (x)) -∇ • Q (x) ∇d (x) = 0, ( 28
)
with

g(d) = (1 -d) 2 , g (d) = -2 (1 -d) .
In the weak coupling scheme, the damage variable at time t n+1 , denoted by d t n+1 , is calculated based on the history field at time t n , denoted by H tn . This H tn remains constant when solving the damage variable d t n+1 . Since Q 0 is constant over the domain, Eq. ( 28) gives:

G c (x) l c (x) + 2H tn (x) d t n+1 (x) -Q 0 ∆d t n+1 (x) = 2H tn (x) + ∇ • Q (x) ∇d t n+1 (x) . ( 29
)
Setting

A tn (x) = G c (x) l c (x) + 2H tn (x) , ( 30a
)
B tn (x) = 2H t n (x), ( 30b
)
D t n+1 (x) = ∇ • Q (x) ∇d t n+1 (x) , (30c)
Eq. ( 29) can be rewritten as:

A tn (x)d t n+1 (x) -Q 0 ∆d t n+1 (x) = B tn (x) + D t n+1 (x). ( 31
)
Polarizing A tn with

A tn 0 = max x∈Ω (A tn (x)) + min x∈Ω (A tn (x)) 2 ,
Eq. [START_REF] Kabel | Mixed boundary conditions for fft-based homogenization at finite strains[END_REF] gives

A tn 0 d t n+1 (x) -Q 0 ∆d t n+1 (x) = τ (x) , ( 32
)
with

τ (x) = B tn (x) -A tn (x) -A tn 0 d t n+1 (x) + D t n+1 (x) .
In Fourier space, the Laplacian operation in Eq. ( 32) can be easily transformed into a multiplication with frequency vector ξ, which gives:

F ∆d t n+1 (ξ) = -(ξ • ξ) d t n+1 ,
where F [-] denotes the Fourier transform operation, and [-] describes components in Fourier space. Hence, at a given fixed-point iteration i, a new damage field can be solved from Eq. ( 32):

d t n+1 i+1 = τ i A tn 0 + Q 0 (ξ • ξ) . ( 33
)
The fixed-point algorithm can be summarized as follows:

Initialization, if t n = 0 :

(a 0 ) d 0 (x) = 0 (b 0 ) Calculate Q (x) and Q 0 if t n = 0 : (a 1 ) A tn (x)
τ i = F[τ i ], i th iteration (b 2 ) d t n+1 i+1 (ξ) = τ i (ξ) A tn 0 +Q 0 .(ξ•ξ) (c 2 ) d t n+1 i+1 (x) = F -1 [ d t n+1 i+1 ] (x) (c 2 ) Calculate D t n+1 i+1 (x), and τ i+1 (x) (d 2 ) e = ||d t n+1 i+1 -d t n+1 i || 2 ≤ e crit , Convergence check (34)
The calculation of D t n+1 is performed in Fourier space because the derivation operation in real space is a simple multiplication in Fourier space:

D t n+1 = F -1 J • ξ • F Q • F -1 J • ξ • d t n+1 (35)
where J = √ -1. Compared to the simplified formulation, the complete formulation adds one more Fourier transform and two more inverse transforms.

Compared to Ref. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the viscous term, which was introduced for numerical stabilization for a rate-independent material, is not presented in this work. This is because from our experience, the viscous regularization has never really been needed in practice, as the solution of the phase-field problem usually converges well with Anderson's acceleration technique.

Concerning the convergence criterion, in this work, the one for phase-field problem is set as ||d i+1 -d i || 2 ≤ e crit = 1.00 × 10 -6 . The convergence criterion of the mechanical problem is set as

( ||∇.σ i+1 || 2 2 ) 1 2 || σ i+1 || 2 ≤ e crit meca = 1.00 × 10 -4
, where the brackets around a variable denote a volume average over the domain, for instance σ = 1

|Ω| Ω σdΩ. The convergence criteria and the tolerance values are the same as in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF].

To speed up the convergence of the fixed-point algorithm, Anderson's acceleration algorithm is incorporated in this work. This technique has been implemented in the software AMITEX [START_REF]AMITEX-FFTP[END_REF], which is also the main environment that is chosen for this work. Note that the grid is distributed in AMITEX, which means all operations are solved in parallel including the mechanical and phase-field solvers. The principal idea of Anderson's acceleration algorithm is to give an improved solution based on the four previous iterations.

A detailed description of this algorithm can be found in Ref. [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF]. Fig. 2 presents a detailed flowchart of the implementation of the complete formulation.

For comparison, the flowchart of the simplified formulation can be found in Appendix B of Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]. As illustrated in Fig. 2, the history field at time t n is applied to calculate the damage field at time t n+1 . Thus, a sufficiently fine time step is required. According to results from Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], a time step δ ε = 5.00 × 10 -7 is utilized in the present work. As for the FFT solver, the work in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] utilized a modified discrete Green operator, the so-called rotated-scheme proposed by Willot [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. In the present work, the basicscheme, the original method proposed by Moulinec and Suquet [START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF], is adopted for the mechanical and phase-field problems for both simplified and complete formulations. A detailed description and quantitative studies of these two schemes can be found in Ref. [START_REF] Ma | Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods[END_REF].

Numerical experiments and discussions 4.1 Unit-cell model with two half fibers on the borders

First, a unit-cell model with two half fibers on the borders is analyzed (as shown in Fig. 3). The model size is 2.00 × 10 -2 mm with a fiber volume fraction of V f = 0.55, and an inter-fiber distance l f = 3.29 × 10 -3 mm, while one voxel size is used in the thickness (z-direction). According to our previous work [START_REF] Ma | Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods[END_REF], for a unit-cell composite model, a resolution of N = 225 is sufficient to get a proper response. This is also used in the present work, which gives a voxel size of h = 8.89 × 10 -5 mm.

The blue part is the fiber for which the material is E-glass with properties [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]: In this series of simulations, the value of l c varies from 5.40 × 10 -4 mm to 3.30 × 10 -3 mm while the inter-fiber distance is fixed. The minimum value of l c satisfies Miehe's criterion, which states that l min c /h = 6.00 > 2.00, where h denotes the voxel size.

E f =
In the following, the model in Fig. 3 is subjected to a macroscopic strain along the xdirection (transverse direction) with a time step δ ε xx = 5.00×10 -7 until the final failure while stress-free conditions are imposed in the other directions ( σ yy = σ zz = σ xy = σ xz = σ yz = 0). These mixed-type loading conditions are applied by using the method presented in [START_REF] Kabel | Mixed boundary conditions for fft-based homogenization at finite strains[END_REF]. Due to stress-free boundaries, Poisson's effect is active during loading.

Simplified formulation

Considering the fiber is much stiffer than the matrix, and the G c value of fiber is much larger than that of the matrix, a logical hypothesis is that the crack would occur in the matrix, and the damage state in the fiber would be close to null in the present case.

However, Then local damage profiles along the center line for different l c values are plotted in Fig.

5.

As shown in the figure, whatever the l c values, the damage profiles are so smooth that we can not even distinguish the transition between matrix and fibers. Thus, another hypothesis can be proposed: the damage state in fibers is not physical damage. They are mainly due to the diffusion of matrix damage into the fibers. To verify the proposed hypothesis, higher critical energy release rate values G f c are tested for the fiber to check whether they can reduce the damage in fibers. The present work takes two groups of l c values for tests (l c = 1.30 × 10 -3 mm and 3.30 × 10 -3 mm). For each group, G f c values vary from 6.00 × 10 1 N mm -1 to 9.00 × 10 20 N mm -1 . Such unrealistically high values are used to justify the proposed hypothesis. Results are plotted in Fig. 6. When G f c increases to a very high value, we would logically expect the damage variable in the fibers to decrease. However, as seen in the figure, whatever the G f c values, the macroscopic stress-strain curves do not change for the two test groups, and the damage profiles remain stable. This phenomenon confirms that the damage in the fiber is artificially caused by the diffusion of matrix damage to the fiber side.

Complete formulation

In this part, the same geometric model, materials and settings are adopted for testing the complete formulation. The results are shown in Figs. 7, 8 and 9. Fig. 7 indicates that the maximum damage state in the fiber is negligible after applying the complete formulation.

In Fig. 8, the boundaries between fibers and matrix can be easily identified for all l c The differences between the simplified and complete formulations on macro-level stressstrain relations are shown in Fig. 10. This figure shows that the complete formulation shows higher peak stress and final failure strain than the simplified formulation. The effect of the damage diffusion on the macroscopic stress-strain curves is not obvious for a smaller l c value (5.40 × 10 -4 mm), while it becomes more and more significant as l c increases. However, a conclusion can be obtained: for heterogeneous models, especially when a large l c is being used, the additional term in the complete formulation has a significant influence on both local crack patterns and macro-level mechanical behaviors.

Study of iteration numbers for convergence

The above results show that the complete formulation eliminates unphysical damage diffusion between phases. However, adding the additional term (shown in Eq. ( 14))

gives a high cost. Fig. 11 provides the total fixed-point iterations of the simplified and complete formulations. As can be seen, the complete formulation needs much larger iteration numbers to get convergence than the simplified formulation. This difference is more than 10 2 times.

Single notch model

To better understand the impact of Anderson's acceleration scheme and the contrast of damage properties, a simpler model (single-notch model) is studied in this work. As shown in Fig. 12, it is a square plate of a length of 1 mm. A straight horizontal notch of 5.00 × 10 -1 mm is located at the mid-height of the left edge. The width of the notch is one voxel size. Because the FFT-based solver is intrinsically periodic, we need to add lateral bands, denoted as M2 (purple material), to break the periodicity and mimic the stress-free boundary conditions. Material M2 is also used for the notch. Note that, in this study, the resolution is given as N = 251 (mesh size h = 4.00 × 10 -3 mm) and 2).

For all three tests, a tensile loading along y-direction with time step δ ε yy = 5.00 × 10 -7

σ xx = σ zz = σ xy = σ xz = σ yz = 0).
The solutions computed for all values of G c and for the two formulations were the same and correspond to the result in Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], with the crack propagating along the same line as the pre-crack. The local damage fields are not reported herein as we are only interested in the influence of the contrast on the iteration numbers.

Fig. 13 shows the total fixed-point iteration numbers with different contrasts. Compared to the simplified formulation, the complete formulation is sensitive to the contrast of damage properties. The higher the contrast is, the more iterations are required. Moreover, the slope of the complete formulation is almost linear. Furthermore, Fig. 13 also indicates that Anderson's acceleration scheme can significantly reduce the iteration numbers. However, for the complete formulation, it is much less effective for a high contrast model. Hence, a more efficient acceleration algorithm needs to be proposed in the future when dealing with high contrast mediums.

Single fiber unit-cell with notch defect under tension

In this third example, a single fiber unit-cell model with a notch defect is studied.

As shown in Fig. 14, a notch of 6.00 × 10 -3 mm is placed at the mid-height of the left edge with one voxel width. The model is a square plate of 2.00 × 10 -2 mm with a thickness of one voxel in the third direction. The circular reinforced fiber is at the mid-height with a diameter of 6.00 × 10 -3 mm and 1.30 × 10 -2 mm away from the left edge. Because the FFT solver is intrinsically periodic, we need to add lateral bands. Each lateral band has a width of five voxels. Meanwhile, the material of the lateral bands is also used for Regarding the material properties, the fiber (blue part) is E-glass with properties: E f = 7.40 × 10 1 GPa, ν f = 0.2, and G f c = 6.00 × 10 1 N mm -1 . These values are chosen to prevent failure of the fibers. The gray region is the matrix for which the material is Epoxy with properties: E m = 4.65 GPa, ν m = 0.35, G m c = 9.60 × 10 -4 N mm -1 . These two materials are the same as in Sec. 4.1. For the lateral bands and notch, where there should be a void, one ten-thousandth of matrix elastic properties are taken to avoid the convergence issue. They are: E e 0 = 4.65 × 10 -4 GPa, ν e = 0.35, and G e c = 9.60 × 10 -4 N mm -1 , the same as the matrix.

The l c values are taken from l c = 2.70 × 10 -4 mm, which satisfies Miehe's criterion l min c /h ≈ 3.00 > 2.00. For all tests, a tensile loading along y-direction with time step δ ε yy = 5.00 × 10 -7 is applied while stress-free conditions are imposed in other directions

( σ xx = σ zz = σ xy = σ xz = σ yz = 0).
The results of local fields for different l c values are shown in Fig. 15. As seen in the figure, both formulations can capture the complex crack propagation pattern. A similar phenomenon in Sec. 4.1 can also be observed herein: the frontiers between fiber and matrix can be observed in the complete formulation, while this transition is smeared out in the simplified formulation. The larger the l c is, the more critical the damage diffusion is. This phenomenon can be shown in profile plots as shown in Fig. 16. In contrast, the complete formulation does not have this issue. In this part of the work, our scope is to analyze the damage state in the matrix and fiber. Thus, the damage in the notch is not our focus point and is not plotted in Fig. 16.

We also observe in Fig. 15(c,f) the diffusion of damage within each phase, as there are non-zero damage values away from the crack in the matrix phase. This is a well-known property of Miehe's phase-field model that has been eliminated in alternative models [START_REF] Wu | Phase-field modeling of fracture[END_REF].

Its analysis is out of the scope of the present work as it is a modeling issue that is not specific to the FFT-based method and heterogeneous materials. Crack propagation patterns predicted using the complete formulation (Fig. 15(a,b,c)) seem challenging to explain from a physical standpoint. After studying the whole damage process, we observe that the two formulations actually share the same first stage. For instance in the simulations with l c = 2.70 × 10 -4 mm, the crack initiates from the notch and propagates until the fiber where it is arrested in the first stage. Then, in the second stage, which is also similar for both formulations, the damage in the matrix continues increasing in the red zones shown in Fig. 15(a,d) rather than at the crack initiated in the first stage. Finally, one or two new cracks initiate in these red zones and propagate until joining the initial crack. Animated pictures showing the second and third stages can be found in the supplementary material.

The difference between the two formulations is that the complete formulation shows later failure, which makes the initiation of cracks in both red zones possible. It is possible that this phenomenon that may seem unphysical is due to the absence of a debonding model at the interface. In this work, indeed, the fracture properties of the interface are simply that of the matrix. This issue might be dealt with by introducing a debonding model in the complete formulation. This will be considered in future work. Fig. 17 shows the macro-level σ -ε relations with different formulations. This figure reemphasizes that the term ∇ (G c l c )•∇d produced non-negligible effects in the high contrast heterogeneous models. The complete formulation offers a higher final failure strain than the simplified formulation. However, the two formulations have similar responses for the low-level damage state, which can also be observed in Sec. 4.1. That is mainly because the ∇d value is small for the low-level damage state, which reduces the effect of the term

∇ (G c l c ) • ∇d.
The tests above (Fig. 15) have shown that the complete formulation can avoid the unrealistic fiber damage for a high contrast model (G f c G m c ). We conduct another series of tests to verify whether the complete formulation can break two phases for low contrast (G f c ≈ G m c ). Hence, in these tests, the matrix properties are unchanged while the fiber properties are reduced to twice those of the matrix (E f = 2E m and G f = 2G m ). Fig. 18 shows that the complete formulation can break two phases as the simplified formulation does. The results of the two formulations are also similar from the point of In conclusion, the lower the contrast is, the more similar the results are.

Single fiber unit-cell with notch defect under shear

To compare the two formulations under different loading conditions, we simulate a model similar to that used in Sec. 4.3 under in-plane shear loading. We delete the two lateral bands, which results in the periodic model shown in Fig. 20.

The materials and properties are the same as in Sec. 4.3. The loading conditions are similar, except that the strain is imposed in direction xy in this section, and the stresses in the other directions are free.

The results are shown in Fig. 21. Both the simplified and complete formulations predict well the in-plane shear crack patterns. Note that the cracks in the red zones are mainly due to the periodic boundary conditions. In these results, because the cracked zones are far from the interface between the fiber and the matrix, it can be seen that the predictions of these two formulations show less difference in terms of local-level crack patterns and macro-level mechanical responses. We can conclude that the difference between the two formulations should be less pronounced when the crack initiates and propagates away from the interfaces.

Multi-fiber model

According to the above tests, we have shown that the simplified formulation can lead to a non-physical damage diffusion for a multi-phase heterogeneous model, e.g., composites, while the complete formulation avoids this. The above tests only show simple examples.

Thus, a more complex model corresponding to a unidirectional fiber-reinforced composite material is tested in this part. This model shown in Fig. 22 is reproduced from Ref. [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF] as well as the material properties: the blue part is fiber, of which the material is E-glass with properties: E f = 7.40 × 10 1 GPa and ν f = 0.20. The value of G f C is not given in Ref. [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF] hence we assume that G f c = 6.00 × 10 1 N mm -1 based on Ref. [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]. The red part is the matrix, of which the material is Epoxy with properties: E m = 4.65 GPa, ν m = 0.35, G m c = 9.60 × 10 -4 N mm -1 . The result for the simplified formulation is shown in Fig. 23(a,b). We observe once again a diffusion of damage from the matrix to the fibers. This re-confirms the damage diffusion between phases for the simplified formulation. As demonstrated previously for a simple model, the result obtained using an unrealistically large value for G f c in Fig. 23(c,d) shows that the micro-level damage fields, especially the fiber damage, do not vary. Fig. 23(e) demonstrates that the macro-level mechanical response does not vary either when increasing G f c to unrealistically large values. The comparison between the two formulations is shown in Fig. 24. It clearly shows that the complete and simplified formulations can give totally different results for microlevel damage fields and macro-level mechanical responses as well. This is because the l c value is large compared to the inter-fiber distance l f in the model. This is challenging because we can neither reduce l c due to Miehe's criterion (l c /h > 2) and its implication on computational cost nor arbitrarily increase l f because it depends on the material.

Meanwhile, it should be admitted that we cannot judge which results are correct without experimental results. Unfortunately, Ref. [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF] also lacks experimental results and uses another type of phase-field that is not comparable. Thus, more studies should be carried out in the future to validate the complete formulation's reliability. This multi-fiber model can be a basis for a thorough validation of the model in a future work. Micro-configuration of the fiber debonding cracks can be observed in-situ, as done for example in Refs. [START_REF] Hobbiebrunken | Evaluation of interfacial strength in CF/epoxies using FEM and in-situ experiments[END_REF][START_REF] Sharma | On the prediction of the bi-axial failure envelope of a UD CFRP composite lamina using computational micromechanics: Effect of microscale parameters on macroscale stress-strain behavior[END_REF]. The important features of this phenomenon are an angle between the debonded surface of fibres and the sequence of linking of the debonding cracks into a continuous transverse crack. The ability of a model to represent such features for a microstructural configuration of fibres taken from in-situ micrographs should be a validation criterion.

Conclusions

In this paper, two formulations of a phase-field model are addressed and compared: one is called the simplified formulation, while another is called the complete formulation. The main difference between these two formulations is that the term ∇(G c l c ) • ∇d is omitted in the simplified formulation, where details can be found in Eq. [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF]. For a homogeneous model, these two formulations are equivalent. In contrast, for a heterogeneous model, the results of the two formulations diverge because the G c values are not constant anymore.

The higher the contrast is, the more critical the divergence is.

The simplified formulation can converge much faster than the complete formulation. For the same problem, the complete formulation requires more iteration numbers for convergence. The higher the contrast is, the more obvious the difference is. One critical aspect, however, can be observed when using the simplified formulation for a heterogeneous model: non-physical damage diffusion occurs between phases. This diffusion is eliminated with the complete formulation. It highly impacts the evaluation of the local micro-level damage state. The difference in macro-level σ -ε relations between the two formulations cannot be neglected either. The larger the l c value is and the closer to the interfaces the cracks are, the more obvious these two phenomena are.

Moreover, for Miehe's phase-field model, the l c value plays a vital role in the mechanical response, which is often considered a material parameter in the literature. This parameter cannot, therefore, be arbitrarily reduced in order to diminish the difference between the simplified and complete formulations. Quantitative comparisons with experimental tests should hence be carried out in future work to find out which formulation is more reliable.

In the present work, Anderson's acceleration algorithm has been applied to speed up the convergence. Despite this, when applying the complete formulation, Fig. 13 shows that the increase in iteration numbers compared to the simplified formulation cannot be neglected. Besides, Anderson's acceleration scheme becomes less efficient for large contrast models. Thus, proposing a sufficiently efficient acceleration algorithm is one of the important perspectives for the next step. Replacing the basic-scheme and its truncated Green operator by more robust alternatives might also alleviate the influence of the contrast on the iteration numbers.

Meanwhile, the authors also noticed that a multi-phase-field technique for heterogeneous mediums has been proposed to solve the above issues. One of the main drawbacks of the multi-phase-field approach is to create one phase-field for each phase. That will increase the complexity when dealing with models with a large number of phases. In contrast, the complete formulation only needs one formulation to describe all phases. Thus, it will also be interesting to compare the efficiency of these two approaches in the future.

Figure 1 :

 1 Figure 1: A solid body Ω with the crack set Γ (or Γ l ): (a) sharp cracks and (b) approximated diffuse crack bands. (c) Illustration of diffusive damage profile of a crack at x = 0 for various length parameters l c .

Figure 2 :

 2 Figure 2: The flow chart of the implementation of the complete formulation into an FFT solver. CV means convergence check.

Figure 3 :

 3 Figure 3: The composite model with two half fibers on the borders (blue part: fiber; red part: matrix).

Fig. 4

 4 Fig. 4 shows the local damage fields. When l c = 5.40 × 10 -4 mm, there is an apparent distinction between matrix and fibers. In contrast, this transition becomes smeared for l c = 3.30 × 10 -3 mm.

Figure 4 :

 4 Figure 4: The local damage fields at the final failure state for different l c values and for the two half fibers model: (a) 5.40 × 10 -4 mm; (b) 8.10 × 10 -4 mm; (c) 1.30 × 10 -3 mm; (d) 2.20 × 10 -3 mm; (e) 3.30 × 10 -3 mm. The orange dashed lines are the fiber contours.

Figure 5 :

 5 Figure 5: The damage profiles along the center line for different l c values and for the two half fibers model. The purple dashed lines are boundaries between fibers and matrix.

Figure 6 :

 6 Figure 6: The local damage profiles and macro-level stress-strain relations for different G f c values and for the two half fibers model: (a), (b): results for l c = 1.30 × 10 -3 mm; (c), (d): results for l c = 3.30 × 10 -3 mm. The purple dashed lines are boundaries between fibers and matrix, and different curves are superimposed.

Figure 7 :

 7 Figure 7: The maximum damage d max f

Figure 8 :Fig. 9 .

 89 Figure 8: The local damage fields at the final failure state for different l c values and for the two half fibers model by applying the complete formulation: (a) 5.40 × 10 -4 mm; (b) 8.10 × 10 -4 mm; (c) 1.30 × 10 -3 mm; (d) 2.20 × 10 -3 mm; (e) 3.30 × 10 -3 mm. The orange dashed lines are the fiber contours.

Figure 9 :

 9 Figure 9: The damage profiles along the center line for different l c values and for the two half fibers model by applying the complete formulation. The purple dashed lines are boundaries between fibers and matrix.

Figure 10 :

 10 Figure 10: The macro-level σ -ε relations for different l c values and for the two half fibers model by applying different formulations: (a) 5.40 × 10 -4 mm; (b) 8.10 × 10 -4 mm; (c) 1.30 × 10 -3 mm; (d) 2.20 × 10 -3 mm; (e) 3.30 × 10 -3 mm.

Figure 11 :

 11 Figure 11: Total fixed-point iteration numbers from the beginning of the loading to the final failure state for the two half fibers model.

Figure 12 :

 12 Figure 12: The geometry of the single notch plate. Material 1 (M1) and material 2 (M2) are the yellow and purple parts, respectively.

Figure 13 :

 13 Figure 13: Total numbers of fixed-point iterations for the phase-field solution from the beginning of the loading to the final failure state for the single notch model.

Figure 14 :

 14 Figure 14: The geometry of composite model with defect (blue: fiber; grey: matrix; red: notch).

Figure 15 :

 15 Figure 15: The local damage fields for different l c values and for the model with notch and fiber under tension by applying the complete (top) and simplified formulations (bottom): (a,d) l c = 2.70 × 10 -4 mm; (b,e) l c = 5.40 × 10 -4 mm; (c,f) l c = 8.10 × 10 -4 mm. The orange dashed lines are the fiber contours, and the red zone is explained in the text.

Figure 16 :

 16 Figure 16: The damage profiles along the center-line for various l c values and for the model with notch and fiber under tension by applying the complete and simplified formulations: (a) l c = 2.70 × 10 -4 mm; (b) l c = 5.40 × 10 -4 mm; (c) l c = 8.10 × 10 -4 mm. The purple dashed lines are boundaries between fibers and matrix.

Figure 17 :

 17 Figure 17: The macro-level σ -ε relations for different l c values and for the model with notch and fiber under tension: (a) l c = 2.70 × 10 -4 mm; (b) l c = 5.40 × 10 -4 mm; (c) l c = 8.10 × 10 -4 mm.

Figure 18 :

 18 Figure 18: The local damage fields for different l c values and for the model with notch and fiber under tension by applying the complete (top) and simplified formulations (bottom): (a,d) l c = 2.70 × 10 -4 mm; (b,e) l c = 5.40 × 10 -4 mm; (c,f) l c = 8.10 × 10 -4 mm. The orange dashed lines are the fiber contours.

Figure 19 :

 19 Figure 19: The macro-level σ -ε relations for different l c values and for the model with notch and fiber under tension: (a) l c = 2.70 × 10 -4 mm; (b) l c = 5.40 × 10 -4 mm; (c) l c = 8.10 × 10 -4 mm.

Figure 20 :

 20 Figure 20: (a) The geometry for the single fiber unit-cell with periodic notch defect test case. (b) Infinite periodic medium from which this unit-cell is extracted (orange square).

Figure 21 :

 21 Figure 21: The local damage fields at final failure for l c = 2.70 × 10 -4 mm and for the model with notch and fiber under shear: (a) the simplified formulation, (b) the complete formulation. (c) The macro-level σ -ε relations. The orange dashed lines are the fiber contours.

Figure 22 :

 22 Figure 22: The complex model with multiple fibers inspired from Ref. [14].

Figure 23 :

 23 Figure 23: The local damage fields at final failure for the complex model using the simplified formulation with l c = 2.90 × 10 -4 mm: (a) G f c = 6.00 × 10 1 N mm -1 , (c) G f c = 9.00 × 10 10 N mm -1 . The fiber damage fields at the same state: (b) G f c = 6.00 × 10 1 N mm -1 , (d) G f c = 9.00 × 10 10 N mm -1 . (e) The macro-level σ -ε relations. The orange dashed lines are the fiber contours.

Figure 24 :

 24 Figure 24: The local damage fields at final failure for the complex model with l c = 2.90 × 10 -4 mm: (a) the complete formulation, (b) the simplified formulation. (c) The macro-level σ -ε relations. The orange dashed lines are the fiber contours.

Table 1 :

 1 's model, which lacks an elastic regime. We demonstrate in the following that this should not be the only reason. The maximum damage d max

	l c (×10 -4 ) mm 5.40 8.10 13.00 22.00 33.00
	d max f	0.15 0.24 0.38	0.56	0.68

Table

.

1 implies a converse condition: the maximum damage in the fiber is not null for all tests and rises with increasing the l c . One may argue that this is the intrinsic feature of Miehef in the fiber at the final failure state for various l c values and for the two half fibers model.

Table 2 :

 2 The damage properties for the three tests on the single notch model. is applied while stress-free conditions are imposed in other directions (

		Test 1 Test 2 Test 3
	G e c (N mm -1 ) 2.70	27	2700
	Contrast	1	10	1000
	G e c /G c


domain: phys spi

text:
Nicolas Vauchel
email: nicolas.vauchel@onera.fr

Eric Garnier

Thomas Gomez

Estimation of Sobol indices

Keywords: Uncertainty Quantification, Sobol indices, Multi-element polynomial chaos, Flight Dynamics, Machine learning

come

Introduction

Global Sensitivity Analysis (GSA) aims at finding the impact of uncertain parameters on a Quantity of Interest (QOI). In the present study, the chosen parameters are called the input parameters, are supposed to be random, independent and to follow a stochastic distribution of known probability density function (pdf). Sobol indices, relying on ANOVA decomposition [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], rank the actions of the input parameters and the interactions between them according to their impacts on the QOI. The indices also give a quantitative number reflecting these impacts. Sobol indices can be estimated with Monte-Carlo methods, but it usually requires a large number of simulations. Approximating the relation between the QOI and the input parameters using a surrogate model in the form of Polynomial Chaos (PC) model permits to establish an analytical link between the chaos coefficients and the Sobol indices [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. Once the model is learned, all the Sobol indices are directly available. However, PC models are not accurate in presence of irregularities or discontinuities. It may need a high truncation degree to achieve accuracy, which becomes problematic with the famous curse of dimensionality. Moreover, Gibbs phenomenon amplifies by increasing the degree of the model in presence of discontinuities.

Multi-elements (ME) models [START_REF] Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF][START_REF] Sargsyan | Uncertainty quantification given discontinuous model response and a limited number of model runs[END_REF][START_REF] Halder | An adaptive minimum spanning tree multielement method for uncertainty quantification of smooth and discontinuous responses[END_REF] were developed to address this issue. The abstract parameter space, containing all possible combinations the input parameters, is partitioned into subsets, called elements, where the QOI is locally smoother. Local PC models are trained on every element. The global model is then a piecewise PC model and contains a method to determine which local model to use when an evaluation is needed. However, with this piecewise model, the classical analytical formula with Sobol indices and PC models does not hold anymore, as the formula relies on orthogonal properties on all the parameter space. In an upcoming paper [START_REF] Vauchel | Estimation of Sobol indices using a multi-element Polynomial Chaos model[END_REF], the analytical link between PC models and Sobol indices is extended in the case of ME models with non time-dependent QOI's.

In the present paper, the approach is extended to the case of a ME model approximating a time-dependent QOI, solution of an Ordinary Differential Equation (ODE). The time-dependent version of a new ME method, called ME-ACD, developed in another upcoming paper [START_REF] Vauchel | A multi-element non-intrusive Polynomial Chaos method using agglomerative clustering based on the derivatives to study irregular and discontinuous Quantities of Interest[END_REF], is used, although any ME polynomial model able to handle time-dependent QOI's would have been appropriate. At every considered instants, the values of all the Sobol indices are obtained with the analytical formula linking them to the coefficients of the local models. The new approach is then able to give an estimation of the temporal evolution of Sobol indices with a relatively low number of runs, even if at a given instant, the QOI is irregular or discontinuous according to the input parameters (traducing the presence of bifurcations).

The method is used in practical applications of Flight Dynamics. Flight Dynamics aims at studying the motion of an aircraft modelled as a solid. The forces and the moments of the fluid on the aircraft rely on a model with a chosen structure and come from data, which usually are data from wind-tunnel experiments. These models are therefore uncertain.

Traditionally, the bifurcation theory is used on the dynamical system and equilibria are linked to behaviours of the aircraft. At low angle of attack, equilibria are usually safe trajectories but in the post-stall domain, equilibria are usually dangerous behaviours, as spins. It has been shown in [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF] that the uncertainties of the model of the aerodynamic forces and moments can have some influence on some predicted spins, notably on their stability. For these behaviours, some combinations of the uncertain parameters can lead to instability and others to stability, which implies that, at a given instant, state variables can present irregularities and/or discontinuities with respect to these input parameters. The method presented in this paper is then suitable to conduct a joint Uncertainty Quantification (UQ) and GSA study on the temporal evolution of the state variables of this system.

Notations

The QOI is denoted ŷ. ŷ depends on d ∈ N * input parameters p = (p 1 , ..., p d ). Every p m follows a stochastic distribution law of known pdf ρ m and all the input parameters are supposed independent. To stay close to the original definition of Sobol indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], every p m is mapped to a corresponding x m following an uniform law in the interval [0, 1]. The mapping takes into account the distribution law of the input parameter

x m = C m (p m ) (1)
where C m is the cumulative distribution function of p m . If ρ m > 0 (and not only ρ m ≥ 0), the mapping is bijective. This hypothesis is supposed verified in the scope of this study. From now on, "input parameters" denotes x = (x 1 , ..., x d ) and "QOI" denotes y(x) = ŷ(p).

GSA aims at ranking the impact of the actions and of the interactions between the input parameters on the QOI. A multidimensional indexing is used to denote these actions/interactions. The subscript i ∈ R j , j ∈ R denotes all items referring to the interactions between the j input parameters x i 1 , ..., x i j . I = {0} contains the index denoting the part of the QOI which does not depend on any input parameter.

For a given vector x ∈ R d , and a given i ∈ I (d) j , x i = (x i 1 , ..., x i j ) is the vector composed of the components of x indexed by the components of i and x ∼i is the vector composed of the d-j other components.

To denote the terms of multivariate polynomial series depending on d variables, another multidimensional indexing logic is used. k = (k 1 , ..., k d ) denotes the term being a coefficient times a product of d univariate polynomials of degree k m depending of x m . Series of polynomials of degree γ ∈ R are then series containing terms denoted by an index k with d m=1 k m ≤ γ. The set containing all these possible indices is denoted K

Sobol indices

In [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], the so-called ANOVA decomposition is introduced

y(x) = d j=0 i∈I (d) j y (j) i (x i ) (2)
with

y (0) 0 = H d y(x)dx (3)
and, ∀j ∈ [|1, d|] ,∀i ∈ I

(d) j , ∀m ∈ {i 1 , ..., i d } 1 0 y (j) i (x i )dx m = 0 (4)
The variance of the QOI is denoted D. With ANOVA decomposition (2), D can be decomposed as a sum of terms called "partial variances"

D = d j=1 i∈I (d) l D (j) i (5)
Every D

(j)

i represents the impact of interaction i of order j on the total variance D. Sobol indices are then defined as

S (j) i = D (j) i D (6)
Partial variances can be estimated with Monte-Carlo methods [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], using random or pseudo-random samplings as Sobol sequences. To estimate all the first order partial variances, (d + 1)N computations of the QOI are needed. To estimate all the second order partial variances, d 2 N additional computations of the QOI are needed and so on. To estimate all the partial variances, 2 d N computations are needed. In [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], a method to estimate the total Sobol indices (being the sum of all the sobol indices involving a given parameter appears) is also available.

Multi-element Polynomials Chaos models

First introduced by Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] to approximate random variables as series of Hermite polynomials depending on Gaussian variables, PC models were later extended in [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF] as series of polynomials depending on random variables having another distribution and having orthogonal properties corresponding to this distribution. PC models' coefficients can be obtained with Galerkin methods or Non Intrusive (NI) methods [START_REF] Xiu | Fast numerical methods for stochastic computations : a review[END_REF]. In the scope of this article, only NI methods are considered. With NI methods, a given number of deterministic computations of the QOI are realised with different combinations of input parameters. These computations are used to get the chaos coefficients. In this paper, the computations are called "training observations". An observation can be seen as a point in an abstract space. The abstract space where the observations are seen as points of coordinates the values of the input parameters is called the parameter space H d . It is the unit hypercube of dimension d.

With the mapped input parameters, PC expansion is

y(x) = k∈K (d) γ α k Φ k (x) (7)
H d Φ a (x)Φ b (x)dx = d m=1 δ am,bm (8)
where (Φ k ) is the normalised shifted multivariate Legendre polynomials basis, products of univariate polynomials orthonormal with respect to the scalar product associated with the uniform law in [0, 1]. δ m,n is the standard Kroenecker symbol.

An analytical formula [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] links Sobol indices of a QOI with the coefficients of its PC model. Conversely to the Monte-Carlo method to get the Sobol indices, all the training observations are used to get the surrogate model and then all Sobol indices are available. A better exploration of the parameter space is therefore possible. However, PC models suffer from the curse of dimensionality, the number of terms of the chaos series being equal to γ+d γ . More than else, accuracy of PC models deteriorates in presence of irregularities and/or discontinuities. ME methods permit to get piecewise polynomial models. The parameter space H d , is partitioned into N E subsets E e called elements where the QOI is locally smoother. ME-ACD, developed in an submitted paper [START_REF] Vauchel | A multi-element non-intrusive Polynomial Chaos method using agglomerative clustering based on the derivatives to study irregular and discontinuous Quantities of Interest[END_REF], permits to get ME polynomials models expressed as

y(x) = N E e k∈K (d) γe c e,k Ψ k (x) I Ee (x) (9)
where Ψ k (x) = d m=1 x km m are the polynomials of the canonical basis of multivariate polynomials depending on d variables and I Ee is the indicator function returning 1 if x ∈ E e or 0 otherwise. With a non time-dependent QOI, ME-ACD algorithm uses agglomerative clustering on the training observations represented as points in a abstract space taking into account the values of the QOI and the values of its derivatives with respect to the input parameters. It permits to partition the parameter space in presence of discontinuities and irregularities with respect to the input parameters. The algorithm tests partitions having between N E min and N Emax elements, and for all these partitions, local polynomial models are trained using Ordinary Least Square method with the training observations located in the elements. Local models of degree between n min and n max are tested, and the optimal degree is selected as the one minimising the predicted residual error sum of squares (PRESS) error. The partition is then selected as the one minimising the maximum of PRESS errors of the local models with optimal degree. This maximum is denoted below "partition error" δ z . If, during the search, δ z ≤ δ min , the partition is kept and there is no need to look for another partition with more elements. If during the search, an element has a number of training observations lower or equal than N outlier , the observations are considered as outliers and are removed from the training set. The training is then conducted once again from the beginning. Computations of the derivatives are done with the use of OLS and first order Taylor theory on every training observations with its nearest neighbours in the space taking into account the values of the input parameters and the value of the QOI. Regressions with a number of neighbours from k min to k max are tested and the one minimising the PRESS error is kept. During the agglomerative clustering, the training observations are connected with their k connect closest neighbours in the parameter space. Once the local models are obtained, a neural network classifier is trained to know which local model to use when an evaluation is needed. In the present paper, a brute force search on a set of learning rates is realised to find the one leading the best accuracy of the classifier on the training set itself.

In ME-ACD expansion [START_REF] Wiener | The homogeneous chaos[END_REF], the canonical basis is used instead of a polynomial basis having orthogonal properties corresponding to the pdf of the input parameters like in classical PC models. The orthogonal properties need information on all the parameter space to be used, and not only on local elements. This explains why the formula linking standard PC models' coefficients to Sobol indices [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] cannot be used in the case of a ME model. In a second upcoming paper [START_REF] Vauchel | Estimation of Sobol indices using a multi-element Polynomial Chaos model[END_REF], an analytical formula linking the coefficients of the local polynomial models of a ME expansion to partial variances is found. In order to establish this link, the elements are first approximated as an union of hyperrectangles. To find the set of hyperrectangles, the parameter space (the d-dimensional unit hypercube) is split in a grid of N h hypercubes in every dimensions (thus being a grid of (N h ) d hypercubes in total). Every hypercube is then labelled with the index of the element it is located in with the use of the classifier of ME method on N random points sampled in the hypercube. Once the hypercubes are labelled, a merging operation with identical hypercolumns is effectuated for every direction. This leads to N r ≤ (N h ) d hyperrectangles. For every hyperrectangle, a change of basis is realised for the local model located in this hyperrectangle from the canonical basis to a mutivariate polynomial basis being orthonormal with respect to a scalar product proper to the bounds of the hyperrectangle. Therefore, the ME model is now a piecewise polynomial model over the hyperrectangles with polynomials having orthonormal properties proper to their corresponding local hyperrectangle. The coefficients in these new basis are then linked with partial variances with thanks to these orthonormal properties. This formula is therefore an extension of the analytical formula linking chaos coefficients of standard PC models to Sobol indices in the case of a ME model.

In the present paper, the QOI is supposed to be the solution of an ODE. It is then time-dependent. The ME model ( 9) is then extended

y(x, t) = N E (t) e k∈K (d) γe(t) c e,k (t)Ψ k (x) I Ee(t) (x) (10)
A discrete set of N t instants {t 1 , ..., t Nt } (indexed by increasing order) is considered. Taking advantage of the potential similarities of the QOI from one instant to the next one, for every instant t z with z ≥ 2, the partition obtained at the previous instant t z-1 is kept and the local models corresponding to its elements are trained with the information of the QOI at instant z. The partition error δ z is calculated. This error is then a metric of accuracy of the global piecewise model at instant t z with the partition corresponding to the previous instant. If δ z ≤ δ max , δ max being a threshold to choose, the partition and the local models are kept for this instant. In this case, the partition is the same for t z-1 and t z , but the coefficients of the local models differ. For the first instant or if δ z > δ max for z > 1, a new suitable partition is obtained with the clustering method of ME-ACD algorithm and local models are found at this instant with the values of the QOI at instant t z . A classifier is then trained for the new partition. Time-dependent ME-ACD leads to N p ≤ N t partitions with their corresponding classifiers. A list T of size N t links the instants to the corresponding partitions. At every instant, a set of coefficients for every elements of the partition considered for this instant is available.

For a dynamical system composed of several coupled ODE's, this algorithm can be repeated for every of the state variables of the system, all state variables being treated independently. As variables can have different scales, a normalised partition error δ z can be used to be compared with the chosen δ min and δ max , these two parameters being chosen for all the variables. In the present study, δ z is normalised by the maximum of the absolute value the considered variable can take in all the scenarios corresponding to training observations, at the considered instant t z .

Once every state variable has its temporal ME-ACD model, temporal evolution of all Sobol indices of all the state variables can be estimated with the method described above.

Applications in Flight Dynamics

Presentation of the system

In this section, the described approach is used to study the influence of some parameters of a Flight Dynamics system. The forces and moments of the fluid on the aircraft are traditionally built on data coming from wind tunnel experiments. In this study, the mathematical models of the forces and moments are built on data coming from experiments conducted at the vertical wind tunnel of ONERA Lille on a light aircraft scale model.

The motion of the aircraft in the terrestrial reference frame, supposed as an inertial reference frame, is modelled by a dynamical system of 8 coupled ODE's. The 8 state variables are :

-V : Norm of the velocity of the aircraft -α : Angle of attack -β : Angle of sideslip -p : roll rate -q : pitch rate -r : yaw rate -Θ : pitch angle -Φ : yaw angle and the parameters are :

-m : mass of the aircraft -ρ : density of surrounding air -S : reference surface of a wing -c : reference chord of a wing -A, B, C : Inertial moments of the aircraft -δ l , δ m , δ n : angular deflections of control surfaces, respectively aileron deflection, elevator deflection and rudder deflection.

Forces applying on the aircraft are -Its weight. Its mass is supposed constant in the scope of the study.

-Propulsion. In the scope of the study, propulsion is not considered.

-Aerodynamic forces from the fluid on the aircraft. In the scope of this study, the density of air is supposed constant.

The global aerodynamic coefficients are the dimensionless aerodynamic forces and moments. In this paper, the model has the following chosen structure

C x = C x0 (α) (11)
C y = C y0 (α, β) + C yδ l (α, δ l ) + C yδn (α, δ n ) + C y pa (α, β, pa ) (12)
C z = C z0 (α, β) + C zδm (α, δ m ) + C z pa (α, β, pa ) (13)
C l = C l0 (α, β) + C lδ l (α, δ l ) + C lδn (α, δ n ) + C l pa (α, β, pa ) (14)
C m = C m0 (α, β) + C mδm (α, δ m ) + C mpa (α, β, pa ) (15)
C n = C n0 (α, β) + C nδ l (α, δ l ) + C nδn (α, δ n ) + C npa (α, β, pa ) (16)
Every coefficient takes the form of a sum of non-linear subfunctions depending of a reduced number of variables. For example, C mδm represents the non-linear effects of the elevator angular deflection on the dimensionless pitch moment. Every subfunction has a corresponding database, data coming from wind tunnel experiments conducted at Onera Lille on a light aircraft scale model. During numerical integration of the dynamical system, evaluations of the subfunctions are done by linear interpolations in the corresponding databases.

pa = c V p cos(α) cos(β) + q sin(β) + r sin(α) cos(β) (17)
is the dimensionless form of the aerodynamic roll, component of the rotation vector along the wind axis (along the velocity of the aircraft). Databases corresponding to subfunctions dependent of pa come from experiments conducted in the vertical wind tunnel of ONERA Lille where rotary balance permits to get a conical motion of the aircraft around the wind axis (which is vertical in these experiments).

The equations of the system are

V = g(-sin Θ cos α cos β + cos Θ sin Φ sin β + cos Θ cos Φ sin α cos β) - 1 2m ρSV 2 C xa (18) α = q a cos β + g V cos β (sin Θ sin α + cos Θ cos Φ cos α) - 1 2m cos β ρSV C za (19) β = -r a + g V (sin Θ cos α sin β + cos Θ sin Φ cos β -cos Θ cos Φ sin α sin β) + 1 2m ρSV C ya (20) ṗ = B -C A rq + 1 2A ρScV 2 C l (21) q = C -A B pr + 1 2B ρScV 2 C m (22) ṙ = A -B C pq + 1 2C ρScV 2 C n (23) Θ = q cos Φ -r sin Φ (24) Φ = p + tan Θ(q sin Φ + r cos Φ) (25) (26)
with q a = p cos α sin β + q cos β + r sin α sin β (27) r a = -p sin α + r cos α (28)

C xa = C x cos α cos β -C y sin β + C z sin α cos β (29) C ya = C x cos α sin β + C y cos β + C z sin α sin β (30) C za = -C x sin α + C z cos α (31)
Traditionally, bifurcation diagrams are drawn with the equilibria of this system [START_REF] Carroll | Bifurcation analysis of nonlinear aircraft dynamics[END_REF]. Equilibria represent safe trajectories in the linear domain (at low angles of attacks), but represent potentially dangerous behaviours, as spins, in the post-stall domain (at high angles of attack). In the present paper, bifurcation theory is not used, but the joint UQ/GSA study is conducted on the temporal evolution of the state variables during the first instants. To integrate the system, a Runge-kutta method of order 2 is employed.

The system has potentially several bifurcations, which is motivating the use of temporal ME-ACD. Every state variable is a QOI ŷ depending on time and on the selected input parameters. Every of these QOIs is approximated by ME-ACD at every considered instant and partial variances are found with the analytical formula linking them to the coefficients of the local polynomial models.

Influence of damping parameters

With a model taking into account rotation effects with data only coming from rotary balance, a lot of equilibria of the system in the post stall region are not stable [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF]. It is shown that it can come from the lack of knowledge on the dependence of the aerodynamic coefficients in terms of q a and r a , components traducing the misalignment between rotation and velocity vectors. Experiments on a vertical wind-tunnel using a rotary balance only give information on dependences of the aerodynamic coefficients on p a , the projection of the rotation vector on the velocity direction. In [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF], authors add the so-called damping terms with information coming from forced oscillations along each of the 3 body axis of a similar aircraft scale model. Damping terms are obtained using a linear hypothesis, as usually q a ≪ p a and r a ≪ p a in spin motions. In this application, damping terms are added artificially in the Onera aerodynamic model. The subfunctions depending on these variables are also supposed linear, and only the dimensionless moments, having more impact, are modified. The added damping terms C lqa , C lra , C mqa , C nqa and C nra are the input parameters of the study. In [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF], the results of the forced oscillatory experiments are added with a change of frame to get the damping terms, as functions depending on α. In the considered approach, they are supposed constant, but taking random values following uniform distribution laws.

C l = C standard l + C lqa qa + C lra ra (32)
C m = C standard m + C mqa qa (33)
C n = C standard n + C nqa qa + C nra ra (34)
-

C lqa ∼ U(-0.1, 0.1) -C lra ∼ U(-0.1, 0.1) -C mqa ∼ U(-10, 0) -C nqa ∼ U(-0.1, 0.1) -C nra ∼ U(-0.1, 0.1)
The intervals of the distribution laws are estimated with the extrema of values of the forced oscillation data provided in [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF] with a wider range for an in-depth theoretical study of the impact of the damping terms.

The chosen initial conditions are the one corresponding to a spin, unstable without the damping terms

Hyperparameters ME-ACD N 10000 N E min 1 N Emax 10 k min 30 k max 60 n min 0 n max 5 δ min 10 -3 δ max 10 -3 k connect 100 N outliers 5 tested learning rates for the classifier [10 -m , 2 × 10 -m , 5 × 10 -m ], ∀m ∈ [1, 7] N h 5 N random 100
Table 1 -Selected hyperparameters for the study of the impact of damping parameters for an unstable spin.

-V 0 = 27 m.s -1 -α 0 = 59 • -β 0 = 5.9 • -p 0 = 1.8 rad.s -1 , q 0 = 0.58 rad.s -1 , r 0 = 3.03 rad.s -1 , leading to a rotation velocity of Ω = 3.6 rad.s -1 = 0.57tr.s -1 -Θ 0 = -32 • -Φ 0 = 13 • and control surfaces angular deflections are set to :

- 1 sums-up the selected hyperparameters for temporal ME-ACD and Figure 2 shows the temporal evolutions of the variance and of the partial variances. High variances at the last considered instants are due to the instability of the equilibrium with some combinations of input parameters and stability with another combinations. As it can be seen on the figure, C mqa and C nqa are the key parameters for the stability of this spin. Figure 3 is a plot of α at instant t = 20s vs the two key parameters. It underlines that for low values of C mqa and C nqa (mapped in [0, 1] on the figure), α keeps the initial value (59 • ), thus a proof of stability, and for higher values (mapped C mqa close to 1 and then real C mqa close to 0), α is not equal to the initial value at t = 20s. Figure 2 underlines the apparition of oscillations with different amplitudes dependent of the unstable scenarios and with pretty similar frequencies, as variances have pseudo oscillations with low minima. The approach permits to verify, with a relatively low number of run, that adding damping terms gives some stability to the studied spin. Moreover, the approach permits to identify that only 2 of the 5 suspected parameters have actually some impact on the considered spin. The limit values of the input parameters changing the stability of the spin are also available with the ME model. Some elements of the ME model at a given instant can be interpreted. For example, the orange element of Figure 3 C nra D (1) C nqa D (2) C nqa, C nra D (1) Cmqa D (2) Cmqa, C nra D (2) Cmqa, C nqa D (3) Cmqa, C nqa, C nra D (1) C lra D (2) C lra, C nra D (2) C lra, C nqa D (3) C lra, C nqa, C nra D (2) C lra, Cmqa D (3) C lra, Cmqa, C nra D (3) C lra, Cmqa, C nqa D (4) C lra, Cmqa, C nqa, C nra D (1) C lqa D (2) C lqa, C nra D (2) C lqa, C nqa D (3) C lqa, C nqa, C nra D (2) C lqa, Cmqa D (3) C lqa, Cmqa, C nra D (3) C lqa, Cmqa, C nqa D (4) C lqa, Cmqa, C nqa, C nra D (2) C lqa, C lra D (3) C lqa, C lra, C nra D (3) C lqa, C lra, C nqa D (4) C lqa, C lra, C nqa, C nra D (3) C lqa, C lra, Cmqa D (4) C lqa, C lra, Cmqa, C nra D (4) C lqa, C lra, Cmqa, C nqa D (5) C lqa, C lra, Cmqa, C nqa, C nra ( to stability.

δ l = -100% -δ m = -20 • -δ n = -20 • Table
Nevertheless, at some instants, the clustering of ME-ACD failed to find a suitable partition, due to the fact that the irregularity is located close to an edge of the parameter space. ME-ACD thus preferred a standard PC model with no partition of the parameter space, thus suffering of a lack of accuracy with the irregularity at these instants. ME-ACD is under development, and ameliorations of the method for concrete applications is the subject of active research. Adding more training observations in the zone of the irregularity might solve this issue.

Sensitivity and Robustness of aircraft's behaviours with respect to the aerodynamic model

An analysis of sensitivity and robustness on stable spins with respect to the aerodynamic model can be conducted with the considered approach. As the aerodynamic model is built on data, experimental uncertainties can affect the predictions of the model. Terms ∆C l , ∆C m and ∆C n can be added to the dimensionless moments. These terms can be taken as input parameters following stochastic distributions in their confidence interval. Variances of the state variables can then give information on the robustness of the behaviour with respect to the uncertainties (a low variance corresponding to a robust behaviour) and partial variances can give a quantitative information of the impact of every uncertainty and every interaction between the uncertainties. Cn D (1) Cm D (2) Cm, Cn D (1) Cl D (2) Cl, Cn D (2) Cl, Cm D (3) Cl, Cm, Cn Cn D (1) Cm D (2) Cm, Cn D (1) Cl D (2) Cl, Cn D (2) Cl, Cm D (3) Cl, Cm, Cn Cn D (1) Cm D (2) Cm, Cn D (1) Cl D (2) Cl, Cn D (2) Cl, Cm D (3) Cl, Cm, Cn The approach is used to study the impact of the uncertainties of the parameters of a Flight Dynamics system. As predicted aircraft behaviours rely on models of forces and moments with a chosen structure and evaluated through interpolation in experimental databases, they suffer from uncertainties. UQ and were conducted on the system and the impact of a modification of the model structure and the impact of the uncertainties of experimental data on a spin behaviour were studied. It permitted to underline, as in [START_REF] Farcy | Sensitivity of Spin Parameters to Uncertainties of the Aerodynamic Model of a Light Aircraft[END_REF], that the behaviour is highly impacted by these uncertainties. The approach permitted to determine that the added damping terms C mqa and C nqa had the most important impact on the stability of the studied spin.

ME method leads to elements which can be interpreted. In the above example, the element where the QOI is locally approximated by a constant polynomial highlights the combinations of values leading to stability of the spin.

The applications highlight that interpretation of some predictions of behaviours must be taken with precaution, as it is not robust to uncertainties and to lack of knowledge on the model of the aerodynamic coefficients. A temporal evolution is considered in this paper. The study of attraction domains [START_REF] Topcu | Robust region-of-attraction estimation[END_REF] of spins, coupled with UQ and GSA, could be another step forward in the study of dangerous aircraft's behaviours with an uncertain model. ME-ACD is under development and may lack of robustness in some practical applications. For example, in the applications of the paper, ME-ACD failed to find a suitable partition at some instants. ME-ACD thus preferred a standard PC model with no partition of the parameter space, thus suffering of a lack of accuracy with the irregularity at these problematic instants. Ameliorations of the algorithm to be more robust with practical applications are under development. Moreover, the approach can currently handle only a moderate number of input parameters (less than 9).

  set containing all the possible indices i indexing interactions of order j. It is the set of combinations without replacement of j elements of [|1, d|]. By convention, I (d) 0

  is the set of combinations with replacement of d elements of [|1, γ|] following the previous inequality.

Figure 1 -

 1 Figure 1 -Picture of the light aircraft scale model fixed on the rotary balance in the vertical wind tunnel of ONERA Lille. ©Onera

  of the coefficients defined at[START_REF] Carroll | Bifurcation analysis of nonlinear aircraft dynamics[END_REF], (15) and (16).

  represents the different combinations of values of the input parameters which lead

Figure 2 -

 2 Figure 2 -Temporal evolution of the variance and of partial variances of (a) V , (b) α and (c) β. (d) Legend.

Figure 4 -

 4 Figure 4 -On the left, temporal evolutions of the variables with all the combinations of input parameters taken as training observations and evolution of the expectation (a) V , (c) α and (e) β. On the right, temporal evolution of total and partial variances of (b) V , (d) α and (f) β

C m q a

Initial conditions are the same from previous section, and damping terms C m qa = -10 and C nqa = -0.1 are added according to the previous study to give some stability to the spin.

Selected input parameters are ∆C l , ∆C m and ∆C n and are supposed to follow a normal distribution law :

-∆C l follows a normal distribution of mean 0 and of standard deviation 0.04 3 . -∆C m follows a normal distribution of mean 0 and of standard deviation 0.08 3 . -∆C n follows a normal distribution of mean 0 and of standard deviation 0.005.

Figure 4 underlines that the considered aircraft behaviour is not robust to experimental uncertainties, even with the added damping terms. Uncertainties on the coefficients C n make the long-time stability of this spin uncertain. More experimental data or more accurate experimental data on this component of the aerodynamic moment are then necessary to be more sure of the actual long time stability of the spin (with the theoretically added damping terms). However, the state variables are constant during the first 10s, which means that it can be ensured that the spin has a short-time stability, and is therefore dangerous, even with the current uncertainties.

Conclusion and discussion

In this paper, an upcoming joint UQ/GSA method to estimate Sobol indices with a ME model is extended to solutions of ODE's. Considered ME method is the time-extension of ME-ACD algorithm [START_REF] Vauchel | A multi-element non-intrusive Polynomial Chaos method using agglomerative clustering based on the derivatives to study irregular and discontinuous Quantities of Interest[END_REF] and the