Ankit Chiplunkar
email: ankit.chiplunkar@airbus.com

Emmanuel Rachelson
email: emmanuel.rachelson@isae-supaero.fr

Michele Colombo
email: michele.colombo@airbus.com

Joseph Morlier
email: joseph.morlier@isae-supaero.fr

Approximate inference in related multi-output Gaussian Process Regression

Keywords: Gaussian Process, Kernel Methods, Approximate Inference, Multi-Output Regression, Flight-test data

In Gaussian Processes a multi-output kernel is a covariance function over correlated outputs. Using a prior known relation between outputs, joint auto-and cross-covariance functions can be constructed. Realizations from these joint-covariance functions give outputs that are consistent with the prior relation. One issue with gaussian process regression is efficient inference when scaling upto large datasets. In this paper we use approximate inference techniques upon multi-output kernels enforcing relationships between outputs. Results of the proposed methodology for theoretical data and real world applications are presented. The main contribution of this paper is the application and validation of our methodology on a dataset of real aircraft flight tests, while imposing knowledge of aircraft physics into the model.

Introduction

The main difference between the physical sciences and machine learning can be explained by the difference between deduction and induction. Physical sciences is deduction: where a very general formula is applied to a particular case. The basics of newtonian physics when applied to a particular aircraft geometry give us inertial loads. The basics of aerodynamics when applied to particular set of aircraft geometry and aircraft states give out aerodynamic pressures. Physical sciences take global rules and apply them to local configurations, whereas machine learning is induction [START_REF] Domingos | A few useful things to know about machine learning[END_REF]. It looks at local features and data, tries to find similarity measures between them and gives a global formula for the environment. For this reason machine learning algorithms perform better in presence of more and more data.

Unfortunately, gathering highly accurate data for physical systems is a costly exercise. Highly accurate CFD simulations may run for weeks and flight test campaigns cost millions of dollars. In this regard there is a dichotomy between these two fields of science, where machine learning needs more data for good performance but procuring data from physical systems is a costly exercise. In this work we consider using prior information of relationships between several outputs thereby effectively increasing the number of data-points.

We use multiple-output Gaussian Process (GP) regression [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] to encode the physical laws of the system and effectively increase the amount of training data points. Inference on multiple output data is also known as co-kriging [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF], multi-kriging [START_REF] Boyle | Dependent gaussian processes[END_REF] or Gradient Enhanced Kriging. Using a general framework [START_REF] Constantinescu | Physics-based covariance models for gaussian processes with multiple outputs[END_REF] to calculate covariance functions between multiple-outputs, we extend the framework of gradient enhanced kriging to integral enhanced kriging, quadratic enhanced kriging or any functional relationship between inputs.

Let us start by defining a P dimensional input space and a D dimensional output space. Such that {(x j i , y j i )} for j ∈ [1; n i ] are the training datasets for the i th output. Here n i is the number of measurement points for the i th output, while x j i ∈ R P and y j i ∈ R. We next define x i = {x 1 i ; x 2 i ; . . . ; x ni i } and y i = {y 1 i ; y 2 i ; . . . ; y ni i } as the full matrices containing all the training points for the i th output such that x i ∈ R ni×P and y i ∈ R ni . Henceforth we define the joint output vector Y = [y 1 ; y 2 ; y 3 ; . . . ; y D ] such that all the output values are stacked one after the other. Similarly, we define the joint input matrix as X = [x 1 , x 2 , x 3 , . . . , x D ]. If Σn i = N for i ∈ [1, D]. Hence N represents the total number of training points for all the outputs combined. Then Y ∈ R N and X ∈ R N ×P .

For simplicity take the case of an explicit relationship between two outputs y 1 and y 2 . Suppose we measure two outputs with some error, while the true physical process is defined by latent variables f 1 and f 2 . Then the relation between the output function, measurement error and true physical process can be written as follows.

y 1 = f 1 + n1 y 2 = f 2 + n2 (1)
Where, n1 and n2 are measurement error sampled from a white noise gaussian N (0, σ n1 ) and N (0, σ n2 ) respectively. While the physics based relation can be expressed as follows.

f 1 = g (f 2 , x 1 ) (2)
Here g is an operator defining the relation between f 1 and an independent latent output f 2 .

While a joint model developed using correlated covariance functions gives better predictions, it incurs a huge cost on memory occupied and computational time. The main contribution of this paper is to apply approximate inference on these models of large datasets and reduce the heavy computational costs incurred. For a multi-output GP as defined earlier the covariance matrix is of size N , needing O N 3 calculations for inference and O N 2 for storage. In this work we compare performance of variational inference [START_REF] Titsias | Variational learning of inducing variables in sparse gaussian processes[END_REF] and distributed GP [START_REF] Deisenroth | Distributed gaussian processes[END_REF] to approximate the inference in a joint-kernel.

The remaining paper proceeds as follows, section 2 provides the theoretical framework for Gaussian Process regression. Section 3 extends the multi-output GP regression in presence of correlated covariances. In section 4 various methods of approximating inference of a multi-output GP are derived. Finally, in section 5 we demonstrate the approach on both theoretical and flight-test data.

Gaussian Process Regression

A gaussian process is an infinite dimensional multi-variate gaussian. Such that any subset of the process is a multi-variate gaussian distribution. A gaussian process can be fully parametrized by a mean and covariance function Eq. 3.

y(x) = GP (m(x, θ), k(x, x , θ)) (3)
A random draw from a Gaussian Process gives us a random function around the mean function m(x, θ) and of the shape as defined by covariance function k(x, x , θ). Hence, Gaussian Process gives us a method to define a family of functions whose shape is defined by its covariance function. A popular choice of covariance function is a squared exponential function Eq. 4, because it defines a family of highly smooth (infinitely differentiable) non-linear functions as shown in Fig. 1(a).

k(x, x , θ) = θ 2 1 exp[- (x -x ) 2 2θ 2 2 ] (4)
A covariance function is fully parametrized by its hyperparameters θ i 's. For the case of Squared exponential kernel the hyperparameters are its amplitude θ 1 and length scale θ 2 .

Given a dataset (x, y) regression deals with finding latent function f between the inputs x and outputs y. While performing polynomial regression we assume that our function f comes from a family of polynomial functions. Since Gaussian Processes are such handy tools to define a family of non-linear functions. In a Gaussian Process Regression (GPR) we start with an initial family of functions defined by a GP called prior Fig. 1(a). Due to the bayesian setting of GP we can calculate the posterior mean and variance as shown in Eq. 6 and Eq. 7. This means that we are effectively eliminating all the functions in the prior that do not pass through our data points Fig. 1(c).

