Maxime Louis

Benjamin Charlier

Stanley Durrleman

Learning Riemannian geometry for mixed-effect models using deep generative networks

We take up on recent work on the Riemannian geometry of generative networks to propose a new approach for learning both a manifold structure and a Riemannian metric from data. It allows the derivation of statistical analysis on manifolds without the need for the user to design new Riemannian structure for each specific problem. In high-dimensional data, it can learn non diagonal metrics, whereas manual design is often limited to the diagonal case. We illustrate how the method allows the construction of a meaningful low-dimensional representation of data and exhibit the geometry of the space of brain images during Alzheimer's progression.

Introduction

There has been a lot of interest into statistical methods for the analysis of manifold-valued data. These methods aim at proposing nonlinear alternatives to usual linear models when the data lies on a Riemannian manifold: the Fréchet mean [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] is an extension of the linear mean, principal geodesic analysis (PGA) [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF] is an extension of Principal Component Analysis (PCA), geodesic regression [START_REF] Fletcher | Geodesic regression on riemannian manifolds[END_REF] is an extension of linear regression and [START_REF] Srivastava | Statistical shape analysis: Clustering, learning, and testing[END_REF] proposes unsupervised clustering of manifold-valued data. In these approaches, the Riemannian manifold and its metric are fixed and typically encode the known constraints on the data, such as positive constraints or positive-definite matrices. These known constraints do not preclude an even more constrained structure of a particular distribution of points on this manifold, that we would like to estimate.