P(f | x, θ) = GP (y|0, K xx ) (5)
m(y * ) = k x * x (K xx ) -1 y (6) Cov(y * ) = k x * x * -k x * x (K xx ) -1 k xx * (7)
We can also improve our predictions by choosing a better prior. This involves optimizing the Marginal Likelihood (ML) P(y | x, θ) calculated as Eq. 8. The probability that our dataset (x, y) comes from a family of functions defined by the prior is called the ML [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. Hence, when we optimize the ML we are actually finding the optimal θ or family of functions that describe our data Fig. 1(b).

P(y | x, θ) = GP (y|0, K xx + σ 2 I) (8)
3 Multi-output Gaussian Process

Given a dataset for multiple outputs {(x i , y i )} for i ∈ [1; D] we define the joint output vector Y = [y 1 ; y 2 ; y 3 ; . . . ; y D ] such that all the output values are stacked one after the other. Similarly, we define the joint input matrix as X = [x 1 ; x 2 ; x 3 ; . . . ; x D ]. For the sake of simplicity, suppose we measure two outputs y 1 and y 2 with some error, while the true physical process is defined by latent variables f 1 and f 2 equation 2. The operator g(.) can be a known physical equation or a computer code between the outputs, it basically represents a transformation from one output to another.

Related Work

Earlier work developing such joint-covariance functions [START_REF] Bonilla | Multi-task gaussian process prediction[END_REF] have focused on building different outputs as a combination of a set of latent functions. GP priors are placed independently over all the latent functions thereby inducing a correlated covariance function. More recently it has been shown that convolution processes [START_REF] Alvarez | Latent force models[END_REF] , [START_REF] Boyle | Dependent gaussian processes[END_REF] can be used to develop joint-covariance functions for differential equations. In a convolution process framework output functions are generated by convolving several latent functions with a smoothing kernel function. In the current paper we assume one output function to be independent and evaluate the remaining auto-and cross-covariance functions exactly if the physical relation between them is linear [START_REF] Solak | Derivative observations in gaussian process models of dynamic systems[END_REF] or use approximate joint-covariance for non-linear physics-based relationships between the outputs [START_REF] Constantinescu | Physics-based covariance models for gaussian processes with multiple outputs[END_REF].

Multi-output Joint-covariance kernels

A GP prior in such a setting with 2 output variables is expressed in equation 9.

f 1 f 2 ∼ GP 0 0 , K 11 K 12 K 21 K 22 (9)
K 12 and K 21 are cross-covariances between the two inputs x 1 and x 2 . K 22 is the auto-covariance function of independent output, while K 11 is the autocovariance of the dependent output variable. The full covariance matrix K XX is also called the joint-covariance. While, the joint error matrix will be denoted by Σ;

Σ = σ 2 n1 0 0 σ 2 n2 (10)
Where, n1 and n2 are measurement error sampled from a white noise gaussian N (0, σ n1 ) and N (0, σ n2 ).

For a linear operator g(.) the joint-covariance matrix can be derived analytically [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF], due to the affine property of Gaussian's equation 11.

f 1 f 2 ∼ GP 0 0 , g(g(K 22 , x 2 ), x 1 ) g(K 22 , x 1 ) g(K 22 , x 2 ) K 22 (11)
Using the known relation between outputs we have successfully correlated two GP priors from equation 2. This effectively means that when we randomly draw a function f 2 it will result in a correlated draw of f 1 such that the two draws satisfy the equation 2. We have effectively represented the covariance function K 11 in terms of the hyperparameters of covariance function K 22 using the known relation between outputs.

Without loss of generality we can assume that the independent output f 2 belongs to a family of functions defined by a Squared Exponential kernel. The joint-covariance between f 1 and f 2 , means that a random draw of independent function f 2 will result in a correlated draw of the function f 1 . The fig: 2(a) shows random draws coming from a differential relationship between f 1 (red) and f 2 (blue) such that f 1 = ∂f2 ∂x . We can see that the top figure is derivative of the bottom one since f derivative goes to zero where f independent goes to maxima or minima. Similarly, the fig: 2(b) shows random draws coming from an integral relationship between f 1 (red) and f 2 (blue) such that f 1 = f 2 . We can see that the top figure is integral of the bottom one since f independent goes to zero where f integral goes to maxima or minima.

GP Regression Using Joint-Covariance

We start with defining a zero-mean prior for our observations and make predictions for y 1 (x * ) = y * 1 and y 2 (x * ) = y * 2 . The corresponding prior according to equation 9 and 10 will be:

Y (X)) Y (X * )) = GP 0 0 , K XX + Σ K XX * K X * X K X * X * + Σ (12)
The posterior distribution is then given as a normal distribution with expectation and covariance matrix given by [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]

m(y * ) = K X * X (K XX ) -1 Y (13) Cov(y * ) = K X * X * -K X * X (K XX ) -1 K XX * (14)
Here, the elements K XX , K X * X and K X * X * are block covariances derived from equations 11. Due to the bayesian setting we have basically eliminated all the functions that do not pass through the points defined by the observed data.

The joint-covariance matrix depends on several hyperparameters θ. They define a basic shape of the GP prior. To end up with good predictions it is important to start with a good GP prior. We minimize the negative log-marginal likelihood to find a set of good hyperparameters. This leads to an optimization problem where the objective function is given by equation 15

log(P(y | X, θ)) = log[GP (Y |0, K XX + Σ)] (15)
With its gradient given by equation 16

∂ ∂θ log(P(y | X, θ)) = 1 2 Y T K -1 XX ∂K XX ∂θ K -1 XX Y - 1 2 tr(K -1 XX ∂K XX ∂θ ) (16)
(a) Joint draws between f1 red and f2 blue such that f1 = ∂f 2 ∂x . The posterior is conditioned such that f1 = 5; x = 0 and f2 = 0; x = 0. The draws from the posterior follow the relationship between outputs and also the conditioning.

(b) Joint draws between f1 red and f2 blue such that f1 = f2. The posterior is conditioned such that f1 = 0; x = 0 and f1 = 0; x = 1. The draws from the posterior follow the relationship between outputs and also the conditioning.

Fig. 3. Multi-Output Gaussian Process Regression Predictions

Here the hyperparameters of the prior are θ = l 2 , σ 2 2 , σ 2 n1 , σ 2 n2 . These correspond to the hyperparameters of the independent covariance function K 22 and errors in the measurements σ 2 n1 and σ 2 n2 . Calculating the negative log-marginal likelihood involves inverting the matrix K XX + Σ. The size of the K XX + Σ matrix depends on total number of input points N , hence inverting the matrix becomes intractable for large number of input points.