Consequently, a number of works [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Kruskal | Multidimensional scaling[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF][START_REF] Maaten | Visualizing data using t-sne[END_REF] offer unsupervised methods to perform manifold learning from data. However, only a few of these methods do obtain a parametrization of a submanifold of the space of observations -most of them do it implicitly such as [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Kruskal | Multidimensional scaling[END_REF][START_REF] Maaten | Visualizing data using t-sne[END_REF]-and almost none of them do estimate a Riemannian metric on the learned manifold: they often use the geometry induced by a standard metric on the space of observations. We argue here that the estimation of a Riemannian metric is an important part of the learning process which complements the manifold learning part. Indeed, there is no reason for the geometry induced by an usual metric on the data to be particularly relevant for machine learning tasks.

To gain a better understanding of the geometry of learning, there has been recent work on deep generative networks and on the geometry of the associated latent spaces [START_REF] Arvanitidis | Latent space oddity: on the curvature of deep generative models[END_REF][START_REF] Shao | The riemannian geometry of deep generative models[END_REF][START_REF] Hauser | Principles of riemannian geometry in neural networks[END_REF][START_REF] Chen | Metrics for deep generative models[END_REF]. Deep generative models -such as auto-encoders [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] or generative adversarial networks [START_REF] Goodfellow | Generative adversarial nets[END_REF]-learn an immersed submanifold of the observation space and a system of coordinates on this submanifold. In these papers, the geometry of the latent space is inferred a posteriori by pulling back the geometry of the observations onto the latent space. In [START_REF] Shao | The riemannian geometry of deep generative models[END_REF] the authors notice, experimentally, that the latent space is approximately flat and linear with this geometry.

Consequently, we propose in this paper to reverse the mechanic. Given a mixed-effect model with a Euclidean structure, we reformulate it on a low-dimensional Euclidean space and use a deep generative network to learn a mapping from this latent space to the observation space. The image of the network is a low-dimensional immersed submanifold that we equip with the push-forward of the Euclidean metric. We illustrate this procedure on two examples:

• An extension of Probabilistic Principal Component Analysis (PPCA) [START_REF] Tipping | Probabilistic principal component analysis[END_REF]: we formulate a nonlinear generalization. We show how this formulation is close to Probabilistic Principal Geodesic Analysis (PPGA) [START_REF] Zhang | Probabilistic principal geodesic analysis[END_REF].

• An unsupervised longitudinal model: we generalize the formulation proposed in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF] to allow inference of a submanifold of the space of observations and of a Riemannian metric on this manifold adapted to the unsupervised task.

In Section 2, we describe the construction of a latent space and of the push-forward of its Euclidean geometry onto the space of observations and we apply it to PPCA and to a longitudinal model. In Section 3, we describe the inference procedure used in the experiments presented in Section 4.

2 Geometrical model

Geometry of generative networks

We describe a generative neural network as a parametric function

Ψ w : U → V where U is an open subset of R d ,
V is an open subset of R D where d, D ∈ N with typically d D and w are the neural network weights. In [START_REF] Shao | The riemannian geometry of deep generative models[END_REF], the authors show that Ψ w (U ) is a d-dimensional immersed submanifold of R D if the activation functions of each layer are smooth and monotonic and if the weight matrix of each layer has maximal rank. The first condition is easy to enforce and we check the second condition after training. Note that an immersed submanifold is not a submanifold in general, but it is locally a submanifold.

A Riemannian metric on a smooth manifold M is a smoothly varying inner product on the tangent bundle T M. In [START_REF] Shao | The riemannian geometry of deep generative models[END_REF][START_REF] Chen | Metrics for deep generative models[END_REF], the authors illustrate how to pull-back a metric on R D to the latent space R d . They note, experimentally, that the induced metric on the latent space is almost flat. Based on this empirical observation, we postulate that there exists a latent space where the pull-back of the Riemannian structure describing the data is Euclidean.

As a consequence, we propose to go the opposite way by pushing-forward the Euclidean metric of a latent space onto a submanifold of the observation space. Doing so enforces by construction the flatness of the latent space during learning. Let g be the Euclidean metric on U , we can define the push-forward of g on Ψ w (V ). For any smooth vector fields X, Y on Ψ w (U ), it is defined as:

Ψ * w (g)(X, Y ) = g((Ψ w ) * (X), (Ψ w ) * (Y ))
where (Ψ w ) * (X) and (Ψ w ) * (Y ) are the pull-back of X and Y on U defined by

(Ψ w ) * (X)(f ) = X(f • Ψ -1 w ) for any smooth function f : U → R. If γ : [0, 1] → U is a geodesic on (U, g), then Ψ w • γ : [0, 1] → V is a geodesic on (Ψ w (U ), Ψ * w (g)).
Finally, for any p on a manifold M, we recall that the Riemannian exponential map Exp p is defined on a neighborhood of 0 in T p M by Exp p (v) = γ(1) where γ is the unique geodesic with γ(0) = p and γ(0) = v.

Note that the function Ψ w parametrizes:

• A submanifold Ψ w (U ) of the space of observations • A metric Ψ * w (g) on this submanifold
We show in the next section how to adapt mixed-effect models which assume a linear structure of the observation space into models which consider a similar structure on a submanifold of the observation space described by a latent space and a neural network Ψ w .

2.2 Generalizing mixed-effect models to Riemannian manifolds.

We consider generative models of the form:

y i = F x i + Gu i + ε i (1)
where y i ∈ R D is an observation, x i ∈ R f is known, u i is a latent random variable in R g , F and G are unknown matrices, and ε ∼ N (0, Σ ε ) is a Gaussian noise. We also suppose a normal prior on the latent variable: u ∼ N (0, Σ u ). We propose nonlinear generalizations of two different models of the form (1): PPCA and a longitudinal model.

Generalizing PPCA. As shown in [START_REF] Tipping | Probabilistic principal component analysis[END_REF], PCA is the Maximum Likelihood (ML) estimate of the PPCA model which writes y i = Gu i + ε i where u is the latent position with prior u ∼ N (0, I d ), G is a D × d matrix containing the principal directions and ε ∼ N (0, σ 2 I D ). This is a model of the form (1) with no covariates x i . We propose to generalize the model to:

y i = Ψ w (u i ) + ε i . (2)
where Ψ w is a parametric family of C 1 immersions. PPCA describes each observation as a point on a linear subspace with coordinates u i while the proposed formulation (2) describes each observation as a point on a submanifold of R D with coordinates u i . G has been removed on the new model, as the neural network is flexible enough to convert the isotropic unit normal distribution of the latent variable u into a distribution on the submanifold Ψ(R d ) which is close to the data distribution.

To illustrate the similarity with PPGA, we denote J = Jac Ψ w (0) the Jacobian of Ψ w and write, for any observation y i close to Ψ w (0):

y i = Ψ w (u i ) + ε i = Exp Ψw(0) ((Ψ w ) * (u i )) + ε i = Exp Ψw(0) (Ju i ) + ε i . (3)
where the first step is by construction of the Euclidean latent space R d and the second is by definition of (Ψ w ) * . If Ψ w was fixed -which amounts to fixing the manifold structure-the estimation of J by ML in ( 3) is almost equivalent to the PPGA proposed in [START_REF] Zhang | Probabilistic principal geodesic analysis[END_REF] (to the detail of the structure of the noise, which we assume isotropic Gaussian in R D when it is assumed isotropic Gaussian on M for the PPGA). One can also note the proximity with a variationnal auto-encoder as described in [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. Finally, model ( 2) should be compared to the Locally Adaptive Normal Distribution (LAND) proposed in [START_REF] Arvanitidis | A locally adaptive normal distribution[END_REF], which proposes to learn a Riemannian metric so that the observations are distributed along a normal distribution for this Riemannian metric. Indeed, after inference of the model ( 2), the Riemannian metric will have been learned so that the observations are approximately distributed according to a normal distribution on the manifold (Ψ w (U ), Ψ * w (g)), as is the case for LAND. But our model goes beyond LAND by learning a submanifold of the space of observations and a metric which is not necessarily diagonal, while staying computationally efficient in high dimension.

Generalizing a model for the analysis of longitudinal data. In [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF], the authors propose a framework for the analysis of manifold-valued trajectories. They assume an a priori Riemannian geometry on the space of observations. We denote (y ij ) j=1,...,ni ∈ R D the observations of the subject i, measured at times (t ij ) j=1,...,ni . In the case of a linear manifold, the model writes:

y ij = v 0 Φ i (t ij ) + As i + ε ij (4)
where

v 0 ∈ R D , Φ i (t) = exp(η i )(t -τ i )
is the time reparametrization of the individual trajectory: for disease modelling, exp(η i ) > 0 controls the pace of progression of the subject and τ i ∈ R controls the time-shift for the subject progression. A is a D × d matrix with d ≤ D -1 which contains directions orthogonal to v 0 , which allow to take into account different positions of the trajectories for different subjects. The latent variables are

u i = (η i , τ i , s i ) with priors η i ∼ N (0, σ η ), τ i ∼ N (0, σ τ ) and s i ∼ N (0, I d-1
). σ τ and σ η are initialized from the data and fixed during the estimation. Ignoring the time reparametrization Φ i , the model ( 4) is of the form (1). We use a similar generalization procedure as for (2) and write the new generative model:

y ij = Ψ w e 1 Φ i (t ij ) + d l=2 s il e l + ε ij (5)
where (e 1 , . . . , e d ) is the canonical basis of R d . t → Ψ w (e 1 t) should be interpreted as the geodesic of mean progression, while the directions (Ψ w ) * (e i ) for i ∈ {2, . . . , d} capture variability between subjects. As for PPCA, v 0 and A have been removed from the new model since the neural network is flexible enough to send e 1 onto a main direction of progression and (e 2 , . . . , e d ) onto vector fields which are orthogonal to this direction of projection. Figure 1 summarizes the construction. Note that as in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF], the subject trajectories are obtained by a parallel shift of a mean geodesic in the directions e 2 , . . . , e n on the manifold described by Ψ w . This can be interpreted as the Riemannian analog of linear translation.

Both generalizations share the similar idea of keeping the structure of the initial model in a lowdimensional space R d and using a neural network to transport this structure onto a submanifold of the space of observations. This general principle of converting the linear structure into an equivalent nonlinear structure can be adapted to numerous models such as Probabilistic Linear Discriminant Analysis [START_REF] Prince | Probabilistic linear discriminant analysis for inferences about identity[END_REF], Random slopes/intercepts models or linear regression.

Finally, in both cases, Ψ w being an isometry between (R d , g) and (Ψ w (R d ), Ψ * w (g)), any computation of Riemannian exponential, logarithm or parallel transport can be done in R d before being pushforward to the observation space. Thus, these operations, which can be costly in a Riemannian setting, are inexpensive in our setting. Besides, the analysis of the data can be done directly in the latent space R d whose Euclidean structure is faithful to the nonlinear geometry of the data.

Inference

We describe here the inference procedure for models [START_REF] Arvanitidis | Latent space oddity: on the curvature of deep generative models[END_REF] and [START_REF] Chen | Metrics for deep generative models[END_REF]. The inference consists in finding the Maximum A Posteriori (MAP) of a directed probabilistic model with latent variables u. A possibility is to use variationnal Bayes approaches such as [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. An alternative is to use the Expectation-Maximization (EM) algorithm. The E step requires the computation of integrals of the form u log (p(y|u, θ)) p(u, θ k )du which are intractable in our case, so we resort the the Stochastic Approximation EM (SAEM) [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] which alternates:

• Simulation. For each observation y i , generate u i , a realization of the hidden variable under the posterior density p(u|y i , θ k ).

• Approximation. Update Q k (θ) = Q k-1 (θ) + γ k ( i p(y i |u i , θ) -Q k-1 (θ)) • Maximization. Set θ k+1 = argmax θ Q(θ).
where during L burn-in iterations γ k = 1 and then for all k > L,

0 ≤ γ k ≤ 1, ∞ i=1 γ k = ∞ and ∞ i=1 γ 2 k < ∞.
Once again, this procedure is intractable in our case, since the maximization step cannot be computed at a reasonnable cost. We therefore replace the Approximation step by simply setting Q(θ) = i p(y i |u i , θ), which can be optimized by stochastic gradient descent. This amounts to keeping only the burn-in phase of the SAEM (γ k = 1) which is, as we noted empirically, the most important phase with respect to space exploration of the individual variables u i .

Simulation-Expectation

To perform the Simulation step, we use the symmetric Hasting-Metropolis sampler [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], a Markov Chain Monte Carlo method. We run 25 iterations of the Markov Chain for each simulation, to limit samples correlation.

Maximization The maximization can be performed by stochastic gradient descent on Q with respect to the neural network weights w. We run ten epochs of gradient descent at each maximization, using Adam [START_REF] Kingma | A method for stochastic optimization[END_REF]. The noise variance σ ε can be updated using a closed-form expression derived from the log-likelihood for both models (2) and [START_REF] Chen | Metrics for deep generative models[END_REF]. Performing the gradient descent repeatedly with slightly changing latent variables seems to have a stabilizing effect on the neural network learning: it fills the latent space around each subjects and can be interpreted as a data augmentation procedure.

Algorithm 1 Inference procedure for models (2) and ( 5) input :Observations (y i ) i , fixed effects (β i ) i , initial parameters θ 0 = (w 0 , σ 0 ε ) and samples u 0 i . output :Estimation of θ MAP and of samples u i distributed following p(u | (y i ), θ MAP ) i . Maximization: Update σ ε using the closed form expression.

Experiments

The network architectures used in the experiments which follow are standard and given in the supplementary materials. We only use C 1 transfer functions to obtain a family of C 1 functions.

A python code for the experiments and the data sets will be made available upon publication of the paper. Note that both the E and the M steps are massively parallelisable: the samples for each subjects are independent and the neural network training can easily be performed on multiple CPUs or GPUs. 200 iterations of the proposed method take an average of 1h30 on 4 CPUs for the longitudinal model, and 10 minutes for the generalized PPCA on a subset of 1000 digits from MNIST.

Generalized PPCA on MNIST

We estimate the model (2) on 1000 randomly selected digits from MNIST [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. The mean squared error (MSE) for different dimensions d are shown on Figure 3. The proposed approach consistently beats PCA on both train and held-out sets. To analyze the effect of the Riemannian metric learning, we compare four dimension reduction methods: the first two PCA components, Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] on the raw images, t-SNE [START_REF] Maaten | Visualizing data using t-sne[END_REF] and MDS [START_REF] Kruskal | Multidimensional scaling[END_REF] performed on the latent positions u of the images after estimation. MDS uses the Euclidean distances on the latent positions u of the observations, which amounts by construction to using the geodesic distances between the observations. It is therefore a direct representation of the learned Riemannian manifold structure. Figure 2 shows the results. All methods are run with default parameters. Interestingly, although the learning is unsupervised, the manifold geometry did capture class distributions in a much more clustered way than the other methods. This underlines the importance of estimating a metric on the manifold and of not relying on usual metrics which may not be adapted to the geometry of the data set.

Generalized longitudinal model

Recovering a synthetic geometry

We perform experiments of model ( 5) on a synthetic data set. Each observation is a 64 × 64 gray-level image of a cross fully described by 3 parameters: the angles of the right and left arm and the length of the arms. The arm angles are drawn along a zero-centered normal distribution. A mean progression scenario is prescribed for the arm lengths. Figure 5 shows synthetic subjects. Note that each image is defined by a finite number of parameters with respect to which it varies smoothly: hence the generated set of images belongs to a 3-dimensional submanifold of the set of 64 × 64 grey level images. We generate thirty different data sets with 200 subjects and 7 observations per subject according to this distribution, 25% of the subjects are kept in a test set. We add a white noise on the images with standard deviation of 0.025, 0.05 or 0.075. We estimate the model with d = 3 on the train set. Figure 6 show the obtained geodesic progression and two different patterns Ψ w (v 0 t + e i ) i=2,...,3 . Figure 4 (left) shows the estimation of the noise variance, averaged over the ten folds for each noise level.

We then fit the test observations onto the obtained model by gradient descent on the individual parameters u. We compare the MSE of reconstruction between train and test set, to check if it generalizes well. We also compute the MSE between the reconstructed images and the original images without the added noise to see that the generative network has a denoising effect. The results are given on Figure 4 (right). Interestingly, for all the noise levels, the model, which was trained on the noisy images, recovers equally well the original images without noise, suggesting a resilience of the model with respect to noise component outside of the data manifold. These experiments indicate the ability of the model to capture the synthetic submanifold and its Riemannian structure.

Cognitive scores

We use the cognitive scores grading the subjects memory, praxis, language and concentration extracted from the ADNI database as in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF] and run the model to estimate the geometry adapted to the progression of the cognitive scores. We emphasize that in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF], and in its derivations for other types of data [START_REF] Bône | Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms[END_REF][START_REF] Koval | Statistical learning of spatiotemporal patterns from longitudinal manifold-valued networks[END_REF], the authors consider user-defined metrics specifically manufactured for the data set. Here, we propose a generic way to learn both the manifold and metric regardless of the kind of data considered. We do recover a geometry which is similar to the one postulated in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF] with logistic-like geodesics, as shown on Figure 7, although this geometry was not prescribed a priori. In addition, our generalization allows more flexibility with different behaviors for each score and non-parallel geodesics. The estimation of a non-diagonal matrix allows a more complicated form of the parallel curves t → Ψ w (e 1 t + e i ) and different shapes of these parallel curves for the different bio-markers as visible for memory and language on Figure 7. We use Magnetic Resonance Images from the ADNI database. We select the images from subjects which ultimately develop the Alzheimer's disease and are at least observed twice. A common 128 × 128 slice is extracted from all rigidly aligned images, chosen so as to contain the Hippocampi and the Amygdala, which are particularly affected during Alzheimer's progression. We perform a 10-fold cross-validation on the subjects. After estimation, the test set of subjects is fitted onto the model by gradient descent on the individual parameters u i .