The fig: 3(a) shows mean and variance for a differential relationship between independent function f 2 (blue) and differential function f 1 (red) and such that f 1 = ∂f2 ∂x . We have conditioned the two functions such that f 1 = 5; x = 0 and f 2 = 0; x = 0. This means that the function f 2 passes through 0 and has a derivative equal to 5 at x = 0. We have not maximized the marginal likelihood for this case and the hyperparameters of K 22 are θ 1 = 1; θ 2 = 0.2. All the corresponding draws from the posterior GP also follow the conditioning.

The fig: 3(b) shows mean and variance for a integral relationship between independent function f 2 (blue) and differential function f 1 (red) and such that f 1 = f 2 .dx. We have conditioned the two functions such that f 1 = 0; x = 1 and f 1 = 0; x = 0. This means that the function f 2 has an integral 0 at the two points [0, 1] . We have not maximized the marginal likelihood for this case and the hyperparameters of K 22 are θ 1 = 1; θ 2 = 0.2. All the corresponding draws from the posterior GP also follow the conditioning. We can observe that the draws will also have an integral 0 between the range [0, 1].

In the next section we describe how to solve the problem of inverting huge K XX + Σ matrices using approximate inference techniques.

Approximating Inference

The above GP approach is intractable for large datasets. For a multi-output GP as defined in section 3.2 the covariance matrix is of size N , where O N 3 time is needed for inference and O N 2 memory for storage. Thus, we need to consider approximate methods in order to deal with large datasets.

Inverting the covariance matrix takes considerable amount of time and memory during the process. Hence, almost all techniques to approximate inference try and approximate the inversion of covariance matrix K XX . If a covariance matrix is diagonal or block-diagonal in nature then methods such as mixture of experts are used eg. distributed GP. Whereas if the covariance matrix is more spread out and has similar terms in its cross diagonals then low-rank approximations are used eg. variational approximation. The remaining section details the two methods for approximating covariance matrix which can be later used to resolve equations 13, 14 and 16.

Variational Approximation on Multi-output GP

Sparse methods use a small set of m function points as support or inducing variables. Suppose we use m inducing variables to construct our sparse GP. The inducing variables are the latent function values evaluated at inputs x M . Learning x M and the hyperparameters θ is the problem we need to solve in order to obtain a sparse GP method. An approximation to the true log marginal likelihood in equation 15 can allow us to infer these quantities.

We try to approximate the joint-posterior distribution p(X|Y ) by introducing a variational distribution q(X). In the case of varying number of inputs for different outputs, we place the inducing points over the input space and extend the derivation of [START_REF] Titsias | Variational learning of inducing variables in sparse gaussian processes[END_REF] to multi-output case.

q(X) = N (X|µ, A) (17)
Here µ and A are parameters of the variational distribution. We follow the derivation provided in [START_REF] Titsias | Variational learning of inducing variables in sparse gaussian processes[END_REF] and obtain the lower bound of true marginal likelihood.

F V = log(N [Y |0, σ 2 I + Q XX ]) - 1 2σ 2 T r( K) (18)
where

Q XX = K XX M K -1 X M X M K X M X and K = K XX -K XX M K -1 X M X M K X M X .
K XX is the joint-covariance matrix derived using equation 11 using the input vector X defined in section 1. K X M X M is the joint-covariance function on the inducing points X M , such that

X M = [x M 1 , x M 2 , ..., x M 2 ]
. We assume that the inducing points x M i will be same for all the outputs, hence x M 1 = x M 2 = ... = x M 2 = x M . While K XX M is the cross-covariance matrix between X and X M .

Note that this bound consists of two parts. The first part is the log of a GP prior with the only difference that now the covariance matrix has a lower rank of M D. This form allows the inversion of the covariance matrix to take place in O N (M D) 2 time. The second part as discussed above can be seen as a penalization term that regularizes the estimation of the parameters.

The bound can be maximized with respect to all parameters of the covariance function; both model hyperparameters and variational parameters. The optimization parameters are the inducing inputs x M , the hyperparameters θ of the independent covariance matrix K 22 and the error while measuring the outputs σ. There is a trade-off between quality of the estimate and amount of time taken for the estimation process. On the one hand the number of inducing points determine the value of optimized negative log-marginal likelihood and hence the quality of the estimate. While, on the other hand there is a computational load of O N (M D) 2 for inference. We increase the number of inducing points until the difference between two successive likelihoods is below a predefined quantity.

Distributed Inference on Multi-output GP

An alternative to sparse approximations is to learn local experts on subset of data. Traditionally, each subset of data learns a different model from another, this is done to increase the expressiveness in the model [START_REF] Rasmussen | Infinite mixtures of gaussian process experts[END_REF]. The final predictions are then made by combining the predictions of local experts [START_REF] Chen | Bagging for gaussian process regression[END_REF].

An alternative way is to tie all the different experts using one single set of hyperparameters [START_REF] Deisenroth | Distributed gaussian processes[END_REF]. This is equivalent to assuming one single GP on the whole dataset such that there is no correlation across experts as seen in figure 5(b). This tying of experts acts as a regularization and inhibits overfitting. Although ignoring correlation among experts is a strong assumption, it can be justified if the experts are chosen randomly and with enough overlap.

If we partition the dataset into M subsets such as The above equation 19 describes the formulation for marginal likelihood. Due to the independence assumption the marginal likelihood can be written as a sum of individual likelihoods and then can be optimized to find the bestfit hyperparameters. After learning the hyperparameters we can combine the predictions of local experts to give mean and variance predictions. The robust Bayesian Commitee Machine (rBCM) model combines the various experts using their confidence on the prediction point [START_REF] Deisenroth | Distributed gaussian processes[END_REF]. In such manner experts which have high confidence at the prediction points get more weight when compared to experts with low confidence.

D (i) = X (i) , y (i) , i = 1, . . . M . log p(y|X, θ) ≈ M k=1 log p k (y (i) |X (i) , θ) (19)
m(Y * ) = (Cov(X * )) -2 β k σ -2 k m k (X * ) (20) (Cov ( (Y * )) -2 = k β k σ -2 k + (1 - k β k )σ -2 * * (21)
In the above equations m k (X * ) and σ k are the mean and covariance predictions from expert k at point X * . σ * * is the auto-covariance of the prior at prediction points X * . β k determines the influence of experts on the final predictions [START_REF] Cao | Generalized product of experts for automatic and principled fusion of gaussian process predictions[END_REF] and is given as

β k = 1 2 (log σ -2 * * -log σ -2 k ).

Experiments

We empirically assess the performance of distributed Gaussian Process and Variational Inference with respect to the training time and accuracy. We start with a synthetic dataset where we try to learn the model over derivative relationship and compare the two inference techniques. We then evaluate the improvement on real world flight-test dataset.

The basic toolbox used for this paper is GPML provided with [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], we generate covariance functions to handle relationships as described in equations 11 using the "Symbolic Math Toolbox" in MATLAB 2014b. Variational inference is wrapped from gpStuff toolbox [START_REF] Vanhatalo | Gpstuff: Bayesian modeling with gaussian processes[END_REF] and distributed GP is inspired from [START_REF] Deisenroth | Distributed gaussian processes[END_REF]. All experiments were performed on an Intel quad-core processor with 4Gb RAM.

Experiments on Theoretical Data

We consider a derivative relationship between two output functions as described in equation 2. Such that

g(f, x) = ∂f ∂x
Since the differential relationship g(.) is linear in nature we use the equation 11 to calculate the auto-and cross-covariance functions as shown in table 1.

Table 1. Auto-and cross-covariance functions for a differential relationship.

Initial Covariance K22 σ 2 exp( -1 2 d 2 l 2 ) Cross-Covariance K12 σ 2 d l 2 exp( -1 2 d 2 l 2 ) Auto-covariance K11 σ 2 d 2 -l 2 l 4 exp( -1 2 d 2 l 2 )
Data is generated from equations 22, a random function is drawn from GP to get f 2 whose derivative is then calculated to generate f 1 . y 1 and y 2 are then calculated by adding noise according the the equations 22. 10,000 points are generated for both the outputs y 1 and y 2 . Values of y 2 are masked in the region x ∈ [0, 0.3] the remaining points now constitute our training dataset.

f 2 ∼ GP [0, K SE (0.1, 1)] σ n2 ∼ N [0, 0.2] σ n1 ∼ N [0, 2] (22)
K SE (0.1, 1) means squared exponential kernel with length scale 0.1 and variance as 1.

Figure 6 shows comparison between an independent fit GP and a joint multioutput GP whose outputs are related through a derivative relationship described in equation 5.1. For figure 6(a) using variational inference algorithm we optimize the lower bound of log-marginal likelihood, for independent GP's on y 1 and y 2 . For figure 6(b) using variational inference we optimize the same lower bound but with a joint-covariance approach as described in section 4.1 using y 1 , y 2 and g(.). We settled on using 100 equidistant inducing points for this exercise [START_REF] Chiplunkar | Sparse physics-based gaussian process for multi-output regression using variational inference[END_REF] and have only optimized the hyperparameters to learn the model. Experiment was run on 10,000 points but only 100 data points are plotted to increase readibility. Inference is performed using variational inference algorithm and equidistant 100 inducing points. We can observe the huge difference between the real data and the predicted mean values at zone with no data.

(b) Joint multi-output GP's for two outputs, mean is represented by solid black line. 2Σ confidence band is represented by light red for f1 and light blue for f2.The dashed black line represents the true latent function values; noisy data is denoted by blue dots. Experiment was run on 10,000 points but only 100 data points are plotted to increase readibility. Inference is performed using variational inference algorithm and equidistant 100 inducing points. We can observe the improved prediction between zone with no data because information is being shared between the two outputs. For the second experiment we compare the Root Mean Squared Error (RMSE) and run-times of distributed GP and Variational Inference algorithms while performing approximate inference. We progressively generate from 10e3 to 10e5 data-points according to the equations 22 and 5.1. We separated 75% of the data as training set and 25% of the data as the test set, the training and test sets were chosen randomly. The variational inference relationship kernel as described in section 4.1 with 100 equidistant inducing points was used. We learn the optimal values of hyper-parameters for all the sets of training data. The distributed GP algorithm as described in section 4.2 was used with randomly chosen 512 points per expert. We learn the optimal values of hyper-parameters for all the sets of training data. The accuracy is plotted as RMSE values with respect to the test set. The runtime is defined as time taken to calculate negative In figure 7(a) the time to calculate negative log-marginal likelihood with increasing number of training points is calculated. As expected the full GP takes more time when compared to variational inference or distributed GP algorithms. The Variational inference algorithm has better run-time till 10e4 data-points after that distributed GP takes lesser time. In figure 7(b) the RMSE error with test set is compared between the variational inference and distributed GP algorithm. Here too the variational inference algorithm performs better lesser number of datapoints but distributed GP starts performing better when we reach more than 10e4 datapoints.

One thing to note is that we have fixed the number and position of inducing points while performing this experiment. While a more optimized set of inducing points will have better results, for datasets of the order 10e5 distributed GP algorithm starts outperforming variational inference. Moreover, upon observing figure 5(b) we can say that the covariance matrix in a joint GP setting is not diagonal in nature and hence an approximation technique which can compensate between low-rank approximation and diagonal approximation should be investigated [START_REF] March | Askit: Approximate skeletonization kernelindependent treecode in high dimensions[END_REF].

Experiments on Flight Test Data

We perform experiments on flight loads data produced during flight test phase at Airbus. Loads are measured across the wing span using strain gauges. Shear load T z and bending moment M x as described in figure 8 are used as two outputs for this exercise. η or point of action of forces and angle of attack α are the two inputs. The aircraft is in quasi-equillibrium in all conditions and there are no dynamic effects observed throughout this dataset. All data is normalized according to airbus policy. The relation between T z and M x can be written as:

M x (η, α) = η edge η T Z (x, α)(x -η)dx (23)
.

The equation 23 is applicable for the η axis. Here, η edge denotes the edge of the wing span. The forces are measured at 5 points on the wing span and at 8800 points on the α dimension. We compare plots of relationship-imposed multioutput GP and independent GP. Then we compare the measures of negative-log marginal likelihood and RMSE for varying number of inducing points.

Figure 9(a) shows the independent (blue shade) and joint fit (red shade) of two GP. The top figure shows T Z with the variance of dependent GP plotted in red and variance of independent GP plotted in blue. Bottom figure shows plots for M X . Since the number of input points is less than 10e5 we use variational inference. 100 inducing points in the input space are used to learn and plot the figure. The variance of red is smaller than that of blue showing the improvement in confidence when imposing relationships in the GP architecture. The relationship between T Z and M X gives rise to better confidence during the loads prediction. This added information is very useful when identifying faulty sensor data since equation 23 will push data points which do not satisfy the relationship out of the tight confidence interval.

Figure 9(b) shows improvement in the negative log-marginal likelihood and RMSE plots upon increasing number of inducing points. 10 sets of experiments were run on 75% of the data as training set and 25% of the data as the test set, the training and test sets were chosen randomly. We learn the optimal values of hyper-parameters and inducing points for all the 10 sets of experiments

Conclusions and Future Work

This paper presents approximate inference methods for physics-based multiple output GP's. We extend the variational inference and distributed GP inference techniques to be applied on physics-based multi-output GP's and reduce the computational load for inference.

Section 5.1 demonstrates the advantage of using multi-output GP in presence of prior known relationships on a theoretical dataset. Significant improvement in prediction can be observed in presence of missing data due to transfer of information occuring due to the prior relationship. Then, we compare the difference in accuracy and run times upon using distributed GP or variational inference to approximate inference. Although variational inference performs better for datasets of size less than 10e4, distributed GP becomes significantly advantageous for datasets greater than 10e5.

Section 5.2 shows real world application of joint kernel on flight loads data. We demonstrate that adding prior physics-based relationships allows us to create a robust, physically consistent and interpretable surrogate model for these loads. Aircraft flight domain consists of various maneuvers, this is further complicated by several relationships between outputs. In the future we wish to develop better strategies to exploit the clustered nature of flight-domain and more than two relationships.

Fig. 1 .

 1 Fig. 1. Gaussian Process Regression

  (a) Joint draws between f1 red and f2 blue such that f1 = ∂f 2 ∂x .(b) Joint draws between f1 red and f2 blue such that f1 = f2.

Fig. 2 .

 2 Fig. 2. Multi-Output Gaussian Process Random Draws

Fig. 4 .

 4 Fig. 4. Variational approximation of covariance matrix for Gaussian Process Regression

( a )Fig. 5 .

 a5 Fig. 5. Distributed approximation of covariance matrix for Gaussian Process Regression

  (a) Independent fit for two GP's mean is represented by solid black line. 2Σ confidence band is represented by light red for f1 and light blue for f2.The dashed black line represents the true latent function values; noisy data is denoted by blue dots.

Fig. 6 .

 6 Fig.6. Experimental results for differential relationship while using variational approximation

  (a) Comparison of time to calculate negative log marginal likelihood for a Full GP, Variational inference and Distributed GP with increasing number of datapoints. We observe for datapoints greater that 10e5 the distributed GP algorithm starts outperforming variational inference (b) Comparison of RMSE for variational inference and distributed GP algorithm. We observe for datapoints greater that 10e4 the distributed GP algorithm starts outperforming variational inference.

Fig. 7 .

 7 Fig. 7. Comparison of run time and RMSE between distributed GP and Variational Inference

Fig. 8 .

 8 Fig. 8. Wing Load Diagram

  (a) 2σ confidence interval and mean of the dependent GP are represented in red shade and solid red line. 2σ confidence interval and mean of the independent GP are represented in blue shade and solid blue line. Experiment was run on 8800 data points Noisy data is denoted by circles only 1 α step is plotted. Confidence interval improves upon adding the relationship kernel. (b) Progression of RMSE and loglikelihood upon increasing number of inducing points. Top plot shows the value of mean and variance of negative logmarginal likelihood. The bottom figure in blue shows the mean and variance of root mean squared error. 10 sets of experiments were run on 75% of the data as training set and 25% of the data as the test set, the training and test sets were chosen randomly.

Fig. 9 .

 9 Fig. 9. Experimental results for aircraft flight loads


Laurent Mezeix

Simon Dols

Christophe Bouvet

Bruno Castanie
email: bruno.castanie@insa-toulouse.fr

Jean-Paul Giavarini

Natthawat Hongkarnjanakul

Experimental analysis of impact and post-impact behaviour of inserts in Carbon sandwich structures

Keywords: Carbon Laminates, impact, insert, pull-through, sandwich

In aeronautics, honeycomb sandwich structures are widely used for secondary structures such as landing gear doors, flaps or floors, and for primary structures in helicopters or business jets. These structures are generally joined by using local reinforcements of the insert type. In the present study, 50 J low velocity impact tests were performed on inserts using a drop-weight device and the impact response and failure patterns were analysed. Impacted specimens were then pull-through tested to failure. Some of the tests were stopped before final failure in order to obtain precise details on the failure scenario. It was shown that, in the cases studied, the residual strength after impact was very high (about 90%) in comparison to the large reductions habitually observed in compression after impact tests.

Introduction

Inserts in sandwich structures are widely used in many fields, such as aeronautics, space and marine or civil engineering. The literature contains a very large number of studies concerning impacts on sandwich structures and their strength after impact (e.g. see literature [START_REF] Castanie´b | Modelling of low energy/low velocity impact on nomex honeycomb sandwich structures with metallic skins[END_REF][START_REF] Castanie´b | Core crush criteria to determine the strength of sandwich composite structures subjected to compression after impact[END_REF][START_REF] Hayman | Approaches to damage assessment and damage tolerance for FRP sandwich structures[END_REF][START_REF] Raju | Impact damage resistance and tolerance of honeycomb core sandwich panels[END_REF][START_REF] Anderson | Experimental investigation of low-velocity impact characteristics of sandwich composites[END_REF]). Sizing methodologies for impact have also been proposed [START_REF] Zenkert | Damage tolerance assessment of composite sandwich panels with localised damage[END_REF] and then evaluated [START_REF] Mcquigg | Compression after impact experiments and analysis on honeycomb core sandwich panels with thin facesheets[END_REF] but the specific problem of impact on junctions has not been addressed. Neither is this issue treated in reference books on impact in sandwich structures [START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Abrate | Dynamic failure of composite and sandwich structures[END_REF]. In general, studies on inserts focus on static sizing methods based on analytical models that are more or less relevant; details can be found in the literature [START_REF] Zenkert | Sandwich construction[END_REF][11][START_REF]Insert design handbook[END_REF]. These methods are effective in practice as long as appropriate allowables are taken. Thomsen et al. [START_REF] Thomsen | Sandwich plates with 'through-the-thickness' and 'fully-potted' insert: evolution of differences in structural performance[END_REF][START_REF] Thomsen | Analysis and design of sandwich plates with inserts -a highorder sandwich plate theory approach[END_REF][START_REF] Bozhevolnaya | Local effects across core junctions in sandwich panels[END_REF] have developed refined analytical models that represent the stress concentrations well at the skin/insert/ core triple point and propose design rules for fully potted or partially potted inserts. An interesting application of these theories is found in the horizontal gradients of foam core that have been proposed to reduce stress concentration [START_REF] Bozhevolnaya | Structurally graded core inserts in sandwich panels[END_REF]. Nevertheless, these methods are limited to the linear domain. Bunyawanichakul et al. [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF] have shown that the failure scenario of a fully potted insert under pullthrough is complex: core shear buckling and tears in the vicinity of the insert, crushing of the potted insert under the screw head and then punching of the laminate skin. These authors propose a nonlinear finite element model to capture this scenario. Other research efforts have focused on various themes such as the strength variability of an insert in relation with the manufacturing process [START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF], the reliability of inserts [START_REF] Smith | Reliability of inserts in sandwich composite panels[END_REF], the development of new technologies [START_REF] Woo | Development of the hybrid insert for composite sandwich satellite structures[END_REF], the development of automated methods of installation [START_REF]New Insert allows automatic placement in composite sandwich panels[END_REF] and the behaviour of inserts under environmental effects [START_REF] Bin | Failure characteristics of carbon/BMI-Nomex sandwich joints in various hygrothermal conditions[END_REF]. Enhanced FE Software and finite element formulation have allowed some authors to propose very refined models. Heimbs and Pein [START_REF] Heimbs | Failure Behaviour of honeycomb sandwich corner joints and inserts[END_REF] have partially modelled insert tearing by representing the honeycomb core with its actual geometry and capturing the buckling of the cells in the vicinity of the insert. Roy et al. [START_REF] Roy | Testing and modeling of nomex TM honeycomb sandwich panels with bolt insert[END_REF] used a similar strategy to model the pull-through of metallic inserts. Their proposed model simulates the linear load/displacement portion of the insert pull-out properly and closely predicts the onset of stiffness reduction. An experimental study was conducted by Song et al. [START_REF] Song | An experimental study of the insert joint strength of composite sandwich structures[END_REF] to investigate the effect of design variables, such as the core height and density, the face thickness and the insert clearance, on the failure loads of a sandwich insert joint. However, despite the obvious interest of the topic, the authors have found no papers on the subject of impact and insert (except their own conference papers [START_REF] Dols | Post-impact behavior of insert in sandwich structures[END_REF][START_REF] Dols | Impact et tenue re´siduelle des inserts de structures sandwichs[END_REF]) and the subject is worthy of investigation.

So, the aim of this paper is first to analyse the damage due to an impact on an insert. There is a wide range of inserts and the choice made here was to use the same technology as Bunyawanichakul et al. [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF]. Impact tests were a priori performed on the insert and some distance away from it. The failure patterns were analysed. Then pull-through tests were conducted in the same way as in Bunyawanichakul et al. [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF] or Adam et al. [START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF]. Stopped tests were also used in order to propose a failure scenario. Finally, some conclusions on modelling strategies are proposed.

Materials and methods

The core was made of Nomex TM honeycomb with a cell size of 4.8 mm, a density of 48 kg/m 3 and a thickness of 10 mm. The potted inserts were made with 3 M Scotch-Weld TM 3500-2B/A. The Nomex cells were filled with an in-house tool as shown in Figure 1. The diameter of the insert was about 30 mm. The shape obtained was almost the same as in Raghu and Battley [START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF] and dispersion could not be avoided. The insert was cured at 175 C for 1 h. The skins were manufactured and cured separately. They were made of unidirectional carbon ply impregnated with dry resin and a hardener mixture (Araldite LY 5052 and hardener HY 5052) baked for 4 h at 80 C under a pressure of 7 bar. Two different skins were made: a thick skin (lower skin [45, 90, À45, 0] 2s , 3.2 mm thick) and a thinner skin (top skin [45, 90, À45.0] s , 1.6 mm thick). This choice was related to an innovative aeronautic floor designed by the industrial partners. The core and the skins were bonded with a redux film under vacuum at 120 C for 1 h to form 140 Â 140 mm 2 sandwiches after trimming. Finally, a hole of diameter 6.35 mm was drilled in the centre of each test piece (Figure 2). Fourteen specimens were manufactured (12 for impact tests and two pristine). Impact tests were performed using a drop tower system (same as in Castanieé t al. [START_REF] Castanie´b | Modelling of low energy/low velocity impact on nomex honeycomb sandwich structures with metallic skins[END_REF]) with a 4 kg spherical impactor 16 mm in diameter. They followed the Airbus Industries Test Method (AITM 1-0010) (Figure 3(a)), except for the size of the specimen (140 Â 140 mm 2 ) which was different because of the post-impact tests. The square specimens simply sat on a 120 Â 120 mm 2 frame. The velocity of the impactor was measured with an optical laser just before its impact with the specimen. A piezoelectric force sensor was placed inside the impactor to measure the acceleration during contact, which allowed the contact force during impact to be computed. The displacements of the rear face were measured by a second laser sensor. The anti-rebound system was manual, economic and efficient; it comprised a simple plank of wood driven by a quick young human neuronal system (Figure 3(b)). The specimen was then positioned so that, a priori, impact occurred on the insert or the honeycomb alone. Several tests at different impact energies were performed and an energy of 50 J was selected so that a residual dent >0.3 mm Barely Visible Impact damage (BVID) remained visible.

The same device as used by Bunyawanichakul et al. [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF] was employed to perform pull-through tests with a 100 kN Instron testing machine (Figure 4). The specimen was simply supported on a square frame with an internal window of 110 Â 110 mm 2 . A machined bolt (NAS 6704-91) was mounted on the specimen as shown in Figures 5 and6. Then the nut was hand tightened and thus the insert was held between the bolt head and the spacer. The grips of the traction machine were applied to the machined part of the screw. A displacement sensor measured the displacement of the screw head during the test and the tensile force was given by the load cell of the machine. The test consisted of pulling the screw until it punched through the lower skin, leading to collapse of the structure. The test was displacement controlled at a speed of 2 mm/min.

Analysis of impact tests Discussion on impact location

Tests showed that there were three impact positions that generated quite different responses and damage scenarios: on the insert (case 1, Figure 7), on the border of the insert (case 2, Figure 7) and on the core (case 3, Figure 7). The geometric dispersion due to the traditional method of filling cells increased the experimental dispersion of this phenomenon, as already reported by Raghu and Battley [START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF]. The force/displacement curves were actually quite different (Figure 7) but the dispersion between the impact locations was low, allowing the impact site to be found a posteriori and the cases to be ranked. The peak force was higher in the case of an impact on the insert, which was to be expected given the greater stiffness of the structure. Moreover, all curves had the same, globally sinusoidal, beginning corresponding to a linear response of the structure, followed by clear force drops indicating that the cells of the honeycomb were first damaged by shear buckling of the cell edges in the vicinity of the potted insert [START_REF] Castanie´b | Modelling of low energy/low velocity impact on nomex honeycomb sandwich structures with metallic skins[END_REF][START_REF] Castanie´b | Core crush criteria to determine the strength of sandwich composite structures subjected to compression after impact[END_REF][START_REF] Abrate | Dynamic failure of composite and sandwich structures[END_REF]. This is consistent with what has already been reported in the literature. On the other hand, a very significant difference was observed between the impacts on the insert and on the core (Figure 8) for the displacements of the rear face, indicating that a very different failure scenario occurred. It was also noted that the peak force was about 20% higher when the impact was on the insert.

Damage analysis for impact on the insert (case 1)

The tests were sorted and collected by impact position category. For three impact tests on the insert, the force/displacement curves are shown in Figure 9. The dispersion is low, which shows that the failure scenario is probably still the same for a given configuration. The shape of these curves is very similar to those observed for impact on laminate [START_REF] Hongkarnjanakul | Validation of low velocity impact modelling on different stacking sequences of CFRP laminates and influence of fibre failure[END_REF], suggesting that the honeycomb plays only a small role. To confirm this analysis, C-scans were carried out and vertical sections were cut (Figures 10 to 12). It appears that the indentation was very small (Figure 10(a)) when the impact was on the insert. This is why the relatively high impact energy of 50 J was used. The damage in the upper skin was located to the right of the impactor. Some damage appeared in the honeycomb but was limited to two or three circumferential cells around the insert. The skin was crushed locally at the impact but the damage was contained. Under the impact location, the insert seemed intact but the lower skin had significant damage. The C-scan (Figure 11) showed very large delaminations oriented on the side of the impact location. These were also visible in the micrographic section (Figure 12). It seems, therefore, that the insert retransmitted most of the energy of the impact to the lower skin. This confirms (i) that the predominant damage occurred in the lower laminated skin and (ii) the a priori analysis of the shape of the force/displacement curve.

Damage analysis for impact on the border of the insert (case 2)

The load/displacement curves for three samples are given in Figure 13. Qualitatively, the shape of the curve is quite different with a fairly large plateau that is typical of the response of honeycomb core sandwich to impact [START_REF] Castanie´b | Modelling of low energy/low velocity impact on nomex honeycomb sandwich structures with metallic skins[END_REF][START_REF] Castanie´b | Core crush criteria to determine the strength of sandwich composite structures subjected to compression after impact[END_REF][START_REF] Hayman | Approaches to damage assessment and damage tolerance for FRP sandwich structures[END_REF][START_REF] Raju | Impact damage resistance and tolerance of honeycomb core sandwich panels[END_REF][START_REF] Anderson | Experimental investigation of low-velocity impact characteristics of sandwich composites[END_REF][START_REF] Zenkert | Damage tolerance assessment of composite sandwich panels with localised damage[END_REF][START_REF] Mcquigg | Compression after impact experiments and analysis on honeycomb core sandwich panels with thin facesheets[END_REF][START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Abrate | Dynamic failure of composite and sandwich structures[END_REF]. There is no peak force but the plateau force is about 20% less than the maximum force for an impact on the insert. The same type of analysis was conducted with postmortem microscopic observation and C-scans (Figures 14 and15). The residual indentation was more marked, with the appearance of broken fibres on the surface of the upper skin. Damage was also different in the core and in the insert. The insert was crushed on half its thickness and the honeycomb core was also crushed locally [START_REF] Aminanda | Experimental analysis and modeling of the crushing of honeycomb cores[END_REF]. A C-scan of the lower skin showed large delaminations similar to that observed in the case of a direct impact on the insert. So in this case, even though the insert was observed to have undergone localized crushing, it remained sufficiently rigid during an impact to transmit the forces to the lower skin and to create large delamination. There was also marked damage to the upper skin (fibre breakage) but it remained localized around the impact. Qualitatively, the energy dissipated (area inside the force/displacement curve) was greater in this case than in the other two configurations, probably because the effects of several damage modes were accumulated.

Damage analysis for impact on the core (case 3)

The load/displacement curves for three samples are given in Figure 16. Here, we find the classic response of sandwich structures with Nomex honeycomb core under impact [START_REF] Castanie´b | Modelling of low energy/low velocity impact on nomex honeycomb sandwich structures with metallic skins[END_REF][START_REF] Castanie´b | Core crush criteria to determine the strength of sandwich composite structures subjected to compression after impact[END_REF][START_REF] Hayman | Approaches to damage assessment and damage tolerance for FRP sandwich structures[END_REF][START_REF] Raju | Impact damage resistance and tolerance of honeycomb core sandwich panels[END_REF][START_REF] Anderson | Experimental investigation of low-velocity impact characteristics of sandwich composites[END_REF][START_REF] Zenkert | Damage tolerance assessment of composite sandwich panels with localised damage[END_REF][START_REF] Mcquigg | Compression after impact experiments and analysis on honeycomb core sandwich panels with thin facesheets[END_REF][START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Abrate | Dynamic failure of composite and sandwich structures[END_REF][START_REF] Zenkert | Sandwich construction[END_REF]. In this case, the plateau area is smaller than in the previous case because there is energy absorption in the core, which explains the decrease of force at the end of the test and the absence of rebound. Because of the high impact energy (50 J), the upper skin was perforated (Figure 17) and the honeycomb significantly crushed [START_REF] Aminanda | Experimental analysis and modeling of the crushing of honeycomb cores[END_REF]. However, as shown in the C-scan of Figure 18, no delamination damage was apparent in the lower skin, which confirms the force transmission role of the insert during the impact and the absorption role of Nomex honeycomb which protects the lower skin.

Summary

A summary of the damage in the three cases is proposed in Table 1. The damage patterns differ significantly from one case to another. Especially when the impact occurs on the insert or even on its edge, the stiffness is sufficient to transfer impact force to the lower skin. So large delamination occurs in this lower skin although only a small part of the upper skin is damaged at the right of the impact. The delamination patterns are original and only one side of the specimen is affected, the rest being somehow ''protected'' by the insert. This is very different when the impact occurs on the honeycomb. In this case, the skin is perforated locally and the honeycomb absorbs a large part of the energy. So the lower skin is protected and no damage is found.

Pull-through after impact tests Test on pristine specimen

It was important to perform a test that would serve as a reference for explaining the failure scenarios and differences between the impact cases. The test curve is given in Figure 19. The curve is partially similar to Bunyawanichakul's results [START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF][START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF]. Part 1 is the linear response without damage. Part 2 is limited by a first change in slope and a slight load jump followed by a recovery of stiffness. To check the damage mode in this phase, a test was stopped and a C-scan taken. It showed that the skins were not damaged (Figure 19). A section also showed no damage in the core, which is consistent with the scenario of post-buckling shear in cells of the honeycomb in the vicinity of the insert presented in the literature [START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF][START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF]. Damage in Zone 3 was not clearly identified as it was not possible to do a stopped test because of a lack of pristine specimens (only 2). The behaviour diverged from the literature [START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF][START_REF] Raghu | Strength variability of inserts in sandwich panels[END_REF]. An explanation will be given by analogy with test impacted specimens later. The shape of the plateau zone (4) clearly showed that the screw punched the lower skin.

Test on impacted specimens

The force/displacement curves for the four configurations are shown in Figure 20. The curves are very similar despite quite different damage patterns. If we look at higher loads, or Zone 4, it is clear that the reduction related to the impacts is weak, especially when compared to compression after impact tests [START_REF] Castanie´b | Core crush criteria to determine the strength of sandwich composite structures subjected to compression after impact[END_REF][START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Abrate | Dynamic failure of composite and sandwich structures[END_REF]. The knowdown is less than 5% in case of impact on the insert and remains below 15% in case of impact on the edge or on the core. In the case of impact on the insert, the main type of damage is delamination in the lower skin. These delaminations have very little influence on the sandwich bending stiffness. The final failure is driven by the punching of the lower skin [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF]. So the impact damage modes on the insert do not have any great influence on the final failure. In the case of impact on the core, the damaged upper skin area is limited and, in the same way, has little influence on the final rupture scenario even though the reduction is more significant here. However, in the case of the pristine sample, the maximum displacement at failure is greater and the punch most progressive. It is possible that the punching of the lower skin interacts with a local buckling of the lower skin occurring due to the large delamination and the compression caused by overall bending. The other possibility is that the upper side impact creates local splitting that generates a sharp break. A post-mortem pattern with such splitting is shown in Figure 21(b) and one of these local buckling events is shown in Figure 21(a). If we look at part 1 of the test curves, we see that the slope of the three damaged inserts is similar while the slope of the pristine sample is higher. Then in part 2, the slopes become identical. As the C-scan of the stopped test (Figure 19) did not show any delamination in the lower skin, the first damage of the pristine sample can be attributed to post-buckling of the cells of the honeycomb in the vicinity of the insert. This scenario has already been demonstrated by Bunyawanichakul et al. [START_REF] Bunyawanichakul | Experimental and numerical analysis of inserts in sandwich structures[END_REF][START_REF] Bunyawanichakul | Nonlinear finite element analysis of inserts in sandwich structures[END_REF] both numerically and experimentally. Since the slope of the pristine sample is identical to those of the impacted specimens in part 2, it can be deduced that the honeycomb cells of these specimens are also sufficiently damaged during the impact to lose their initial stiffness. This observation is consistent with the qualitative observations of Figures 10, 14 and 17. The behaviour in zone 3 is more difficult to explain. The four types of specimens have the same slope, which would mean a priori that the damage is identical and that the slight jump up from the areas 2 and 3 of the pristine specimen corresponds to the creation of identical damage. Nevertheless, it has been shown that the damages of the three impacted specimens are quite differents. Only one stopped test was conducted in this area, on the specimen impacted on the honeycomb, and cracks were found near the drilled hole in the lower skin (Figure 22). This shows that the damage occurring during the tests was probably more complex and coupled.

Conclusion

An experimental study was conducted on impacts on an insert in sandwich structures. The behaviour and damage varied according to the position of impact: on the insert, on the edge of the insert or on the core. In the first two cases, the insert transmitted the impact force to the lower skin of the sandwich, creating large delaminations in it. In the case of an impact on the core, classical damage was retrieved and the lower skin and the insert were little affected. Pull-through tests were carried out and, whatever the impact location, the strength after impact remained high, losses of only 5-15% being observed. The full damage scenario is difficult to identify because many modes interact: buckling and fracture of the honeycomb cells, punch of the skin, delaminations, local buckling, fibre breakage near the hole edge. Only a dialogue between testing and calculation and more stopped tests would clarify the failure scenario. The modelling is challenging and only very advanced models, e.g. the discrete ply model used by the authors on similar problems [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], would capture damage in the skin. A fully nonlinear model able to capture compression/shear buckling of Nomex honeycomb, like the one recently developed by Seemann and Krause [START_REF] Seemann | Numerical modeling of nomex honeycomb sandwich cores at meso-scale level[END_REF], is also necessary.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors wish to thank Sogeclair Aerospace for their financial and operational support.

Figure 1 .

 1 Figure 1. Manufacturing of the potted insert.

Figure 2 .

 2 Figure 2. Design and photo of finished sample.

Figure 4 .

 4 Figure 4. Overview of pull-through tests.

Figure 3 .

 3 Figure 3. Impact device (a) details and (b) overview.

Figure 6 .

 6 Figure 6. Specimen ready for pull-through tests.

Figure 5 .

 5 Figure 5. Machined bolt for loading the insert.

Figure 8 .

 8 Figure 8. Rear face displacement versus time for impact on the insert and on the core.

Figure 7 .

 7 Figure 7. Influence of impact location on the impact response of the insert.

Figure 9 .

 9 Figure 9. Force/displacement responses for impact on the insert.

Figure 10 .

 10 Figure 10. (a) Typical dent after impact on the insert and (b) post-mortem observation of a specimen.

Figure 11 .

 11 Figure 11. C-scan of the lower skin after an impact on the potted insert.

Figure 13 .

 13 Figure 13. Force/displacement responses for impact on the border of the insert.

Figure 12 .

 12 Figure 12. Micrographic section of lower skin after an impact on the potted insert.

Figure 14 .

 14 Figure 14. (a) Typical dent after impact on the border of the insert and (b) post-mortem observation of a specimen.

Figure 15 .

 15 Figure 15. C-scan of the lower skin after an impact on the border of the insert.

Figure 16 .

 16 Figure 16. Force/displacement responses for impact on the honeycomb core.

Figure 17 .

 17 Figure 17. (a) Typical dent after impact and (b) post-mortem observation of a specimen impacted on the core.

Figure 18 .

 18 Figure 18. C-scan of the lower skin after an impact on the honeycomb core.

Figure 19 .

 19 Figure 19. Pull-through test on a pristine specimen.

Figure 20 .

 20 Figure 20. Comparison of pull-through response for the different cases.

Figure 21 .

 21 Figure 21. (a) Local buckling of the lower skin after pull-through, post-mortem view, case of impact on insert. (b) Failure pattern of the upper skin with splitting of the upper ply near the impact.

Figure 22 .

 22 Figure 22. Local crack around the hole drilled in the lower skin, impact on the honeycomb core specimen.

Table 1 .

 1 Summary of damages occurring in the impact on a potted insert.

		Damages
	Impact
	location	Upper skin	Core	Lower skin
	Potting	Slight indentation	Cracking 45 of the honeycomb,	Large delamination
			local crushing
	Border	Indentation,	Cracking 45 of the honeycomb,	Large delamination
		slight punching	crushing
	Honeycomb	punching	Crushing of the honeycomb	No delamination