Disease model progression

A mean trajectory and the patterns of progressions Ψ w (v 0 t + e i ) i=2,...,d are given on Figure 8. Those different patterns correspond to parallel trajectories to the reference geodesic on the learned submanifold as shown on Figure 1. The mean trajectory recovers grey matter loss during the disease progression, and especially grey matter loss near the Hippocampi and ventricles.

The different patterns of progression on Figure 8 illustrate slight differences in the global pattern, which correct for the different individual anatomies as well as for different progression patterns.

Reconstructing observations. To illustrate the reconstruction of the images by the generative model, we provide on Figure 9 the original images for a subject as well as its reconstructions. The reconstructed images are smoother and blurrier than the original images, since the model discarded what it sees as noise in the manifold estimation and since the cost is an 2 distance. To check whether the model generalizes to unseen observations, we fit, for each fold, the test set of observations to the model by gradient descent on the latent variables u i . The MSE on the train set is of 0.087 ± 0.0008 while it is of 0.088 ± 0.003 on the test set, which shows that the model did not overfit.

Changing the latent space dimension.To analyze the effect of the dimension d of the latent space, we repeat the estimation of the model with d varying between 2 and 22. Figure 10 shows the results. The model MSE decreases steadily until d ∼ 18, where it stagnates. This may indicate that the MRI data lies on a 18-dimensional submanifold of the whole space of observations. Analyzing the model output. We compare the distributions of the latent variables between subpopulations having different ages of diagnosis, alleles of the APOE gene -which indicates strong predisposition to Alzheimer's disease-or gender. For all the folds, we found a negative correlation between the pace of progression of the disease and the age of diagnosis, as observed in the literature [START_REF] Bigio | Synapse loss is greater in presenile than senile onset alzheimer disease: implications for the cognitive reserve hypothesis[END_REF], a pace of progression higher when at least one allele of the APOE gene is present, as described in [START_REF] Sadigh-Eteghad | Association of apolipoprotein e epsilon 4 allele with sporadic late onset alzheimer's disease. A meta-analysis[END_REF], and a higher pace of progression for female subjects than for male subjects [START_REF] Farrer | Effects of age, sex, and ethnicity on the association between apolipoprotein e genotype and alzheimer disease: a meta-analysis[END_REF]. This shows that the Riemannian modelling proposed in the paper does yield informative results on real data.

Discussion

We show how to generalize two linear mixed-effect models to learn a Riemannian manifold of the space of observations and a Riemannian metric on this manifold which optimize the likelihood of the model. We show the importance of learning a Riemannian metric that is adapted to the model on top of learning a submanifold of the space of observations. This generalization procedure may be applied to a variety of generative linear models. The generalization of the mixed-effect model for longitudinal data is promising in the perspective of unsupervised disease progression model, since it identifies a geometry of progression for the images and a metric which may be of relevance for tasks such as early diagnosis. Further analysis of the obtained geometry such as classification with the learned metric are natural continuations of this work. Because of the huge parametric family of considered manifolds and metrics parametrized by the neural network, the identifiability of the parameters of the proposed models remains an open question, shared by all current deep learning methods.

Figure 1 :

 1 Figure 1: Transporting the Euclidean geometry onto a submanifold of the observation space using the neural network Ψ w . t → Ψ ( e 1 t) is the geodesic of mean progression.

1 repeat 2 3

 23 Simulation-Expectation: Draw a candidate u c i from the proposal distribution p b (.|u k i , y i ) for all iCompute the acceptance ratio τ =p(yi|θ k ,u c i )p(u c i |θ k ) p(yi|θ k ,u k i )p(u k i |θ k )and accept u c i with probability min(τ, 1).

4

 4

5 Update w k by stochastic gradient descent: 6 for epoch ← 0 to 10 do 7 for 9

 5679 batch b in (y i , u i ) i do 8 Compute Q = y,u∈b p(y|u, θ k )Update w by stochastic gradient descent w = w -α • ∇ w Q or using Adam.

Figure 3 :

 3 Figure 3: MSE with standard PCA and model (2) for varying dimension d.

Figure 5 :

 5 Figure 5: Each row is a synthetic subject.

Figure 6 :

 6 Figure 6: Top row is t → Ψ w (e 0 t). Following rows are t → Ψ w (e 0 t + e i ) for i ∈ {2, 3}.

Figure 7 :

 7 Figure 7: Learned main progression of the cognitive scores. Dashed curves are t → Ψ w (e 0 t ± e 1 ): the first estimated parallel direction.

Figure 8 :

 8 Figure8: Top row is the reference geodesic t → Ψ w (e 0 t). Following rows are t → Ψ w (e 0 t + e i ) -Ψ w (e 0 t) for i ∈ {2, . . . , d}. Grey matter loss is visible in the main progression, with some variations in the parallel curves.

  0.1 0.2 0.3 0.4 0.5

Figure 9 :

 9 Figure 9: Top row: original trajectory. Bottom row: reconstruction by the model.

Figure 10 :

 10 Figure 10: MSE with varying latent space dimension.

Preprint. Work in progress.


Joël Serra

Jean-Emmanuel E Pierré

Jean-Charles C Passieux

Jean-Noël N Périé

Christophe Bouvet

Bruno Castanié
email: bruno.castanie@insa-toulouse.fr

Caroline Petiot

Validation and modeling of aeronautical composite structures subjected to combined loadings: The VERTEX project. Part 2: Load envelopes for the assessment of panels with large notches

Keywords:

come

Validation and modeling of aeronautical composite structures subjected to combined loadings: The VERTEX project. Part 2: Load envelopes for the assessment of panels with large notches

Introduction

The general context of the validation of composite aeronautical structures was recalled in the first part of this article [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF]. Certification is still largely based on an experimental, semi-empirical approach called the ''pyramid of tests" (Rouchon [START_REF] Rouchon | Certification of large airplane composite structures, recent progress and new trends in compliance philosophy[END_REF]). This approach is applied to structural sizing in general (common areas, structural details, junctions, etc.), taking into account nondetectable low energy impact damage (below the Barely Visible Impact Damage -BVID) or exceptional damage. In the latter category, aircraft manufacturers must justify their structure under static loads with a large crack that could occur as a result of an exceptional event. Moreover, for certification purposes, this crack must pass through a fuselage frame and is commonly named a ''two-bay crack". In general, the certification tests are carried out at full scale with internal pressure as shown in the Airbus test (Fig. 1). As part of aeronautical programs, tests of sandwich panels with large cuts have also been carried out by means of the FASTER Fixture Test [START_REF] Leone | Structural testing and analysis of honeycomb sandwich composite fuselage panels November[END_REF]. These tests were widely instrumented by a combination of DIC (Digital Image Correlation) analysis, infrared thermography and acoustic emission. They were carried out under various combinations of biaxial tension and internal pressure. The size of the panels was several square meters. Cross-cutting tests were also carried out on PRSEUS type structures on the same system [START_REF] Bergan | Full-scale test and analysis of a PRSEUS fuselage panel to assess damage-containment features[END_REF].

However, the majority of academic studies on notches are performed at the coupon scale for practical reasons. From an application point of view, due to the complexity of crack tip failure modes, many holes or cracks are sized using point-stress, by calculating the damage at a distance d 0 from the stress concentration [START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments -a review[END_REF][START_REF] Tan | Stress concentrations in laminated composites[END_REF]. For this model, the stress is averaged over this distance d 0 , called the ''average distance", in order to smooth the stress field in the surroundings of the hole/notch. One of the important issues generated by the study of holes and cracks is the scaling effect [START_REF] Bazant | The scaling of structural strength[END_REF][START_REF] Bazant | Size effect and fracture characteristics of composite laminates[END_REF][START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF][START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF][START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF], which can be rendered by ''averaged stresses" type models even if the distance d 0 can depend on several factors [START_REF] Tan | Stress concentrations in laminated composites[END_REF], such as stacking sequence, material, and hole/notch size. The size effect leads, in particular, to a small hole being less penalizing than a large hole (for homothetic size of structures). Numerous experimental and numerical investigations have been conducted in order to explain the physical meaning of the size effect [START_REF] Bazant | The scaling of structural strength[END_REF][START_REF] Bazant | Size effect and fracture characteristics of composite laminates[END_REF][START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF][START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF][START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF]. For example, Wisnom et al. [START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF] showed that the strength of notched composite laminates decreased when the notch size increased and concluded that the hole size effect was triggered by the presence of non-critical damage at the ply level, such as fiber failure, ply splitting, delamination or matrix cracking in the vicinity of the hole, smoothing the stress concentration. The ratio of the size of this ''fracture process zone" to the specimen size explains the strength difference between small and large specimens. Recently, by using our own modeling strategies known as ''DPM" for Discrete Ply Modelling, we succeeded in capturing this size effect on open hole tensile specimens [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF] and also on notched specimens [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. This strategy was initially developed to simulate low energy impact on laminates [START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF][START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF] or pull-through [START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF]. Thanks to its unique capability to capture the failure scenario, the discrete nature of matrix cracking and inter-and intra-laminar damage coupling, this strategy was recently extended to in-plane issues [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF]. In the present paper, this strategy will be extended to the issue of large notches under complex loading and will be recalled only very briefly below.

The VERTEX test method and the FE-SDIC (Finite Element-Stereo Digital Image correlation) measurement strategy have been presented in the first part [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] and are based on the authors' past experience on complex loadings and structural testing [START_REF] Castanié | Theoretical and experimental analysis of asymmetric sandwich structures[END_REF][START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF]. Details of the VERTEX collaborative project are given in [20]. Tests carried out on technological specimens and under complex tensile, shear and combined loadings using the VERTEX method are presented in the first part of this paper. A first approach to the calculation/tests dialog on double thickness specimens under tensile and shear loading is then developed. Eventually, new tests called ''envelope" tests are proposed and performed. They allow a load path to be followed with the aim of validating the behavior of the notched structure. Prospects regarding research directions and in terms of approaches for the certification of composite structures are also proposed.

Materials and methods

Description of specimens and tests

The stacking sequences were the same as in [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. The plies were made of Hexcel's T700-M21 carbon epoxy unidirectional laminate with a nominal thickness of 0.125 mm. Two specimens were also manufactured with a double thickness ply of 0.250 mm. Three symmetrical stacking sequences of 13 plies were studied. The three different stackings (C3-1, C3-2, C3-3) presented the same number of plies in each direction (0�, 90� and ±45�), and only the relative position of the plies changed between layups:

-C3-1 [45/45/X/X/X/90/0/90/X/X/X/45/45] -C3-2 [X/X/X/X/0/90/0/90/0/X/X/X/X] -C3-3 [X/X/X/X/X/0/0/0/X/X/X/X/X] Specimen dimensions were 560 540 mm 2 with a central area of interest of 400 400 mm 2 [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF]. A central crack 100 mm long and 2 mm wide was milled in the center of the specimen [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. The preliminary test carried out in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] made it possible to demonstrate that the fasteners sustained the loads without the need for reinforcement plies and therefore the test specimens had a constant thickness. However, unlike in the preliminary test, aluminum plates were added here to maximize load transfer through friction. The loading principle of the VERTEX assembly is shown in Fig. 2. A view of the test specimen mounted on the VERTEX assembly is shown in Fig. 3.

The list of tests carried out is presented in Table 1. Four tests were carried out for each laminate: tensile, shear, combined and envelope. Two additional tests were carried out on C3-1 with a double thickness ply (0.25 mm instead of 0.125 mm) and were used for the dialog between tests and calculations. These specimens are called ''C3-1d". Tensile and shear tests were controlled by imposed displacements as indicated in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF]. For the combinedloads test, the four actuators were driven simultaneously in pairs. Based on the tensile and shear displacements, an average linear loading path was imposed (Fig. 4). The case of envelope loads was more complex and is detailed in Section 5.

Monitoring strategy

A speckle was printed on the specimens and was digitally defined to optimize the resolution (Fig. 5). For the tests carried out, three pairs of cameras were installed (Fig. 5): a pair of 5 Mpx Pike cameras acquired images of the whole plate (red), a second pair focused on one end of the cut (green) and the third pair, of Prosilica 29 Mpx cameras, acquired images of the entire cut (blue). Given the differences in sensors and lenses used for notch edge Fig. 6. Tension and Shear Shell Forces in the specimen (Cases 2, 3 and 6, see Table 1).

C3-1 C3-2 C3-3 Traction X + XX X X Shear X + XX X X Tension + Shear X X X Envelope Curves X X X
measurements, the resolutions of the zones shown in blue and green were approximately the same. The ''green" pair therefore merely offers redundancy in the useful measurement for the dialog between the test and the calculation. As explained in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] and in [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF], the numerous structural redundancies render an in situ measurement of the specimen inevitable. The mean stresses were determined by the theory of the laminates from the strains obtained by SDIC (here VIC3D), averaged over the red zones (Fig. 6) at the edge of the zone of interest of the specimen. The parameters used for stereocorrelation are those commonly used. The displacement was measured using subsets of 29 29 px 2 . The distance between two consecutive measurements was 7 px. To compute strains, 15 15 displacement measurements were used and the strain was thus estimated using an optical gauge of size J v = N step (N filter 1) + N subset = 127 px. Knowing the resolution of the cameras (5 Mpx), the size of the field of view (350 350 mm 2 ) and the size of the virtual gauge: 127 px, the surface in the coordinate system of the specimen whose projection corresponds to this virtual gauge was estimated to be 33 33 mm 2 . It seems sufficiently small to account for the (small) variation of the gradients in the observed zones (red in Fig. 6). The ) were determined by dividing the displacements measured at extreme points (in blue in Fig. 6) by the distance separating them. The imposed shear was calculated in the same way (Fig. 6). The behaviors of the three laminates studied (C3-1, C3-2 and C3-3) were analyzed according to the same methodology for all types of loading: Tension, Shear, Tension/Shear and Envelope. Another method for determining the stress resultants involved in the loading of a specimen by the VERTEX assembly would have been, firstly, to measure the displacement field on

Table 2

Comparison between VERTEX tests and open-hole and notched specimen under tension.. the edge of the plate by stereocorrelation, secondly, to impose these boundary conditions on a model and, thirdly, to recover the reaction forces associated with each constrained node in motion. After some preliminary tests, this procedure was not implemented in this first phase of the VERTEX studies, mainly because of problems of unresolved border regularization in the FE-SDIC method at that time. The other method described, more robust, was therefore preferred.

Specimens C3-1 C3-2 C3-3 � r 1 (%) CV (%) � r 1 (%) CV (%) � r 1 (%) CV (%) ''

Tests results and analysis

Tensile loading tests

Stress/strain curves obtained in tension using the methodology explained in the previous section are shown in Fig. 7. The vertical arrows represent the final failure of the laminates, corresponding to cuts that propagated to the edges of the test piece. The displacement control of the actuators associated with the low acquisition frequency used (1 image/s for the stereocorrelation cameras) generated oscillations on the curves in Fig. 7 and on the others obtained using the same method presented hereafter (Fig. 8 and followings). Using the in situ measurement pair (r xx and e xx ), linear curves were obtained until the damage occurred. In this way, the phases of initial positioning and the possible non-linearities of the assembly were filtered. The specimens were of the same stiffness, 38 GPa, close to the 43 GPa determined experimentally on a notched coupon [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. This decrease in stiffness was probably due to the ratio of notch width to specimen width, which was higher for VERTEX (1/4) specimens than for coupon-scale (1/6) specimens. This reduction in stiffness was counterbalanced by an increase due to the specific boundary conditions of the VERTEX assembly: the 4 sides of the test piece were clamped, which, by the Poisson effect, theoretically increased its stiffness by approximately 10% (1/ (1 m 2 ) with m 0.3).

Localized bending was observed near the notch (Fig. 8, Layout C3-1). It accelerated the propagation of the crack and thus reduced the maximum failure stress. Points F, indicated by the name of the draping shown in Fig. 7, corresponded to a maximum deflection of the notch border approximately equal to the thickness of the laminate (single thickness: 1.625 mm). The same reference (1.625 mm) was used (in value) for laminates C3-1d. The double thickness laminate, C3-1d, being stiffer in bending, required a higher load to reach similar bending. The failure stress observed for the double thickness specimen was therefore be greater than those observed for the ''single thickness". The direction of the local bending observed for the C3-1 laminate (deflecting towards the outside of the box), whatever the thickness of the plate, was not the same as that observed for laminates C3-2 and C3-3 (towards the inside of the loading box). The phenomenon thus seems very unstable and may depend on the initial imperfections of the plate. To verify the latter hypothesis, several tension tests should be carried out on the same laminate and the imperfections measured. It is also possible that this nonlinear geometrical phenomenon was due to the introduction of a local bending moment at the junction between the frame and the upper part of the central box formed by the composite specimen, a problem already identified in [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF]. Nevertheless, this behavior was not detected on preliminary tests on aluminum plates [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF].

The first failures visually identified are indicated by the points R in Fig. 7. For the C3-1 laminate, the failure was initiated at the notch tip and appeared just before total laminate failure. Cracking of the matrix oriented at 45� on the right side of the notch (Fig. 9), was detected (R C31 ). The main difference between the rupture scenario of the VERTEX test specimens and those of the simple test specimens studied in [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF] were the influence of boundary conditions. The tension applied in the Vertex test bench is ''quasi pure" because the plate is bolted on the 4 sides. Also, the edges of the notch are relatively close to the bolted areas. In this case, the principle of Saint-Venant is not satisfied, and this is what happens in most ''structural" tests. Nevertheless, this type of complex test is more representative of the real structures, such as a stiffened fuselage and provides more information than classical tests performed on coupons (Fig. 1). The final failure patterns of the VERTEX specimens are shown in Fig. 10. For laminate C3-1, the pattern is similar to those observed on notched coupons [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. It resembles a pull-out type with a crack propagation direction inclined at a few degrees from the horizontal. The laminates C3-2 and C3-3 appear to be damaged along different directions: the crack propagates in a direction oriented at 45� and then almost orthogonally to the direction of tension, which is unusual in structural tests [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF] but is similar to what is observed in open hole tensile tests [START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF]. It is interesting to compare the VER-TEX tests with the traditional coupon tests carried out in [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. The different tensile stresses for all the tensile specimens studied in this work are summarized in Table 2. The standard test curves are shown in Fig. 11. For the VERTEX tests, we find the following relative order of stresses at failure: r 1 C32 ffi r 1 C31 < r 1

C33 . The failure stresses of laminates C3-1 and C3-2 are very close. This confirms that the C3-3 laminate has the best tensile strength, as already observed for the coupons. For the large cut issue, the VER-TEX test exhibits structural effects: influence of boundary conditions and probably of initial imperfections of the specimens, which generate localized bending around the notch. These effects produce a significant decrease in the failure strength, which must be taken into account when designing the structure. In addition, the difference between the failure stresses obtained for the ''Notched" and the ''VERTEX notched" specimens is increased by the difference between the ratios of ''failure length to width of the specimen" (respectively 1/6 and 1/4). It also should be noted that the stiffnesses of VERTEX specimens are very similar to those of the ''coupon" type specimens (Fig. 11) thus validating the method of determining the stress / strain curves from the DIC.

Shear loading tests

The shear stress/strain curves are shown in Fig. 12. The vertical arrows represent the final failure of the laminates corresponding to cuts that have propagated to the edges of the test piece. In the same way as in the study of laminates subjected to tension, specific experimental points, noted F, denote a deflection of the edge of the notch here, equal to approximately 1.625 mm, and the points R represent the first damage visually observed (on the surface). Because of the large area of the specimens and their small thickness, the shear buckling occurs quite early and explains the position of points F and the first inflection of the test curves. Buckling occurs a little later in C3-3 and C3-2 than in C3-1 due to the presence of external plies at ±45�. The ''double thickness" plate C3-1d logically buckles even later because of its greater stiffness. The shear buckling can be clearly seen at 45� in Fig. 13. From the points R, the behavior of the test pieces shows a strong interaction between the propagation of the cracks and the postbuckling behavior which is very difficult to describe (snap-through). Reference may be made to one of the films of the shear tests (see [21]). As for tension, the first visible damage is a crack at 45� at the notch tip (Fig. 14). The crack propagates slightly but interacts with the postbuckling [21] until sudden final failure. The shear failure patterns for the three laminates are shown in Fig. 15. The failure scenarios and patterns are substantially identical. Crack propagation at 45� is expected in the case of a shear test [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF], which is the case for C3-1 and C3-2 -but not for C3-3, where propagation is rather in the 0� direction. This can be explained by the presence of the three 0� plies cluster in the center of the laminate.

Combined tension/Shear loading tests

The curves of the combined tensile and shear tests are presented in Fig. 16. The vertical arrows here also represent sudden, and loud, propagation of the damage, associated with a significant loss of stiffness. There is a rapid change of slope due to shear buckling despite the presence of a tension stress. The direction of the buckling is, however, deflected with respect to the shear (45�) with a slightly different direction (about 30�). The failure initiations (points R) are very similar to those observed previously except that the cracks at 45� on the surface are not observed here. The film of one of the tests is visible in Ref. [21,[START_REF]Vertex Tension + Shear Test[END_REF]. The final failure patterns for the three laminates are shown in Fig. 17. The patterns are quite different: for C3-1 and C3-2, the propagation is horizontal, as in pure tension and over practically the entire width of the plate. For C3-3, the failure pattern is completely asymmetrical: the propagation of the crack occurs almost horizontally to the left of the original notch, and is oriented at +45� on the other side. In addition, the combined stress strength of this laminate is much greater than for the other two.

Thanks to the VERTEX tests, it is possible to draw the response of the notched structures under under the three loading conditions. In Fig. 18, the three experimental points recorded for each configuration (laminate/loading) studied are represented: F corresponds to the local bending (or buckling), R to the first damage visually observed on the surface and the total breakage is defined when a first significant decrease in stiffness is observed. The points corresponding to the tensile stresses are the closest to the abscissa axis. Those associated with shear stresses are closest to the ordinate axis. The markers representing the tensile/shear stresses are located between the points mentioned above. The markers are connected by straight line segments to ease understanding. In the case of tensile and shear loading, the uniaxial character of the test degrades as a coupling mode exists in connection with the control displacement driving explained in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] and the local nonlinear phenomena (postbuckling and local bending). Overall, C3-1 and C3-2 show similar behavior. The crack propagation for C3-3 across the R-points and up to total failure is observed to be much more progressive than for the other two lay-ups. Thus the results on larger samples appear to confirm that the C3-3 laminate has the best resistance, as already noted on the coupon samples subjected to tension [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. In general, at this scale, the different interactions between the various phenomena make it difficult to give simple predictions of the quadratic criterion type.

Modeling of VERTEX tests

Discrete ply model strategy

The Discrete Ply Model (DPM) has already been described in several papers [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF][START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF][START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF][START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF][START_REF] Castanié | Theoretical and experimental analysis of asymmetric sandwich structures[END_REF][START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF] and the same characteristics were kept here as the ones used by Serra et al. for scaled open-hole tension test modeling [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF] and notched tensile tests [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF]. The main features are briefly recalled below. The principle of the DPM model is to use a dedicated mesh able to represent the discrete character of composite damage like matrix cracking, delamination and naturally, thanks to the geometry of the mesh their coupling (Fig. 19). This approach requires a complex mesh with coincident nodes at the interface (Fig. 19) and a special analysis to mesh holes or notches [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF][START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF]. The damage laws act as follows:

The delamination is taken into account using cohesive elements between two consecutive plies (or groups of plies) of different orientation, each ply being modeled with one volumic finite element through the thickness. Then the damage in the delamination interface elements is classically driven by cohesive law with the critical Strain Energy Release Rate (SERR).

The matrix cracking is taken into account using interface elements normal to the transversal direction. The damage of the matrix cracking interface elements is driven using quadratic Hashin fracture criteria evaluated in the neighboring volumetric finite elements. This interaction is known to be crucial to account for the complex damage morphology observed in composite structures.

Fiber fracture is taken into account using continuum damage mechanics with an original formulation linking the integration points of the volumetric finite element to impose a constant SERR per unit area [START_REF] Ostré | Edge impact modeling on stiffened composite structures[END_REF]. This approach can be compared to methods using the characteristic element length, which allow mesh-size independent modeling [START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF][START_REF] Marcin | Development of a macroscopic damage model for woven ceramic matrix composites[END_REF][START_REF] Laurin | A multiscale hybrid approach for damage and final failure predictions of composite structures[END_REF].

Another strong advantage of this approach is the limited number of parameters, all of which are provided by experimental tests [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF][START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF].

Test/computation dialog

Only the test specimens with double thickness plies C3-1d are modeled here. The general method of data exchange between a finite element model and the FE-SDIC measurements has been explained in the first part of this publication [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] and [START_REF] Pierré | Finite element stereo digital image correlation: framework and mechanical regularization[END_REF][START_REF] Pierré | Unstructured Finite Element-based Digital Image Correlation with enhanced management of quadrature and lens distortions[END_REF]. It gives access to an experimental displacement field directly expressed on the finite element mesh of the model. This avoids complex projections and interpolations between classic DIC and FEA tools. Doing so, experimental displacement-driven simulations are very easy to perform. In the case of the notched composite specimens, it is described in Fig. 20. The measurement is made on the red zone (350 350 mm 2) , which is smaller than the area of interest (400 400 mm 2 ) in order to limit edge effects due to bolting. FE-SDIC [START_REF] Pierré | Finite element stereo digital image correlation: framework and mechanical regularization[END_REF][START_REF] Pierré | Unstructured Finite Element-based Digital Image Correlation with enhanced management of quadrature and lens distortions[END_REF] measurements are provided in a 300 300 mm 2 area (where the material is assumed to behave elastically) minus a central rectangle around the notch (in which the evolution of the response can become strongly nonlinear). The boundary conditions are introduced in linear increments as in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF] and following the method described by Sztefek and Olsson [START_REF] Zagainov | Composite materials in aerospace design[END_REF]. To optimize the computation time, the Abaqus mesh uses three different zones (Fig. 21), generally following the same principle as in [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF].

The three components of the displacements measured by FE-SDIC method are first imposed on the edge of zone 1, 300 300 mm 2 (white square, Fig. 20). Then, in order to be able to simulate the rotation better, an out-of-plane displacement is imposed on 3 rows of additional nodes, (green area Fig. 20). In Zone 1 (Fig. 21), the size of the FE-SDIC and Abaqus meshes are identical. The boundary conditions are then simply prescribed dof by dof. The Abaqus mesh is made with SC8 R linear elements with one element in the thickness, size 10 10 3.25 mm 3 . In Zone 2, a mesh is made for each ply with C3D8R linear brick elements. The size of the elements are equal to twice the thickness, i.e. 0.5 0.5 0.25 mm 3 . In Zone 3, at the edge of the crack, the DPM method is applied with C3D8/COH3D8 damaging elements of the same size (0.5 0.5 0.25 mm 3 ). In contrast to the modeling of the notched specimens studied at the coupon scale [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF], and given the importance of bending/buckling phenomena and the absence of axial symmetry of the boundary conditions imposed, no symmetry is used in the case of the ''VERTEX" tests. Therefore, all the 13 plies are modeled. Since DPM imposes the density of the mesh by the thickness of the plies, the number of elements of Zones 2 and 3 is divided by 4 compared to the ''single thickness" layups. The calculation takes approximately 4 days (20 CPUs) for a total of 1.7 10 6 elements. For this reason, at this stage, the other tests were not modeled. As in the case of coupons [START_REF] Serra | Experimental and numerical analysis of CFRP notched coupons under tensile loading[END_REF], the phenomenon of ''shear locking" is not preponderant in this work because the elements located away from the cut use reduced integration (C3D8R and SC8R). However, this type of elements can often lead to poor modeling of bending phenomena because of hourglass modes. To remedy this problem, an artificial energy control of these modes is added to allow the most consistent representation possible. Since this artificial energy is negligible compared to the total internal energy (<1%), the quality of the results is not altered.

Results

Tensile loading test

The numerically determined stress/strain curve is very similar to the experimental one (Fig. 22). In order to maintain a reasonable calculation time, only 9 of the 44 DIC images used to determine the experimental curve are kept to determine the 8 steps of numerical simulation used. The overall behavior of the laminate appears to be correctly represented. The first surface-detected damage on the composite laminate (R), which is similar to that shown in Fig. 9, is well transcribed by the numerical model (Fig. 23), which correctly simulates matrix cracks oriented at 45� at the notch tip. A comparison with the observed displacement field (Fig. 24) shows that the DPM gives a very accurate simulation of the behavior observed experimentally, in particular the initiation of local bending around the notch. Fig. 25 compares the numerical and experimental failure patterns. On the right side of the specimen, the propagation of damage seems to occur initially in a direction oriented at 45� and then horizontally (orthogonal to the direction of tension). A similar scenario is obtained using the numerical model.

Shear loading test

A first computation was made without introducing any damage, in order to verify that the postbuckling response was correctly cap-tured, and then a second DPM computation was performed. Fig. 26 shows that the numerically determined stress/strain curve is very similar to the curve determined experimentally up to point F in the case of the shear test. Despite the fact that the experimental and numerical constraint curves are close, the numerical buckling is symmetrical to the real situation at the beginning of the computation and then shifts to the real solution after a certain time, with or without simulated damage. This shows that the method for introducing the displacements chosen can still be enhanced, especially through better consideration of initial imperfections and rotations, which is crucial for good modeling of postbuckling phenomena. Nevertheless, the DPM model correctly identifies the initiation of damage observed at point R (Fig. 27), deformed shapes at point F 0 (Fig. 28) and the failure pattern (Fig. 29). However, a divergence in the sense of a probable stiffness increase of numerical origin appears in the explicit DPM calculation beyond the point F 0 .

Envelope curves

The certification of composite structures is based on the pyramid of tests (see the first part of this paper [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project[END_REF]) and essentially on uniaxial tests brought to total failure. The general sizing of aeronautical composite structures in civil aircraft is complex and depends on the areas considered [START_REF] Niu | Composite airframe structures[END_REF][START_REF] Zagainov | Composite materials in aerospace design[END_REF][START_REF] Kassapouglou | Design and analysis of composite structures[END_REF][START_REF] Abrate | Dynamic Failure of Composite and sandwich structures[END_REF]. Generally speaking, the basic principles of sizing are postbuckling, damage tolerance, and repair and material issues: aging and material scatter (A and B values). In addition, the limit loads (LL) represent the maximum flying loads that an aircraft can handle (hard landing, avoidance maneuvers), for which the aircraft must be sized in fatigue, and ultimate loads (UL), the static limit loads multiplied by a factor of 1.5. The cumulative knockdown factors are represented in an empirical and qualitative way in Fig. 30. The blue curve is taken to represent the profile of the stresses during a flight. For a given problem, it can therefore be estimated that current composite structures will not exceed 33% of their average static failure value for limit loads and 50% for extreme loads. This shows both the margin of progression of composites in aeronautics and the robustness of current structures. For example, some aeronautical structures have demonstrated their ability, with impact and manufacturing defects in pristine condition, to withstand 3 times the lifetime in flight, and then withstand an extreme load to finally satisfy a crash test. The idea of envelope curves is that, in order to validate the sizing of an area or to answer a given problem (large notches, impact damage), it is possible to perform a load path using the VERTEX methodology at a percentage of the first failure load (R Points) for which the structure remains undamaged. In our case, the percentages used are 0.33, 0.5 and 0.66 and the results for the C3-1 laminate will be presented here. The loading path required by the test results is bilinear (Fig. 31). At this stage of the development of VER-TEX, the path will be determined from the displacements of the jacks obtained for the first failures.

The results of the envelope tests are presented in Fig. 32 in a strain space and Fig. 33 In Fig. 31, the segments drawn represent the load displacement paths of the actuators of the Vertex machine. These figures show the influence of the nature of the structural tests and of the ''blind" control by the displacement of the actuators at this stage of the use of the VERTEX methodology. The average stresses and strains are not pure because of secondary strains or stresses. A sudden unbuckling is at the origin of the strong discontinuities observed during the paths C ? D. It is important to note that a state of stress (or strain) is influenced by the ''path" followed to establish it. For the coefficient 0.33, the points A of the envelope curves are very close to the tensile curves while the points B and C are respectively distant from the curves of tension/shear and shear. These deviations between the points obtained by non-proportional loading and those obtained by proportional loading are accentuated when the coefficients 0.5 and 0.66 are used. It is clear that, because of the nonlinear nature of postbuckling, the loading path alters the behavior of the structure (buckling and unbuckling). The objective was nevertheless achieved, that is to say part of the Tension-Shear domain was swept in a single test on a notched specimen. X-ray tomography (Fig. 34) confirmed the absence of damage. The sample was shown to be healthy for all cases of extreme loads in the Tension-Shear range, and even beyond, as shown by the use of the factor 0.66.

Conclusions and perspectives

In these two papers, a methodology for testing composite structures under complex loadings has been developed. A specific test device for performing tensile/compression, shear and internal pressure tests has been designed, tested and validated. This type of test requires both control and monitoring of the resultant loadings applied to the specimen under test. Thus, a specific full-field measurement method based on a finite element approach to stereo image correlation (FE-SDIC) has been developed. The measurement mesh is shared with the finite element model, making data exchanges much easier. A specific protocol has been created and allows for effective and direct test/computation communications between simulated and experimental quantities. This methodology was then applied to the large cuts issue. The experiments revealed very complex behavior with coupling between postbuckling and crack propagation. DPM modeling was then used to successfully predict the first damage and crack propagation onset. The very long computation times (4 days) have allowed only two tests to be modeled for the moment, with plies 0.25 mm thick. Efforts therefore need to be made in terms of optimization of computation time and the damage model. A concept of envelope curves has also been presented. By means of a suitable load path, it enables validation of the structure integrity by means of a swept zone. At this stage, the boundaries of the domain are the experimental results obtained in tensile, shear and combined loadings. The next step will be to compute these boundaries numerically by the DPM method or by the VERTEX partners' approaches [START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF][START_REF] Marcin | Development of a macroscopic damage model for woven ceramic matrix composites[END_REF][START_REF] Laurin | A multiscale hybrid approach for damage and final failure predictions of composite structures[END_REF]. In the long term, certification of structures should evolve towards even more Predictive Virtual Testing (Fig. 35). The VERTEX approach and the principles of envelope test methodology open the way towards multiple perspectives leading to savings in cost and lead time in the certification process for aeronautical products. They contribute to perspectives of savings in testing, especially at expensive subcomponent and component levels, and to the de-risking at lower pyramid scale level, earlier in the process and at cheaper price, of design principles by intensive usage of numerical predictions supported by a smaller number of test validations. The contribution of VERTEX should ultimately lead to a significant reduction in the number of design loops and the cost of certification of aeronautical structures (Table 3).

Fig. 1 .

 1 Fig. 1. Typical ''Two-bay crack" test (Courtesy of Airbus).

Fig. 2 .

 2 Fig. 2. Vertex test Rig.

Fig. 3 .

 3 Fig. 3. View of a notched specimen bolted on the test rig.

Fig. 4 .

 4 Fig. 4. Combined loading actuators driving policy.

Fig. 5 .

 5 Fig. 5. Set of 6 cameras and areas covered on the specimen.

'

  'global" strains (e global xx and e global

  xy

Fig. 7 .

 7 Fig. 7. Notched tension, stress vs strain, averaged and normalized. Four specimens.

Fig. 8 .

 8 Fig. 8. Local bending (in red) around the notch, Point F C3-1 .

Fig. 9 .

 9 Fig. 9. View of a matrix crack at 45 �for the C3-1 Laminate at point R.

Fig. 10 .

 10 Fig. 10. Final failure patterns (a) C3-1, (b) C3-2, (c) C3-3.

Fig. 11 .

 11 Fig. 11. Comparison between coupon tests and VERTEX tests for notched specimens under tensile load.

Fig. 13 .

 13 Fig. 13. Out-of-plane displacement field for C3-1 shear test, at points F and R and just before failure.

Fig. 12 .

 12 Fig. 12. Notched shear, stress vs strain, averaged and normalized, four specimens.

Fig. 14 .

 14 Fig. 14. Cracks visible at the tip of the notch at point R and just before failure.

Fig. 15 .

 15 Fig. 15. Final failure patterns for shear tests.

Fig. 16 .

 16 Fig. 16. Notched tension + shear, shear stress vs tension stress, averaged and normalized, three specimens.

Fig. 18 .

 18 Fig. 18. Summary of the behavior of the three notched laminates under tension, shear and combined loading.

Fig. 19 .

 19 Fig. 19. Discrete ply model main principle and mesh coincidence at interface.

Fig. 20 .

 20 Fig. 20. Method of exchange of boundary conditions between FE-SDIC and numerical model.

Fig. 21 .

 21 Fig. 21. Abaqus Mesh for analysis of VERTEX tests.

Fig. 22 .

 22 Fig. 22. Experimental/numerical comparison for VERTEX tension test.

Fig. 23 .

 23 Fig. 23. DPM modeling of the first damage (point R).

Fig. 24 .

 24 Fig. 24. Comparison of the displacement fields resulting from measurement (DIC) and the numerical model between B and R (W) for maximum tension.

Fig. 25 .

 25 Fig. 25. Failure pattern: numerical vs experiment.

Fig. 26 .

 26 Fig. 26. Stress / strain curves for plate C3-1d, ''double thickness" under shear-Numerical / experimental comparison.

Fig. 27 .

 27 Fig. 27. First damage in the upper ply for double thickness shear test, comparison between experiment and DPM modelling.

Fig. 28 .

 28 Fig. 28. Displacements at point F 0 , shear test double thickness, comparison between DPM and test.

Fig. 29 .

 29 Fig. 29. Final failure for shear test, double thickness, comparison between experimental and numerical result.

Fig. 30 .

 30 Fig. 30. Qualitative overview of sizing policy for aeronautical composite structure.

Fig. 31 .

 31 Fig. 31. Loading paths for Envelope Curves.

  in a stress space. These two figures show the behavior of the envelopes performed for the different coefficients (0.33, 0.5 and 0.66) and those of the previous tests (Tensile, Shear, Tensile/Shear) for the C3-1 laminate. The continuous cycle corresponds to: O-A-B-C-D-A-B-C-D-A-B-C-D. The points associated with the first damage to the surface (R) visually observed are shown on the three test curves (Tensile, Shear, Tensile/Shear).

Fig. 33 .

 33 Fig. 33. Envelope tests in the stress space.

Fig. 32 .

 32 Fig. 32. Envelope tests in the strain space.

Fig. 34 .

 34 Fig. 34. Micro tomographic sections revealing no damage around the tip of the notch.

Fig. 35 .

 35 Fig. 35. Perspective for VERTEX: new pyramid of tests with more Predictive Virtual Testing.

Table 3

 3 Material properties of the T700/M21 used in the Discrete Ply Model. Tensile strain in fiber direction at damage initiation 1.70% e C 0 (%) Compressive strain in fiber direction at damage initiation Interface fracture toughness for shear mode (II & III) 1.6 N/mm

Table 1

 1 Tests performed.

Acknowledgments

The research that led to the results presented above received funds from the French National Research Agency under the VER-TEX project MATETPRO program (ANR -12 -RMNP-0001).The academic authors gratefully acknowledge the support provided by Airbus Group Innovation. The test rig was funded by Region Occitanie through a CPER 2010-2014 grant. The authors also wish to thank Sogeclair Aerospace for their work as the contractor on the design and manufacture of the test rig.