Bijan Mohammadi
email: bijan.mohammadi@um2.fr

Keywords: Uncertainty quantification, non-intrusive, principal angles between subspaces, worst-case analysis, Value at Risk, parallelism, robust design

We propose a systematic procedure for both aleatory and epistemic uncertainty quantification of numerical simulations through geometric characteristics of global sensitivity spaces. Two mathematical concepts are used to characterize the geometry of these spaces and to identify possible impacts of variability in data or changes in the models or solution procedures: the dimension of the maximal free generator subspace in vector spaces and the principal angles between subspaces. We show how these characters can be used as indications on the aleatory and epistemic uncertainties. In the case of large dimensional parameter spaces, these characterizations are established for quantile-based extreme scenarios and a multi-point moment-based sensitivity direction permits to propose a directional uncertainty quantification concept for directional extreme scenarios (DES). The approach is non-intrusive and exploits in parallel the elements of existing mono-point gradient-based design platforms. The ingredients of the paper are illustrated on a model problem with the Burgers equation with control and on a constrained aerodynamic performance analysis problem.

Introduction

The classification of uncertainties in epistemic (or reducible) and aleatory (or stochastic and non reducible) categories is now well established [START_REF]AIAA Guide for the Verification and Validation of Computational Fluid Dynamics Simulations[END_REF][START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF][START_REF] Iaccarino | Quantification of Uncertainty in Flow Simulations Using Probabilistic Methods[END_REF]. These uncertainties have different origins. In particular, the former can be reduced improving our models and solution procedures. Several Uncertainty Quantification approaches exist. They belong to two main categories. Without being exhaustive, one can use either statistical techniques directly through, for instance, Monte Carlo simulations, knowing the Probability Density Function (PDF) of our uncertain parameters (together with variance reduction techniques), or use different intelligent sampling (e.g. stratified, latin hypercube). Or, one can try to recast a new mathematical model which is then solved using deterministic methods. Methods such as perturbation expansion methods for random fields or stochastic operator expansions belongs to this category. Non-sampling expansions, such as the polynomial chaos methods [START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF], also belong to this class. In all cases, one can also introduce low-order models (e.g. response surface methods) [START_REF] Schilders | Model order reduction: Theory, research aspects and applications[END_REF][START_REF] Obinata | Model reduction for control system design[END_REF][START_REF] Qu | Model Order Reduction Techniques with Applications in Finite Element Analysis[END_REF] to a priori assimilate available information and use these reduced order models in the methods above for uncertainty quantification instead of the initial high-fidelity models.

This work proposes a systematic way to quantify epistemic and aleatory uncertainties through geometric characterizations of global sensitivity spaces. We consider two very generic situations. In the first case we would like to analyze the sensitivity of state variables with respect to a given functioning parameter (e.g. the flow Mach number in an aerodynamic problem). In the second case we are interested in a situation where the simulation aims at predicting a given quantity of interest (e.g. the maximum value of a variable in a given area) and there are several control parameters involved in addition to the functioning parameter mentioned before. The approach for these uncertainty quantification is through the analysis of geometric characteristics of two global sensitivity spaces defined by the gradients of the quantities of interest for a sampling of the functioning parameter range. The sensitivity spaces are built using adequate gradients of either the state variables or a functional with respect to different functioning parameters.

Once these global sensitivity spaces built, we analyze the dimensions of their respective free generator subspaces. Previous works have shown how to use these dimensions for adaptive sampling [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Gallard | Robust parametric design by multipoint optimization[END_REF]. In the presence of different modeling or solution methods, principal angles between these sensitivity subspaces permit to measure the deviation due a change in the modeling. Indeed, different complexity can be envisaged for the governing equations and the equations can be solved with different numerical schemes. Each situation provides a different set of global sensitivity spaces. The dimensions of the corresponding sensitivity spaces are therefore interesting measures for both the epistemic and aleatory (at given modeling procedure) uncertainties. Indeed, if at given modeling procedure, the dimensions of the sensitivity spaces remain unchanged when enriching the sampling of the functioning parameter range, this would be a first indication of the level of sensitivity of the simulations with respect to this parameter. Once this is established, principal angles between subspaces permit to analyze both the impact of a given evolution of the modeling on the sensitivity spaces or an enrichment of our sampling. Eventually, constant dimension and low angles will clearly indicate a situation of low uncertainty.

These ingredients can be used in a context of multi-point robust analysis of a system to define worst-case scenarios for its functioning. To this end we combine a multi-point search direction with the probabilistic features of the control parameters through their quantiles and, in particular, their Value at Risk [START_REF] Ph | Value at Risk: The New Benchmark for Managing Financial Risk[END_REF][START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF]. These ingredients permit to define the concept of directional uncertainty quantification and directional extreme scenarios (DES). Global sensitivity spaces are then built for these extreme scenarios and the above geometric characteristics permit again to measure the impact of the variability of the parameters of the problem. The main interest of this construction is to account for the variability of the parameters in large dimension without a sampling of the corresponding parameter space.

Concerning the cost of these analysis, one can say that, when using the same calculation ingredients that for a high-fidelity simulation (i.e. without calling for low-order models or cheaper discretizations), the best calculation complexity one might think of for a simulation under uncertainty is when its cost is comparable to the deterministic situation. This is clearly unreachable. Now, suppose the proposed approach for uncertainty quantification is such that all the extra effort can be achieved in a fully parallel manner and parallel to the initial deterministic calculation. This permits for the time to solution to remain unchanged for a simulation when accounting for the presence of uncertainties. This is the case with the sampling based approaches mentioned above, but those are quite expensive in large dimensions and do not take advantage of all the available simulation environments. In particular, when an adjoint-based optimization environment exists. This work proposes an original way to upgrade existing platform without abandoning what has been built for the deterministic situations and with keeping the time to solution unchanged in the presence of uncertainties.

The ingredients of the paper are illustrated on a model problem with a Burgers equation with a distributed control and on a three dimensional flow problem in a context of robust shape sensitivity analysis and design involving a complex computational environment both in terms of state equation and geometry manipulation.

Geometric quantification of variability in global search spaces

Let us denote by α ∈ I ⊂ IR a scalar control parameter (we call this a functioning parameter such as the Mach or the Reynolds numbers) and by x ∈ IR n the remaining set of control parameters (e.g. parameters in a shape parameterization model, etc). Together, x and α describe all the functioning conditions of our system. The size n of x is usually large. We also consider a functional j(U (α, x)) and a state variable U (α, x) ∈ IR N with N >> n. U is solution of a state equation F (U (α, x)) = 0, also called the governing equations. The solution of F is usually expensive.

The first situation we consider is when we would like to quantify the uncertainties on the state U with respect to α given through its interval of variation. The work takes place in an environment where the sensitivities ∇ α U (α, x) ∈ IR N of the state with respect to the functioning parameter α are available. Because α is scalar, these sensitivity evaluations are easy and can be carried out by finite differences or using the complex variable method noticing that:

U (α + iε, x) = U (α, x) + iε∇ α U (α, x) + o(ε 2 ),
which permits to access simultaneously both the state and its sensitivity:

U (α, x) = Re(U (α + iε, x)), ∇ α U (α, x) = ε -1 Im(U (α + iε, x)).
The complex variable method avoids the difficulty in finite differences with the choice of ε but requires the source code to have its variables in complex instead of real. Now, consider a multi-point evaluation of U (α + iε, x) for a sampling of size m of I denoted by I m ⊂ I.

A first vector space of global sensitivity can be defined as:

s m = Span(∇ α U (α k , x), α k=1,...,m ∈ I m ) ⊂ IR N ×m .
The geometry of this space indicates the global sensitivity of the state U with respect to α.

The second situation we consider is when the simulations aim at predicting some quantities of importance. These can be, for instance, extreme values reached in some area (e.g. maximum temperature on a blade or a measure of performance). These can be seen as the evaluation of a scalar functional j.

This time the work takes place in an environment where the sensitivities ∇ x j(α, x) ∈ IR n of the functional with respect to the control parameters are available. This is, for instance, the case if one uses a gradient based optimization platform with, ideally, an adjoint formulation to make the cost of the evaluation of ∇ x j(α, x) independent of the size n of the control space. We give an example of this in section 4.

Again, we consider a multi-point evaluation of the functional j and its sensitivities for the sampling I m . And, a second vector space of global sensitivity can be defined as:

S m = Span(∇ x j(α k , x), α k=1,...,m ∈ I m ) ⊂ IR n×m .
To simplify the notations, we have considered a same sampling for the two situations. But this is not a prerequisite. Two quantities of interest are the dimensions p and q of the subspaces s m and S m .

In particular, a system can be said robust if the free subspaces generating s m and S m have dimensions p, q << m.

These dimensions depend on the sampling I m of the range of variation of α. They also depend on the modeling (the governing equations) and the way the equations are solved (the numerical schemes). Each situation provides a different set global sensitivity spaces s m and S m . p and q are therefore interesting primary quantificators for both the aleatory and epistemic uncertainties.

2.1. Algorithm for the construction of the global sensitivity spaces. The calculations necessary for the construction of s m and S m can be carried out in a fully parallel manner. Suppose we have a direct simulation chain linking the parameters (α, x) to the state U solution of a state equation F (U (α, x)) = 0 and to a functional j. The following parallel algorithm describes the different steps one needs to take to build the different ingredients we need for our analysis.

                       Given x, I m ,                    -m parallel state evaluations U (α l , x), α l ∈ I m ,
-m parallel state sensitivity with respect to α: ∇ α U (α l , x), α l ∈ I m , -m parallel evaluations of j(α l , x), α l ∈ I m , -m parallel solutions of the adjoint state V equation:

V t F U (U (α l , x)) = j t U , α l ∈ I m , -m parallel evaluations of ∇ x j(α l , x) = j x + (V t F x ) t , α l ∈ I m .
As previously mentioned, α being scalar, ∇ α U (α l , x) is built with finite differences or with the complex variable method. Once the gradients available we can evaluate the actual dimensions of the two vector spaces s m and S m . This analysis is carried out with an incomplete Gram-Schmidt orthonormalization algorithm.

2.2.

Incomplete Gram-Schmidt orthonormalization. Dimensions p and q can be estimated by Gram-Schmidt orthonormalization. Having a set of vectors {v 1 , ..., v m }, the dimension d = dim(Span{v 1 , ..., v m }) is the rank at which the rest u l after successive projections is below a threshold value T OL:

(1)

                   {v 1 , ..., v m }, T OL = given, u 1 = v 1 , iterations l = 1, ..., m u l = v l -l-1 j=1 <v l ,u j > <u j ,u j > u j , until u k < T OL, d = l.
The algorithm above gives p (resp. q) with v

l = ∇ α U (α l , x) (resp. v l = ∇ x j(α l , x)).
The choice of an approximate or incomplete Gram-Schmidt procedure is deliberate and is made to avoid pollution by numerical noise or artifacts [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. Indeed, a full orthonormalization might bring extra dimensions not containing any pertinent information. In the sequel we give several examples of such over-estimation of the free subspace generator and also discuss this issue in section 2.3.1 in a context of multi-point robust optimization.

But, the threshold T OL can also be seen as a multi-scale cut-off where small scales in the fields falling below it are not taken into account in the definition of the global sensitivity spaces. Finally, this construction also gives a level of confidence on the sampling I m of I: if a new direction built for a new α is not linear combination of previous directions, this new sample point should be accounted for in I m+1 . Gram-Schmidt orthonormalization has been previously used in uncertainty quantification [START_REF] Witteveen | Efficient quantification of the effect of uncertainties in advectiondiffusion problems using polynomial chaos[END_REF].

2.3. Link with multi-point minimization. The spaces s m and S m represent the global sensitivity of the problem with respect to α. In the context of multi-point minimization this α-dependency is accounted for in the definition of the descent direction using a vector of S m (e.g. through a linear combination of ∇ x j(α k , x) involving weights). Let us propose a suitable choice for this construction which will be used in section 2.4 to account for directional variability in x using the value at risk of these control parameters.

Consider the following constrained minimization of a functional µ(x) encapsulating this dependency using a weighting which can be non uniform and target-based [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF]:

(2) min

x µ, such that σ * (x) ≤ σ * 0 , where

µ = 1 ω α l ∈Im ω l j(α l , x), ω = m l=1 ω l , (3)
the weights ω l account for the kind of performance we eventually want for the design: constant performance over the functioning parameters ranges. For the sake of simplicity and without loss of generality one can consider here both the sampling and the weighting to be uniform. As in First-Order Second Moment (FOSM) methods [START_REF] Melchers | Structural Reliability Analysis and Prediction[END_REF], σ * is to monitor the regularity of the final performance which should be as regular as possible with respect of α:

(4) σ * = 1 2 α l ∈Im < ∇ α j(α l , x), ∇ α j(α l , x) >,
where <, > is the Euclidean scalar product. σ * 0 is the level of variability for the initial design or of a reference configuration. Here also, as the size of α is usually small (here one), ∇ α j is evaluated by finite differences on I m . Problem (2) can be solved using a descent method with the descent direction d given by:

(5) d = ∇ x µ-< ∇ x µ, ∇ x σ * > ∇ x σ * + η∇ x σ * ,
where 0 < η 1 and a = a/ a is the normalized vector a and ( 6) < ∇ α j(α l , x), ∇ αx j(α l , x) > .

∇ x µ = 1 ω α l ∈Im ω l ∇ x j(α l , x),
∇ xα j(α l , x) is derived from ∇ x j(α l , x) by finite differences on I m , component by component. The definition of the descent direction permits to make sure that both µ and σ * decrease for small descent steps. Indeed, a first order development in x gives:

σ * (x + δx) -σ * (x) = ∇ x σ * ∇ x σ * .(δx) = -ρη ∇ x σ * ≤ 0,
and we have

µ(x + δx) -µ(x) = ∇ x µ ∇ x µ.(δx) = -ρ ∇ x µ (1 -ζ 2 + ηζ),
where ζ =< ∇ x µ, ∇ x σ * >. Therefore, µ is also decreasing as 1 -ζ 2 + ηζ ≥ 0 because ζ is a cosine function and therefore |ζ| ≤ 1 and η is chosen such that 0 < η 1 as shown in figure 1.

2.3.1.

Reducing the size of the search space. The descent direction d permits a control of both moments µ and σ * but to make the minimization efficient one should not consider all elements ∇ x j(α l , x) at α l ∈ I m , but only pertinent ones. These are the q vectors (q < m << n) in the maximal free generator in S m at given T OL [START_REF] Gallard | Robust parametric design by multipoint optimization[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. This gain in efficiency is related to the fact that working with only q vectors instead of m drastically reduces the dimension of the search space, making the underlying optimization problem simpler. Figure 2 illustrates this idea in the context of alternate minimization by successive projections when all the m directions are considered (split in two sets of q and m -q directions) versus a situation where the minimization takes place only on the subspace of dimension q. This obviously affects the convergence in a gradient based minimization.

(m-q) directions X 0 Figure 2. Impact of a reduction in the dimension of the search space by considering the 'pertinent' q directions generating S m at given T OL.

2.4.

Value at Risk for worst-case sensitivity spaces. We showed how to introduce some geometric characterizations of the global sensitivity spaces built with the gradients of the state ∇ α U (α l , x) or a functional ∇ x j(α l , x) for α l ∈ I m a sampling of a functioning parameter variation range. This was feasible as there are few of such parameters in a simulation procedure (usually one or two of importance). Now, we would like to account for the variability in the control parameter x ∈ IR n with 1 n. Therefore, a sampling of the variability ranges of the components of x is out of question because computationally too expensive.

In [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF] we show how to introduce the Value at Risk (VaR) concept in optimization algorithms with the aim of quantifying our confidence on the optimal solution at low complexity without a sampling of the control space. We would like to use this concept here and define, in combination with the direction d (5), directional uncertainty quantification analysis for extreme scenarios.

In financial engineering, the Value at Risk (VaR) is a widely used a-quantile to measure the risk of loss on a given asset [START_REF] Ph | Value at Risk: The New Benchmark for Managing Financial Risk[END_REF]. It defines, for a given probability level (0 < a < 1) and time horizon (typically one day), a threshold value for the loss X on the asset:

VaR a = inf{l ∈ IR : P (X > l) ≤ 1 -a}. A given control parameter x i=1,...,n represents, with a confidence level of a, the interval [x i + VaR - a , x i + VaR + a ] around x i and VaR - a ≤ 0 ≤ VaR + a . One expresses then possible deviations from x via two probability density functions (PDF) and the corresponding values at risk. One can make the hypothesis that the upper and lower bounds of the variations are symmetric, in which case VaR - a = -VaR + a . This is typically the case when uncertainties on a parameter follow a Gaussian distribution. Gaussian distributions are interesting as their Values at Risk are explicitly known: VaR 0.99 = 2.33 and VaR 0.95 = 1.65 for N (0, 1) and VaR α (N (0, σ)) = σVaR α (N (0, 1)).

To summarize x + VaR a (x) is a closed domain in IR n around x:

B a (x) = {z ∈ IR n : x i + VaR - a (x i ) ≤ z i ≤ x i + VaR + a (x i ), i = 1, .
., n}. Now, we would like to use this information together with the descent direction d given by ( 5) to define extreme values for the local variability of the functional along the direction d. We make the assumption that the functional is monotonic over B a (x), continuous and bounded. It is therefore uniformly continuous. These are

X d X+VaR a +/-along d B a (x) Figure 3. Sketch of x ± = d ∩ ∂B a (x).
reasonable requirements for a robust design to be possible. The intersection of the boundary of B a (x) and d permits to define two extreme values for x as worst-case scenarios:

(8)

x ± = d ∩ ∂B a (x).
A sketch of this construction is given in figure 3. The assumptions above on the functional make that its extreme values over B a (x) are reached on the boundary ∂B a (x). This requirement fails near critical points of the functional. In [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF] we showed how this can be linked with the notion of over-solving where it becomes useless to solve accurately near an optimum when the variation in a control parameter by a descent method falls below the level of the variability for this parameter. In other words, in the presence of uncertainties, one should not distinguish between the points in B a (x). The next step is to define the global sensitivity spaces s m and S m at points x ± following the procedures described in the previous sections using ∇ α u(α l , x ± ) and ∇ x j(α l , x ± ) for α l ∈ I m . Let us call these s ± m and S ± m . Their respective dimensions, denoted by p ± and q ± , can be found at given T OL by incomplete Gram-Schmidt orthonormalization of section 2.2. These give very precious information on the impact of local variability on x on our global search spaces. And from the calculation complexity issue, it is important that these have been built without a sampling of our parameter space of dimension n for x.

2.5. Angles between subspaces . We discussed how to evaluate the dimension of different sensitivity spaces in the presence of uncertainties. Another interesting information comes with their respective positions or orientations. A natural tool to quantity these is given by the principal angles between vector subspaces. We use the mathematical concept of 'principal angles' between subspaces in the Euclidean space IR n initially introduced by C. Jordan [START_REF] Jordan | Essay on geometry in n dimensions[END_REF]. If the maximum principle angle between the two subspaces is small, the two are nearly linearly dependent. Geometrically, this is the angle between two hyperplanes embedded in a higher dimensional space.

Let us briefly recall the concept of principal angles and how to practically compute them [START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF][START_REF] Jiang | Angles between Euclidean subspaces[END_REF]. We present the approach for subspaces of the same dimension k, but it is not necessary for the two subspaces to be of the same size in order to find the angles between them. We need n ≥ 2k to be able to exhibit two orthogonal subspaces. If n < 2k, some principal angles necessary vanish and for n = k they all vanish.

The k principal angles θ i , i = 1, ..k are between 0 and π/2. This is an important point and will be used later to take advantage of the positivity of the cosine of the angles. The principal angles are defined as:

cos(θ i ) = < a i , b i > a i b i = max{ < a, b > a b : a ⊥ a m , b ⊥ b m ; m = 1, .., i -1},
where a j ∈ A and b j ∈ B. The procedure finds unit vectors a 1 ∈ A and b 1 ∈ B minimizing the angle θ 1 between them. It then takes the orthogonal complement of a 1 in A and b 1 in B and iterates. This procedure is not useful in practice as computationally inadequate. We would like to be able to find the angles θ i from the inner products < a i , b j > of the elements of two basis of A and B [12]1 Suppose A and B are two subspaces of dimension k of IR n , n ≥ 2k and let {a i , i = 1, .., k} and {b i , i = 1, .., k} be two arbitrary orthonormal basis for A and B. Orthonormal basis are easy to obtain through the Gram-Schmidt orthonormalization procedure. Consider M being the matrix of the projection operator P r A of B on A defined by:

P r A (b i ) = k j=1 < b i , a j > a j , M = (< b i , a j >) i,j .
The principal angles can be linked to this operator [START_REF] Shonkwiler | Poincare duality angles for Riemannian manifolds with boundary[END_REF] through:

M = GΣH t ,
where G and H are orthogonal matrices and Σ = diag(cos(θ i )).

As G and H are orthogonal matrices, this is a Singular Vector Decomposition (SVD) of M . G and H are unknown at this point. But, we will show that we do not need them to get the θ i . Otherwise, the approach would be be again computationally useless.

We recall that the columns of G are the left-singular vectors of M and eigenvectors of M M t and the columns of H are the right-singular vectors of M and eigenvectors of M t M . And most important that cos 2 (θ i ) are the eigenvalues of P r t A P r A which writes in matrix form as:

M t M = (GΣH t ) t (GΣH t ) = HΣ 2 H t with Σ 2 = diag(cos 2 (θ i )).
Therefore, to find the principal angles between subspaces A and B, knowing an orthonormal basis in each subspace, one should calculate M and find the eigenvalues of M t M and take the square root of them. This last operation is valid as the angles are between 0 and π/2, and their cosine therefore always positive. Finding the eigenvalues of M is easy because the size of M is small. This procedure is still valid if the subspaces have different dimensions. The projection operator can be defined as well as its transpose and the eigenvalues of M t M are real as this is a symmetric square matrix. In our case, n being the size of the control parameter x, it is usually much larger than k which is the dimension of the generator subspace of S m at given T OL.

2.6. Definition of subspaces A and B in our multi-point formulation. Let us define subspaces A and B in the context of our multi-point formulation. In our problem the size n is very large and is given by the dimension of the control space. Subspaces A and B are of dimension k and we have 2k << n. For the sake of simplicity we consider that both subspaces are of the same dimension.

Suppose we have two sets of vectors G(α l ) and g(α l ) for α l ∈ I m and A = Span(G(α l ), α l ∈ I m ) and B = Span(g(α l ), α l ∈ I m ) two subspaces. The orthonormal basis {a i , i = 1, .., k} and {b i , i = 1, .., k} are defined by our incomplete Gram-Schmidt orthonormalization procedure of section 2.2 applied to {G(α l ), α l ∈ I m } and {g(α l ), α l ∈ I m }. We have therefore A = Span(a i , i = 1, .., k) and B = Span(b i , i = 1, .., k) but with k < m (the sampling size for α).

2.6.1. Principal angles as a measure of epistemic uncertainties. Suppose we have two solution procedures (u, j) and (ũ, j) leading after linearization to two sets of gradients:

(

G(α l ), g(α l )) = (∇ x j(α l , x), ∇ x j(α l , x)), (G(α l ), g(α l )) = (∇ α u(α l , x), ∇ α ũ(α l , x)
). The previous procedure for each set leads to two pairs of spaces A = s m (resp. A = S m ) and B = sm (resp. B = Sm ). Principal angles between the couples of subspaces (A, B) give quantitative indications of the impact of a change in the solution procedure on the global search spaces. In particular, putting extra effort in a modeling or solution procedure complexity would be useless if the angles remain small as the outcome, even if more expensive, will have no impact on a design using this simulation procedure. 2.7. Deviation in multi-point direction . An important situation where the calculation of the eigenvalues of M t M is not necessary is with multi-point analysis described in section 2.3 where the direction d has been introduced. Suppose one has two sets of gradients vectors G(α l ) and g(α l ) for α l ∈ I m , from the linearization of different state equations or numerical schemes leading to d and d using formula [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. We said that we do not consider what falls below the tolerance threshold T OL. Our incomplete Gram-Schmidt orthonormalization leads to subspaces A(G) or B(g) of dimensions k < m (again, the dimensions can be different). As {a i , i = 1, ..., k} and {b i , i = 1, ..., k} are orthonormal basis of A and B, projecting d on these spaces gives:

P r A (d) = i=1,...,k < d, a i > a i and P r B ( d) = i=1,...,k < d, b i > b i .
In a context of optimization, these are eventually the most important quantities and the impact of a change in the modeling should be measured through the deviation between d and d instead of through the angles between A and B. Therefore, extra modeling or numerical efforts has an impact if the angle ∠ P r A (d), P r B ( d) is large. This indicator is easy to estimate during optimization iterations [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. Even this is only an a posteriori indicator, it is interesting as in practice a given design procedure is often applied several times in a same context. The indicator is evaluated for the first design and provide a confidence level for possible use of low-order models in the further designs.

2.8. Link with classical UQ quantificators. Before illustrating these ingredients on our examples, let us summarize the philosophy of our approach and see where classical UQ quantificators, through intervals or moments of the outputs, appear.

In a simulation chain, the PDF of the outputs can be uniform in which case it can be described by an interval. It can be Gaussian and defined by its mean and variance. Or, it can be of another type needing a posteriori identification and involving more moments for its description.

Mean and variance have been defined for the aleatory variability due to the functioning parameter α as described in section 2.3 on the multi-point optimization algorithm. The multi-point direction d given by ( 5) permits their control. We get different values for the moments and directions for different models or solution procedures. And their variations provide a measure for the epistemic uncertainties.

The directional extreme scenarios are built using the PDF of the control parameters x in large dimension. The variability in x is a posteriori and aleatory. Indeed, once the design is made and an optimum found, one cannot be sure, for instance, that it can be manufactured exactly following its definition and we have an uncertainty here which appears a posteriori and which is not epistemic as it cannot be reduced improving our modeling or solution procedures. But, the PDF of the uncertainty is usually known, through manufacturing margin for instance. This permits to define directional margin intervals for each of the components of x as described in section 2.4. What makes things complex is that because the optimum depends on the modeling and solution procedures, these intervals will be different for different procedures. It is therefore impossible here to dissociate the aleatory and epistemic uncertainties.

Hence, the splitting of the control parameters into α and x is not motivated by their aleatory versus epistemic features. The splitting is motivated by the calculation complexity: α is of small dimension while x belongs to a large dimensional space. Therefore, whatever the PDF of α is, the moments of the state U and the functional j can be computed as a sampling of the parameter space for α (following its PDF) can be envisaged. Here we have illustrated the approach with uniform PDF but this is not a perquisite. On the other hand, it is not possible to proceed in the same way to account for the variability in x because of its dimension.

The paper proposes to go beyond simple mean/variance or intervals for the quantification of the outputs uncertainties. It analyzes the geometry of the global sensitivity spaces. This is because one cannot be confident, even with low aleatory and epistemic uncertainties, if the global sensitivity spaces are very different and, in particular, if the sensitivity spaces built around the extreme scenarios deviate (e.g. a situation of small variance with also small gradients but pointing in very different directions).

Finally, the different estimations proposed in the paper can be made in a fully non-intrusive way. We use in parallel the ingredients available in an existing monopoint optimization platform where the gradient of the functional is available through an adjoint formulation.

Uncertainty quantification for the Burgers equation

Let us apply the concepts presented above for quantification of some of the uncertainties for the solution of the following Burgers equation with a distributed control in its right-hand side:

∂ t u + ∂ y u 2 2 = x(y)u, (9) u
(α, t, -1) = α, u(α, t, 1) = -0.8, α ∈ I = [0.9, 1.1], u(α, t = 0, y) = u 0 (α, y) = - α + 0.8 2 y + α -0.8 2 ,
where the control x(y) = 0.3y. This simple model gathers two sources of uncertainties often present in practice. α control the inlet boundary condition and can be seen as the inflow Mach number in a flow simulation and this quantity might only be known through a probability function (PDF) or an interval (with uniform PDF). Uncertainties might also be on distributed data over the domain as in x(y) = 0.3y(1 + N (0, σ)) where a Gaussian PDF is assumed. Uncertainties over the definition of a shape studied in section 4 is a particular case of this latter. For σ = 0, the steady solution of ( 9) is piecewise parabolic and has a jump at y = s: u(y) = 0.15y 2 + α -0.15 for y < s, u(y) = 0.15y 2 -0.95 for y > s, [START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF] and the shock position is found by asking for the flux to have no jump:

(11) u - s = -u + s therefore s = - 1.1 -α 0.3 .
We solve equation ( 9) with a centred scheme for first and second order spatial derivatives after introduction of a numerical viscosity given by ν h = 0.5hu max where

u max = max(|u i-1 |, |u i |, |u i+1 |
) on a 100 points uniform mesh (i.e. h = 0.02). Time integration is based on a third order Runge-Kutta scheme with local times step given by h/u max . These are very basic choices and much sophisticated schemes could have been envisaged but this is not central to our discussion. Figure 4 shows the solution by this scheme versus the analytic expressions [START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF] for the two extreme values in α in its domain of variation I = [0.9, 1.1]. We called this the state equation solver and its solution the state variable U .

3.1. Global sensitivity space s m . To build the global sensitivity space s m we need u α (α l , x) for α l ∈ I m . But u has a shock. We need, therefore, to address the question of sensitivity evaluation in the presence of discontinuities [START_REF] Majda | The stability of multi-dimensional shock fronts[END_REF][START_REF] Godlewski | On the linearization of hyperbolic systems of conservation laws. Application to stability[END_REF][START_REF] Bardos | A Formalism for the Differentiation of Conservation Laws[END_REF][START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. Suppose the initial condition u(α, t = 0, y) = u 0 (y) has a discontinuity at s = 0 satisfying the entropy condition u -(α, 0, 0) > u + (α, 0, 0) (the velocity decreases after a shock). Then u(α, t, y) has a discontinuity at s(t) which depends on α and propagates at a velocity given by the Rankine-Hugoniot condition: ṡ = ū := (u + + u -)/2 where u ± denote its values before and after the shock.

Let us denote by H the Heaviside function and by δ the Dirac function, its derivative such that H α (y -s(t)) = -s α (t)δ(y -s(t)). u can be rewritten as u(α, y, t) = u -(α, y, t) + (u + (α, y, t) -u -(α, y, t))H(y -s(t)) which after differentiation in α gives: [START_REF] Shonkwiler | Poincare duality angles for Riemannian manifolds with boundary[END_REF] u α = u - α -s α (t)[u]δ(y -s(t)), where [u] = u + -u -is the jump of u across the shock and u - α the point-wise derivative of u -with respect to α. In practice we get u α from a linearized version of our simulation code for the solution of the Burgers equation ( 9) (see section 2). The question is then how to link this numerical derivative to the previous analysis ?

Differentiating [START_REF] Jordan | Essay on geometry in n dimensions[END_REF] in α we have:

(13) ∂ t u α + ∂ y (uu α ) = x(y)u α , u α (α, y, 0) = u 0 α (α, y) = - y 2 , u α (α, y = ±1, t) = 0.
From what said above one should not be allowed to do that as uu α in [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF] has no meaning at s(t) because it involves the product of a Dirac function with a discontinuous function. One can notice however for steady solutions (i.e. the time derivative vanishes) and outside shocks that the equation reduces to uu y = ux(y) which gives u y = x(y). After derivation in α we obtain u αy = 0. So u α must be constant outside shocks and as we have u 0 α = -y/2 + 1 from (9), we find u α (y) = 1 for -1 ≤ y < s and u α (y) = 0 for s < y ≤ 1 as shown in figure 4.

In shock regions, the classical solution is to say that ( 12) is valid at all points except at (t, s(t)), and that the Rankine-Hugoniot condition, differentiated, gives: [START_REF] Godlewski | On the linearization of hyperbolic systems of conservation laws. Application to stability[END_REF] ṡα (t) = ūα (s(t), t) + s α (t)∂ y ū(t, s(t)).

We see that the derivative of the Burgers equation has two unknowns: u α and s α .

The entropy condition insures uniqueness of u α given by ( 13) and a jump condition across the shock is not necessary because the characteristics left and right of the shock point left and right too. Once u α is computed, s α is given by the linear ordinary differential equation [START_REF] Godlewski | On the linearization of hyperbolic systems of conservation laws. Application to stability[END_REF]. But, again, in practice we never proceed in that way. In particular, s α is never evaluated and, as we said, u α is estimated, for instance, using finite differences or the complex variable method applied to the solution of the Burgers equation. One should therefore not have a complete confidence on u α in these regions. Again, one situation of interest is when the solution is steady in time and because s α = 1/2s from (11), we recover from ( 14) that in shock regions we have:

(15) ūα (s) = - 1 2s ∂ y ū.
Figure 4 shows the comparison between this expression and the outcome by finite differences for α = 0.9 and α = 1. One sees that the numerical gradient produces the right sign and also nearly the right amplitude for u α . Such sensitivity analysis can therefore be somewhat trusted in the construction of the global sensitivity space s m , at least for steady solutions. This space is built using the algorithm in section 2.1 and the incomplete Gram-Schmidt orthonormalization described in section 2.2. Figure 5 shows the outcome of this analysis for a m = 50 points sampling of I and for a 100 points uniform sampling of T OL ∈ [10 -6 , 1] in the incomplete orthonormalization.

For this 100 points mesh, the gradients u α are vectors of size 100 and one saw from our analysis and shown in figure 4 that these take for value 0 or 1, except in shock regions. The orthonormalization algorithm detects that and one sees that the 9) vs. the analytic expressions [START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF] for α = 0.9 and 1.1. Lower-right: u α for the two simulations versus u α (y) = 1 for -1 ≤ y < s and u α (y) = 0 for s < y ≤ 1. Lower-left: u α for α = 0.9 and α = 1 in shock regions versus expression [START_REF] Bardos | A Formalism for the Differentiation of Conservation Laws[END_REF] in dashed.

dimension of s m is two over a wide range of larger values of T OL. Increasing the accuracy (i.e. reducing T OL) in the orthonormalization permits to detect new free dimensions up to a maximum of 9. The sampling I 50 is therefore a safe choice. The new dimensions incrementally detected by the algorithm also characterize this multiscale feature of solution. One difficulty here is that these extra dimension can have unphysical origins and, in particular related to numerical noises for instance. We discussed this in sections 2.1 and 2.3.1 and also reported in [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. We will illustrate it in the next section for the definition of S m . The numerical noise might artificially increase the size of the global sensitivity space and should also be removed when the sensitivities are used for the definition of descent directions in a multi-point optimization context.

3.1.1. Multidimensional α. We worked with only one functioning parameter α. In the presence of more than one functioning parameters α i=1,... , a similar analysis can be carried out for each of them leading to spaces s i=1,... m and a ranking can be delivered on the variability each of these parameter can induce on the state through dim(s i=1,... m ). This information is interesting as the norms of the gradients alone would not necessary imply large variability. Actually, gradients can be large and still the global sensitivity spaces of low dimension. These constructions can be carried out efficiently if an adjoint formulation is available as we will see in the analysis of spaces S m in the next section. J(α, x(y), t, y)dydt,

where the control x(y) ∈ IR 100 is defined by its 100 values of the mesh points. One typical situation in inverse problems is with J = 1 2 (u(α, x(y), t) -u d (α, y)) 2 to measure the distance between the state variable and a desired state u d . We choose for u d the piecewise parabolic solution [START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF]. This is therefore an inverse problem and the existence of the solution is guaranteed with x d = 0.3y for all α.

Of course, the control can be recovered knowing u d in regions where the solution u has no shocks and where we can therefore write:

(17) u d x d (y) = ∂ t (u d ) + u d ∂ y (u d ) = u d ∂ y (u d ).
The second equality is because we are looking for a control independent of the time. So x d = ∂ y (u d ) which is known in regions where u d is smooth and we recover the target control which is independent of α. But, in the context of robust design, as described in section 2.3, this a priori information is usually unavailable.

In our model problem the size n of the control parameter x and the size N of the state variable u are the same (i.e. ∇ α u and ∇ x j are both vectors of size 100). Identification of space S m follows the same procedure than for s m in the previous section except that we need to introduce an adjoint formulation as the size of the control space is n while α was a scalar and could be accessed in realistic time with finite differences or the complex variable method. Figure 6 shows an example of the distribution of the adjoint based ∇ x j(α l , x) for α l ∈ I 10 (a 10 points uniform sampling of I) as described below.

3.2.1. The adjoint formulation. The derivative of j = (-1,1)×(0,T ) J with respect to x writes:

∇ x j = (-1,1)×(0,T ) (J x + J u u x )dtdy.
For our functional targeting u d , J x = 0. From the linearized Burgers equation we have (same linearization procedure than for u α in the previous section):

(18)

∂ t u x + ∂ y (uu x ) -x(y)u x = 0, u x (y, t = 0) = - α 0.6 , u x (y = ±1, t) = 0.
But, unlike with u α ∈ IR N , we need to avoid working with u x ∈ IR N ×N . One introduces therefore an adjoint state v solution of the backward equation:

(19) ∂ t v + u∂ y v -x(y)v = J u , v(y, t = T ) = v(y = ±1, t) = 0,
and write

(-1,1)×(0,T )

J u u x = (-1,1)×(0,T ) (∂ t v + u∂ y v -x(y)v)u x = - 1 -1 (u 0 ) x v(y, t = 0) = -v(y, t = 0),
where (u 0 ) x = δ(y) formally indicates the sensitivity with respect to the independent variable x of the initialization of the dependency chain (α, x) → u 0 → u → J. The adjoint state v has no shock because its time boundary condition is continuous and the characteristics integrated backward never cross the shock [START_REF] Giles | Analytic adjoint solutions for the quasi-one-dimensional euler equations[END_REF]. Lower:

d = ∇ x µ-< ∇ x µ, ∇ x σ * > ∇ x σ * +0.01∇ x σ * defined in section 2.3.
Having in hand the set of {∇ x j(α l , x), α l ∈ I 100 }, subspace S 100 can be built and its dimension analyzed. Figure 7 shows the history of the Gram-Schmidt orthonormalization for different values of T OL. One sees that the generator subspace in S 100 is at most of dimension four.

3.2.2.

Comparative studies between s m and S m . We notice that q = dim(S 100 ) < dim(s 50 ) = p. This characteristic is often encountered in optimization problems. It means that the functional is less sensitive to the variability of the functioning parameters than the state. It gives a justification for the use of reduced order models in the solution of the optimization problems. It tells, indeed, that it should be easier to build a low-order model for the functional than for the state. This indicator also permits to discriminate between several outputs of a system, giving a ranking in terms of their sensitivity to the variability of the functioning parameters of the system. Having said that, one notices however that we needed to use a 100 points sampling of the range of α for the analysis of S 100 . This was necessary to be on the safe side as the fourth dimension has been identified after more than 50 vectors considered. Hence, the 50 points sampling used in the s 50 analysis would have been insufficient. This is puzzling. Indeed, dim(s 50 ) > dim(S 100 ) but we see that we need a finer sampling to identify the free generator of the sensitivity space built with ∇ x j(α l , x) than the space built with ∇ α u(α l , x). This is an indication of the care one needs to take when building a low-order model to be used in a minimization problem for j.

In each case, we have verified that enriching more the sampling in α does not produces any new independent direction with respect to s 50 and S 100 respectively. But, it would be interesting to know if there are any a priori upper bounds for the size of the sampling in α. Here, the control space for x is of dimension 100 (x(y i ) on y i=1,...,100 the 100 points spatial mesh). This example shows a situation where the necessary sampling is the one given by the worst-case theoretical bound given in [START_REF] Li | Robust Airfoil Optimization to Achieve Consistent Drag Reduction Over a Mach Range[END_REF][START_REF] Gallard | Robust parametric design by multipoint optimization[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF]. This bound is related to the fact that one needs n+1 independent vectors to generate a vector space of dimension n. Fortunately, this situation is quite rare and we show in section 4 an example where dim(S m ) << n. However, this example suggests that when defining an optimization problem and if several formulations of the minimization problem are available, one should advantage the functional (even through a reformulation of the problem) for which the global sensitivity space S m is less sensitive to the size of the sampling. Figure 8 shows the same analysis than in picture 7 but for a different functional where J only involves the local derivative in space at one of the boundaries: J = (∂ y u| y=-1 ) 2 . This is similar to the situation we will have with the aerodynamic coefficients in example 4. So, not only the worst-case upper bound can be avoided (which is fortunate), but it should be. Indeed, here the size of the two gradients ∇ x j(α l , x) and ∇ α u(α l , x) are the same and equal to the size of the mesh because the control x is distributed. But, in a three dimensional simulation, for instance, the upper bound would suggest again a size of the sampling proportional to the size of the spatial mesh which is simply out of question because of the related calculation complexity. Shape optimization problems, such as in section 4, are particular situations where the control is only distributed on a portion of the domain: a surface plunged in it. Still the upper bound would be too constraining.

Let us make a final remark on the comparative costs of the two constructions s 50 and S 100 . The former requires 50 solutions of our Burgers solver and finite differences to build gradients ∇ α u(α l , x). The latter 100 solutions of our Burgers solver and its adjoint. Without the adjoint, the construction of S 100 would have requested 10 4 solutions of the Burgers equation which is unreasonable even for this 1D case.

We could have avoided all this discussion working with a 100 points sampling and s 100 . But, we think it is important to bring into light such differences. It is clear that we will not be able to use exhaustive sampling of the functioning parameters ranges in real life applications (like in section 4). It is therefore important to be aware of the risks related to a partial sampling. This weakness is present in any sampling based statistical approach and the mentioned worst-case upper bound is useless when the size of the control space is large.

Worst-case sensitivity spaces s ±

m and S ± m . We showed how to introduce dimensional geometric characterizations of the global sensitivity spaces obtained using the gradients of the state ∇ α u(α l , x) or a functional ∇ x j(α l , x) for α l ∈ I m . Now, let us introduce the Value at Risk (VaR) concept to define bounds for the dimensions of our global sensitivity spaces in a context of directional worst-case analysis as described in section 2.4.

Consider our Burgers equation with control where we assume the control uncertain and described through a same Gaussian probability density function for all the components of x ∈ IR n where n = 100 for our 100 points mesh: x(y i ) = 0.3y i + N (0, σ) with σ = 0.1.

We use this information together with the direction d given by ( 5) to define extreme values for the local variability of the functional along d as shown in figure 3 and described in section 2.4. Global sensitivity spaces s ± 50 and S ± 100 are then at points x ± using ∇ α u(α l , x ± ) and ∇ x j(α l , x ± ) for α l respectively in I 50 and I 100 . Figure 9 shows the dimensions p ± and q ± at different T OL of the free subspaces generating s ± 50 and S ± 100 . In this case, p + and q + appear often larger but this is not necessary the case in general.

3.4. Variability on subspaces s m and S m due to a change in the modeling or solution procedure. We saw how to build subspaces s m and S m for a given sampling I m of the range of a functioning parameter and also how to define directional extreme scenarios (DES) leading to spaces s ± m and S ± m . And, we saw how to quantify the variability of these global search spaces through the dimension of the free generator subspaces.

But these depend on several parameters such as the numerical methods used in the calculation of the state variables or for the definition of the gradients. Indeed, one simple question, for instance, can be the impact on the global sensitivity spaces of . Dimensions p ± (upper) and q ± (lower) at different T OL of spaces s ± 50 and S ± 100 , respectively generated by ∇ α u(α l , x ± ) and ∇ x j(α l , x ± ) for α l respectively in I 50 and I 100 . the accuracy of the gradients ∇ α u(α l , x) and ∇ x j(α l , x). Also, there is often large discussions on the approach which should be adopted for their definition (often discussed under the items linearize and discretize or discretize and linearize).

Our aim is to give a quantitative way to see how much deviation a change in the modeling complexity or numerical scheme for larger accuracy, and therefore generally computational complexity, brings into these global sensitivity spaces. This measure of the deviation comes in addition to possible variations in the dimension of the respective sensitivity spaces. Principal angles between subspaces are an interesting instrument to quantify such deviations.

Following the procedure presented in section 2.6, figure 10 shows the angles between global search spaces s m and sm with first and second order accurate approximations of gradients ∇ α u(α l , x) ∈ IR 100 for α l ∈ I 50 . After orthonormalization with T OL = 10 -4 , the two subspaces free generators appeared being of the same dimension 7. This is not necessarily the case as we will see for the sensitivity spaces S 100 and S100 . We see that the two subspaces have two common axis but are sometimes nearly orthogonal. This is a measure of the error one commits in the global sensitivity space definition if one uses lower accuracy in the gradients. This also means that these directions are very sensitive to small perturbations of the functioning parameters.

Figure 11 shows the angles between global search spaces S 100 and S100 built with the gradients ∇ x j(α l , x) respectively with a discrete adjoint approach and with first order accurate gradients by finite differences. As discussed in section 3.2.2, this latter is very costly even for this model problem as it requires here 10 4 solutions of the direct problem (x is distributed over the 100 mesh points and we have 100 sampling points in α) while the adjoint solutions cost about 100 solutions (for 100 points sampling in α) of the direct problem and this cost is independent of the number of the control points (here mesh points). After orthonormalization with T OL = 10 -4 , the dimension of the S 100 is 4 but S100 appears being of dimension 18. This difference is related to the regularity of the two gradients as shown in the figure for α = 0.9. The fluctuations in the gradient by finite differences can be reduced with a second order finite difference approximation or even removed with the complex variable method, but our aim is to show a common situation where they exist. This is an example of principal angles calculation between subspaces of different dimension. Again, this permits to measure the impact of a loss in accuracy in the solution procedure. This is also an illustration of the mentioned numerical artifact justifying the use of incomplete Gram-Schmidt procedures [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. One also sees that despite the global sensitivity space S 100 is of smaller dimension than s 50 , Upper: ∇ x j(α l , x) for α l = 0.9 respectively with a discrete adjoint approach (continuous) and with first order accurate gradients by finite differences. Lower: principal angles between global search spaces S 100 and S100 built with these gradients. its sensitivity to a loss of accuracy is larger as there is, for instance, no common axis between S 100 and S100 . This also suggests that a robust optimization problem with an admissible solution in the search space defined by the adjoint will not probably be reached by a minimization method using a descent direction in S100 .

Uncertainty quantification for a flow solution over an aircraft

After the model problem with the Burgers equation, we would like to apply the different geometric uncertainty quantification ingredients to a realistic 3D simulation for a full aircraft in transonic cruise condition. Several sources of variability exist for these problems due to a change in the weight of the aircraft during the flight or due to variability in the flight conditions. One example of the latter is when the aircraft cruises against transverse winds which are very common. Usually aircraft are designed for a range of angle of incidence. These designs are usually realized with the sideslip angle set to zero. 4.1. Using ingredients of a shape optimization platform. We work in the framework of an existing shape optimization platform which has been previously presented in [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF]. We use, in particular, three of its simulation codes for the shape parameterization, for the fluid dynamic calculations and for the shape sensitivity analysis of aerodynamic coefficients. This is a very standard and generic situation and these particular choices are not central to our discussion.

Let us briefly recall our direct dependency chain linking independent variables (α, x) to dependent variables (Q, U ) describing geometrical entities and state variables and to the cost function j and constraints C i , i = 1, ..., n c .

(21) (α, x) → (α, Q(x)) → U (α, Q(x)) → (j, C i )(α, x, Q(x), U (α, Q(x))).
x denotes a CAD-free parameterization [START_REF] Mohammadi | Shape Optimization in Fluid Mechanics[END_REF] which does not require a priori local regularity assumptions on the shape as it is implicitly the case in CAD-based shape definitions. More precisely, x represents shape deformations along the normal to the triangular faces of the surface mesh as shown in figure 13. For the problem discussed here this search space has a dimension n ∼ 5000. This parameterization receives different denominations and belongs to the same class as node-based or free-form shape definitions. In all these approaches the regularity of the deformation needs to be controlled [START_REF] Firl | Regularization of shape optimization problems using fe-based parametrization[END_REF][START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. α is the sideslip angle inducing fully 3D effects on the flow around the plane making the consideration of a full aircraft necessary for the analysis. However, because the airplane geometry is symmetric spanwise, it is not necessary to consider a symmetric range for the transverse wind. We consider the sideslip angle α in I = [0, 10 o ]. In our multi-point analysis, we consider a uniform sampling I m of [0, 10 o ] with m = 30 points. The other flow conditions are given by a Mach number of 0.8 and zero inflow incidence. Together with the sideslip angle these parameters fully describe a 3D inviscid flow around the aircraft.

Q(x) denotes the auxiliary unstructured mesh related geometrical quantities. U (α, Q(x), ) denotes flow variables depending also on the extra parameter α not part of those involved in the definition of the shape. More precisely, U (α, Q(x)) = (ρ, ρ u, ρE) t represents the conservative flow variables solution of the Euler equations in conservation form. Our flow solver is based on a finite volume Galerkin method on unstructured tetrahedral meshes. The details of the implementation are available in [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF][START_REF] Mohammadi | Fluid dynamics computation with NSC2KE[END_REF]. Of course, other choices are possible for the flow solver and the literature on numerical methods for compressible flows is huge. Let us briefly describe the different ingredients in our solver. The Roe [START_REF] Roe | Approximate Riemann solvers, parameters vectors and difference schemes[END_REF] flux is used for the approximation of the advection operator together with MUSCL reconstruction with Van Albada limiters [START_REF] Van Albada | Flux vector splitting and Runge Kutta methods for the Euler equations[END_REF] in presence of shocks. We target steady solutions and use time marching with local time steps to reach these. The local time steps are proportional to h/( u + c) with h the local mesh size and c the local speed of sound. The time integration procedure is explicit and is based on a low-storage Runge-Kutta scheme (RK4). To illustrate the level of convergence to the steady state we use here a typical convergence history is shown in figure 12. As we are interested by aerodynamic coefficients, the sufficient level of convergence retained for the flow is when these coefficients are considered to be converged. This is an important point as one needs to particularly monitor the calculation complexity for our multi-point shape optimization algorithm to be efficient, in particular if one intends to quantify the uncertainties in reasonable time. We will see that limiting the number of forward iterations is also useful for efficient adjoint calculation. We consider a classical aerodynamic problem where two main quantities of interest are the drag C d and lift C l coefficients:

(22) C d (α, x) = 1 2ρ ∞ (α) u ∞ (α) 2 shape(x) p(α, Q(x))(u ∞ (α).n(Q(x))dγ,
where superscript ∞ indicates inflow conditions. The lift coefficient is evaluated with formula [START_REF] Hascoet | Tapenade 2.1 user's guide[END_REF] where u ∞ (α) is replaced by u ⊥ ∞ (α) in the boundary integral. Aircraft performance analysis concerns its payload and range. These are directly linked to the aerodynamic coefficients of the aircraft called the lift (conditioning the payload) and drag (conditioning the fuel consumption) coefficients. The lift coefficient often appears through an inequality C l -C target l ≥ 0 or equality constraint C 1 = (C target l -C l ) 2 with C target a target performance. Let us consider this second situation.

Structural efficiency and necessity of useful free volume also implies the consideration of geometric criteria such as a constraint on the volume V of the aircraft or its by-section definition. The volume of an object Ω (here the aircraft) is expressed through the boundary integral formula: V = Ω 1 = Ω 1 3 ∇.( X) = ∂Ω X. n, where X = (x 1 , x 2 , x 3 ) t is the local coordinate over the shape.

We would like to use the ingredients of the paper and the global sensitivity space S 30 for the analysis of this aerodynamic performance problem. Two approaches are possible. Either introduce a functional j involving the drag coefficient and the constraint on the lift coefficient accounted for by penalty: j = C d + aC 1 , a ∈ IR + and define the sensitivity space S 30 for j (as in section 3.2). Or define two global sensitivity spaces S C d 30 and S C l 30 for each of the aerodynamic coefficients. An important remark is that there is no need to account for geometric constraints in the analysis as those are independent of α.

The second approach requires the solution of two adjoint problems to get ∇ x C d and ∇ x C l but it gives valuable information on each of the sensitivity spaces S C d 30 and S C l 30 and their relative position. A lower dimension for one would mean that the corresponding coefficient is less sensitive to α. In the same way larger or smaller principal angles would tell us how antagonist the behaviors of the two coefficients are.

In all case, because these are vector spaces and Therefore, the first approach gives upper uncertainty quantification bounds and will be considered in the sequel.

∇ x j = ∇ x C d +2a(C target l -C l )∇ x C l is linear combination of vectors in S

4.2.

Global sensitivity space S 30 . To proceed with our global sensitivity space S j 30 (the j upper-script will be omitted in the sequel) analysis we need to provide ∇ x j(α l , x) for α l ∈ I 30 as in our model problem with the Burgers equation. Functional sensitivities with respect to the shape ∇ x j are computed by automatic differentiation in reverse mode [START_REF] Griewank | Computational derivatives[END_REF] using tapenade [START_REF] Hascoet | Tapenade 2.1 user's guide[END_REF]. As we said, our direct Euler solver uses time marching to steady solutions. An interesting optimization of the reverse mode comes from the fact that, our situations of interest being stationary in time, there is no need in storing the forward states for backward integration [START_REF] Christianson | Reverse accumulation and implicit functions[END_REF][START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. All intermediate states can be replaced by the final converged one. Still, it is good to keep the number of forward time iterations low as discussed earlier as this gives also an indication of the number of backward iterations in the adjoint code.

Figure 13 shows the shape and a view of the triangular surface mesh for flow calculations. This mesh is also considered as the parameter x. The figure also shows four snapshots of ∇ x j(x, α l ) for four values of α l = 2, 5, 7 and 10 o . These are some of the gradients used to define the multi-point direction d using [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF].

Figure 14 shows the impact of the sorting (from 1 to m and from m to 1) of the vectors ∇ x j(α l , x) during the incomplete Gram-Schmidt orthonormalization on the dimension of the subspace generator in S 30 for different values of the threshold T OL. The figure also shows the histories of all orthonormalizations for the different thresholds. As we would like to avoid the directions related to possible numerical artifacts, the threshold retained for the analysis is T OL = 10 -4 for which we have dim(S 30 ) = 10. These dimensions have been found rapidly during orthonormalization. We expect therefore to be in a situation discussed in section 3.2.2 making the size of the sampling sufficient to identify the global sensitivity space at this threshold.

CAD-free parameterization (x) &

CFD surface mesh (on which Q(x) is defined) Grad x j(x, alpha l ) for 4 alpha l CAD-based aircraft shape Figure 13. Aircraft shape and a view of the triangular surface mesh for flow calculation which also serves for the definition of the CAD-free parameterization. ∇ x j(α l , x) for four values of α l = 2, 5, 7 and 10 o .

4.3. Worst-case sensitivity spaces S ± 30 . We assume a priori local Gaussian Value at Risk information for the shape of the aircraft x = x 0 + N (0, σ(x 0 )) with σ(x 0 ) variable in space over the shape as shown in figure 15. We recall that this is in a space of dimension 5000 (along local unit normal vectors to the surface mesh). We assume that this distribution is symmetric spanwise and that larger uncertainties are along the wings and increasing spanwise. This variability can be, for instance, related to in-flight icing or elastic deformation of the shape during flight. Hence, at given probability, this VaR models the maximum local shape variation normal 5) quantify the directional sensitivity of S 30 at given T OL through q ± = dim(S ± 30 ) and the principal angles between S ± 30 .

to the shape with respect to a reference configuration x 0 . We follow the approach presented in sections 2.4 and 3.3 to see the bounds for the dimension of spaces S ± 30 defined using ∇ x j(α l , x ± ) for α l ∈ I 30 where x ± is defined from (8) along the direction d given by d = ∇ x µ-< ∇ x µ, ∇ x σ * > ∇ x σ * + η∇ x σ * . This is shown in figure 15. Because we want the shape to remain symmetric spanwise, this constraint has been included in the definition of d.

One can notice that unlike in figure 9 where q + = dim(S + 30 ) appeared often larger that q -= dim(S - 30 ) this is not the case here. Once the spaces S ± 30 identified, principal angles between them give a second set of quantitative information on the impact of the geometric variability on the global sensitivity space S 30 . Here, one sees that the sensitivity spaces have different dimensions and that they share only one common direction. The other directions deviate with angles up to 22 degrees. This analysis shows the importance of including possible variability of a shape (e.g. due to its elastic features or in-flight icing) during its design. It highlights the necessity of proceeding with multi-disciplinary optimization in such situations [START_REF] Martins | Multidisciplinary design optimization: A Survey of architectures[END_REF]. These probabilistic extreme values for shape deformations permit worst-case analysis for the sensitivity of a design to possible perturbations.

4.4.

Variability on S 30 due to mesh refinement. One interest of the ingredients of the paper is to permit to analyze the impact of different space or time discretizations on the solution of the state equations. Indeed, the evolution of the dimension of S 30 and the principal angles will tell us how useful the introduction of a more sophisticated and costly numerical scheme or discretization is. Indeed, if one sees no change in these parameters, this would mean that it is useless to spend more effort in that direction because the outcome will have no impact on the global sensitivity space where the design takes place.

To illustrate this discussion, we analyze the impact on the global sensitivity space S 30 of two levels of discretization. The meshes have respectively 200.000 and 400.000 elements. These are quite coarse meshes and we are probably far from mesh independence for the flow solution. Still, this analysis permits to quantify how much a given mesh refinement is useful for our aerodynamic performance analysis problem. This permits to define a new mesh independence concept which is actually the truly important one: the impact of mesh refinement on the global search space where the design takes place. Figure 16 shows the multi-point direction d built for our two meshes using ∇ x j(α l , x) for α l ∈ I 30 evaluated on each mesh. One sees that d is much richer for the fine mesh and features small scale details which are not present in the direction built for the coarse mesh. This explains why very large principal angles are detected as the corresponding directions do not exist in the global sensitivity space for the coarse mesh.

Concluding remarks

Geometric characteristics of global sensitivity spaces have been used for uncertainty quantification addressing both issues of aleatory and epistemic uncertainties. The sensitivity spaces have been built for the state variables and also when specific quantities of interest are targeted. The impact of the variabilities on the data, the model or the solution procedures have been measured through the variation of the dimension of the global sensitivity spaces built for each situation. Beyond this dimensional analysis, the relative positions of the different spaces have been analyzed Angles at given TOL Angles (Deg) q = dim ( S m (fine mesh) )

Small scale differences for d between the coarse and fine meshes.

Figure 16. Impact of mesh refinement on S 30 . Upper-left: d on the fine mesh features small scales. At the threshold chosen the dimension of the global sensitivity space is larger and large principal angles are detected.

and possible deviations between sensitivity spaces quantified using the concept of principal angles between subspaces. To provide a viable solution in term of computational complexity, a sampling of the range of variation of the design parameters has been avoided when working on large dimensional control space. This has been made possible by the introduction of quantiles, such as the Value at Risk, and use of directional information from a moment-based search direction. Together these ingredients provide a framework for directional uncertainty quantification of extreme (worst-case) scenarios. This approach appear being suitable for the evaluation of the pertinence of a given modeling, or an increase in the modeling or discretization complexity, in the context of simulation under uncertainty. It is also suitable because it is non-intrusive and uses existing ingredients in mono-point optimization platforms in a fully parallel manner.

Figure 1 .

 1 Figure 1. Domain of variation for 1 -ζ 2 + ηζ when ζ ∈ [-1, 1] and η ∈ [0, 1].

Figure 4 .

 4 Figure 4. Upper: u(α, t → ∞, y) solution of (9) vs. the analytic expressions[START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF] for α = 0.9 and 1.1. Lower-right: u α for the two simulations versus u α (y) = 1 for -1 ≤ y < s and u α (y) = 0 for s < y ≤ 1. Lower-left: u α for α = 0.9 and α = 1 in shock regions versus expression[START_REF] Bardos | A Formalism for the Differentiation of Conservation Laws[END_REF] in dashed.

Figure 5 .

 5 Figure 5. Upper: dim(s 50 ) during incomplete orthonormalization for different T OL. Lower: dim(s 50 ) vs. T OL.

3. 2 .

 2 Global sensitivity space S m . Now suppose the quantity of interest j is of the form:

Figure 6 .

 6 Figure 6. Upper: ∇ x j(α l , x) with adjoint formulation for α l ∈ I 10 . Lower:d = ∇ x µ-< ∇ x µ, ∇ x σ * > ∇ x σ * +0.01∇ x σ * defined in section 2.3.

Figure 7 .

 7 Figure 7. Cost function (16) with J = 1 2 (u(α, x(y), t) -u d (α, y)) 2 . Upper: dim(S 100 ) during incomplete orthonormalization for different T OL. Lower: dim(S 100 ) vs. T OL.

Figure 8 .

 8 Figure 8. Cost function[START_REF] Griewank | Computational derivatives[END_REF] with J = (∂ y u| y=-1 ) 2 . Upper: dim(S 100 ) during incomplete orthonormalization for different T OL. Lower: dim(S 100 ) vs. T OL. The global sensitivity subspace is less sensitive to the size of the sampling than in figure7.

Figure 9

 9 Figure 9. Dimensions p ± (upper) and q ± (lower) at different T OL of spaces s ± 50 and S ± 100 , respectively generated by ∇ α u(α l , x ± ) and ∇ x j(α l , x ± ) for α l respectively in I 50 and I 100 .

Figure 10 .

 10 Figure 10. Principal angles between global search spaces s 50 and s50 built respectively with first and second order accurate approximations of ∇ α u(α l , x) for α l ∈ I 50 .

Figure 11 .

 11 Figure 11. Upper: ∇ x j(α l , x) for α l = 0.9 respectively with a discrete adjoint approach (continuous) and with first order accurate gradients by finite differences. Lower: principal angles between global search spaces S 100 and S100 built with these gradients.

Figure 12 .

 12 Figure 12. Convergence histories during time marching to steady state. Norm of the vector of conservation variables and C d and C l aerodynamic coefficients.

Figure 14 .

 14 Figure 14. Upper: impact on dim(S 30 ) of the order in which the vectors ∇ x j(α l , x) have been considered in the incomplete Gram-Schmidt orthonormalization vs. T OL. Lower: histories of dim(S 30 ) during orthonormalizations for the different threshold.

Figure 15 .

 15 Figure[START_REF] Bardos | A Formalism for the Differentiation of Conservation Laws[END_REF]. VaR(x) and direction d given by (5) quantify the directional sensitivity of S 30 at given T OL through q ± = dim(S ± 30 ) and the principal angles between S ± 30 .

This is interesting in our multi-point optimization context where we can exhibit orthonormal basis of the global search space for the multi-point optimization problem using Gram-Schmidt orthonormalization.

Acknowledgements: Adjoint codes for both the model problem with the Burgers equation and with the Euler equations have been obtained using Tapenade AD tool developed at INRIA-Sophia Antipolis by L. Hascoet and his team Tropics.


Franz Chouly
email: franz.chouly@univ-fcomte.fr

Patrick Hild
email: patrick.hild@math.univ-toulouse.fr

Yves Renard
email: yves.renard@insa-lyon.frurl

A NITSCHE FINITE ELEMENT METHOD FOR DYNAMIC CONTACT

Keywords: unilateral contact, elastodynamics, finite elements, Nitsche's method, timemarching schemes, stability

We present a new approximation of elastodynamic frictionless contact problems based both on the finite element method and on an adaptation of Nitsche's method which was initially designed for Dirichlet's condition [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF]. The corresponding space semi-discretized weak form reads:

               Find a displacement u h : [0, T ] → V h such that for t ∈ [0, T ] : ρü h (t), v h + A Θγ h (u h (t), v h ) + Γ C 1 γ h [P γ h (u h (t))] + P Θγ h (v h ) dΓ = L(t)(v h ), ∀ v h ∈ V h , u h (0, •) = u h 0 , uh (0, •) = uh 0 . (1)
In the above formulation V h is a finite element space built from standard Lagrange finite elements, piecewise linear or quadratic,

A Θγ h (u h , v h ) := Ω σ(u h ) : ε(v h ) dΩ - Γ C Θγ h σ n (u h )σ n (v h ) dΓ and P Θγ h (v h ) := v h n -Θγ h σ n (v h ).
The linear form L(•) stands for prescribed body and boundary forces. The domain of the elastic body is denoted by Ω and the contact boundary by Γ C ; T is the final time of simulation ; ρ is the density of the elastic material ; the notation v h n stands for the normal component on Γ C of v h , and σ n (v h ) is the normal stress on Γ C ; u h 0 (resp. uh 0 ) is an approximation of the initial displacement u 0 (resp. the initial velocity u0 ). The notation [•] + stands for the positive part of a scalar quantity, and •, • stands for the L 2 (Ω) inner product. The parameter γ h is a positive piecewise constant function on the contact interface Γ C (γ h (x) = γ 0 h T (x), where γ 0 is a positive constant and h T (x) is the size of the mesh element T ). Note the additional numerical parameter Θ which can be freely chosen in R. As in Nitsche's method for (static) unilateral contact, values of interest for Θ are -1, 0, 1 [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF].

A main interesting characteristic is that this approximation produces well-posed space semi-discretizations, for any value of γ 0 > 0, contrary to standard finite element discretizations. We study associated energy conservation properties and manage to prove that :

d dt E h Θ (t) = (Θ -1) Γ C 1 γ h [P γ h (u h (t))] + uh n (t) dΓ where E h Θ (t) := E h (t)-ΘR h (t)
is an augmented energy associated to the semi-discrete solution u h (t). The mechanical energy is given by

E h (t) := 1 2 ρ uh (t) 2 0,Ω + 1 2 a(u h (t), u h (t)) and R h (t) := 1 2 γ h 1 2 σ n (u h (t)) 2 0,Γ C -γ h -1 2 [P γ h (u h (t))] + 2 0
,Γ C is an extra term, which represents, roughly speaking, the non-fulfillment of the contact conditions at the semidiscrete level. Note in particular that the symmetric variant (Θ = 1) conserves the discrete augmented energy E h Θ .

Various time-discretizations for (1) are then considered, with the families of θ-schemes and Newmark schemes. We also introduce a new hybrid scheme which is second-order accurate and numerically stable without any restriction on the time-step in the case Θ = 1. This new scheme is inspired from [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], and introduces much less dissipation than unconditionally stable variants of θ-schemes and Newmark schemes. We study theoretically the well-posedness of each discrete scheme as well as its energy conservation properties. We finally achieve the corresponding numerical experiments.


C Casenave
email: casenave@laas.fr

E Montseny
email: emontseny@laas.fr

Dissipative state formulations and numerical simulation of a porous medium for boundary absorbing control of aeroacoustic waves

Keywords: porous medium, aeroacoustic waves, boundary control, diffusive representation, state representation, non rational operators, energy functional, energy balance

The problem under consideration relates to a model of porous wall devoted to aircraft motors noise reduction. For such a medium, the parameters of the propagation equation depend on the frequency, so the corresponding time-model involves non rational convolution operators. Consequently, the impedance of the wall is a non rational function of the frequency. On the basis of complex analysis and causality properties, we introduce infinite dimensional state formulations of these operators. The coupling of the so-obtained model with a standard aeroacoustic one then leads to a time-local system whose analysis is simplified thanks to the existence of an energy functional in the sense of which the global dissipativity is insured. Some numerical results are given to illustrate the theoretical results.

INTRODUCTION

Aircraft motors noise reduction is currently an important challenge for aerospace industry. In the case of hot zones such as exhaust nozzles, a specific porous wall was proposed in [START_REF] Gasser | Etude des propriétés acoustiques et mécaniques d'un matériau métallique poreux à base de sphères creuses de nickel[END_REF] for absorption of a wide part of the energy of incident acoustic waves. The frequency model of such a material which describes the wave propagation inside the porous medium, is given on (ω, z, x) ∈ R×]0, 1[×]0, X[ by: e(x) iω ρ eff (iω) û + ∂ z P = 0 e(x) iω χ eff (iω) P + ∂ z û = 0,

where u and P designate the Fourier transforms of the velocity and the pressure in the porous medium, e(x) denotes the thickness of the porous wall and the parameters ρ eff (iω) and χ eff (iω) are respectively the so-called effective density of Pride et al. [START_REF] Pride | Drag forces of porous-medium acoustics[END_REF]) and the effective compressibility of Lafarge [START_REF] Lafarge | Propagation du son dans les matériaux poreux à structure rigide saturés par un fluide viscothermique[END_REF]). These parameters are expressed [START_REF] Gasser | Etude des propriétés acoustiques et mécaniques d'un matériau métallique poreux à base de sphères creuses de nickel[END_REF]):

⎧ ⎨ ⎩ ρ eff (iω) = ρ (1 + a (1+b iω) 1 2 iω ) χ eff (iω) = χ (1 -c iω iω+a (1+b iω) 1 2 ), (2)
with:

ρ = ρ 0 α ∞ , χ = 1 P 0 , a = 8μ ρ 0 Λ 2 , a = 8μ ρ 0 Λ 2 , b = 1 2a , b = 1 2a , 0 < c = γ -1 γ < 1,
where the physical parameters ρ 0 , P 0 , μ, γ, α ∞ , Λ, Λ are respectively the density and pressure at rest, the dynamic viscosity, the specific heat ratio, the tortuosity, the high frequency characteristic length of the viscous incompressible problem and the high frequency characteristic length of the thermal problem.

The model ( 1) is completed by the following limit conditions: u |z=1 = 0 (wave reflection at z = 1) P |z=0 = w.

(3) When coupling this system with the fluid medium, two connection conditions at the interface Γ are necessary:

P fluid |Γ = P |z=0 and u fluid |Γ • n = φ u |z=0 ,
where n is the outgoing unit normal on Γ and φ > 0 is the porosity coefficient of the material [START_REF] Gasser | Etude des propriétés acoustiques et mécaniques d'un matériau métallique poreux à base de sphères creuses de nickel[END_REF]). Then, w and y := u |z=0 (4) can be interpreted respectively as the input and the output of (1), the correspondence w → y defining the impedance operator of the porous wall.

Parameters ρ eff and χ eff depend on the frequency ω, so the time model obtained by inverse Fourier transform of (1), is written:

e(x) (∂ t • ρ eff (∂ t )) u + ∂ z P = 0 e(x) (∂ t • χ eff (∂ t )) P + ∂ z u = 0, (5)
where ρ eff (∂ t ) and χ eff (∂ t ) are the (causal) convolution operators associated to the symbols ρ eff and χ eff respectively. Model ( 5) is not time-local because of the presence of convolution operators ρ eff (∂ t ) and χ eff (∂ t ).

We propose in section 2 a new formulation of such operators by means of a suitable diffusive input-output dif-ferential system. This formulation allows us to obtain, in section 3, a time-local augmented input-output model equivalent to (1,3,4). We show that in the sense of an explicitly known functional, this model is consistent from the energy dissipation point of view. Then in section 4, the dissipativity of the whole problem is easily deduced from a suitable global energy functional. In section 5, we compute the analytic expression of the porous wall impedance from which we build a reduced diffusive state formulation of the associated operator. Finally, some numerical results and simulations are presented in section 6 .

DIFFUSIVE REALIZATION OF CAUSAL CONVOLUTION OPERATORS

In this section, we present a particular case of a methodology introduced and developed in [START_REF] Montseny | Représentation Diffusive[END_REF] in a general framework.

We consider a causal operator defined, on any continuous function w : R + → R, by

w → t 0 h(t -s) w(s) ds. ( 6
)
We denote H the Laplace transform of h and H(∂ t ) the convolution operator defined by (6).

Let w t (s) = 1 ]-∞,t] (s) w(s) and w t (s) = w t (ts). From causality of H(∂ t ), we deduce: H(∂ t )(ww t ) (t) = 0 for all t; then, we have for any continuous function w:

(H(∂ t )w)(t)= L -1 (H Lw) (t)= L -1 H Lw t (
H(p) Ψ w (t, p) dp. ( 10
)
We denote Ω the holomorphic domain of H. Let γ a closed 1 simple arc in C -; we denote Ω + γ the exterior domain defined by γ, and Ω - γ the complementary of Ω + γ . By use of standard techniques (Cauchy theorem, Jordan lemma [START_REF] Lavrentiev | Méthodes de la théorie des fonctions d'une variable complexe[END_REF])), it can be shown:

Lemma 2. For γ ⊂ Ω such that H is holomorphic in Ω + γ , if H(p) → 0 when p → ∞ in Ω + γ , then: (H(∂ t )w) (t) = 1 2iπ γ H(p) Ψ w (t, p) dp, ( 11
)
where γ is any closed simple arc in Ω + γ such that γ ⊂ Ω - γ .

We now suppose that γ, γ are defined by functions of W 1,∞ loc (R; C), also denoted γ, γ. Under hypothesis of lemma 2, we have [START_REF] Montseny | Représentation Diffusive[END_REF]):

1 Possibly at infinity Theorem 3. If the possible singularities of H on γ are simple poles or branching points such that |H • γ| is locally integrable, then:

1. with ν = γ 2iπ H • γ and ψ(t, .) = Ψ w (t, .) • γ: (H(∂ t )w) (t) = R ν(ξ) ψ(t, ξ) dξ; (12) 2. with γn → γ in W 1,∞ loc and ν = γ 2iπ lim H • γn in the sense of measures: (H(∂ t )w) (t) = ν, ψ(t, .) , ( 13
)
where ψ(t, ξ) is solution of the evolution problem on (t, ξ) ∈ R * + ×R:

∂ t ψ(t, ξ) = γ(ξ) ψ(t, ξ) + w(t), ψ(0, ξ) = 0. ( 14
)
Remark 4. In the sequel, we will indifferently denote ν, ψ or ν ψ dξ the duality scalar product between a continuous function ψ and a measure ν.

Finally, in the particular case γ(ξ) = -|ξ|, we deduce from symmetry of the problem that there exists a measure μ such that

+∞ -∞ ν ψ dξ = +∞ 0 μ ψ dξ.
The state equation ( 14) is infinitedimensional. To get numerical approximations, we consider a discretization (ξ l ) l=1:L of the variable ξ and approximations μ L of the γ-symbol μ using atomic measures and of the form:

μ L = L l=1 μ l L δ ξ l .
Let us denote M L the space of atomic measures on the mesh {ξ l } l=1:L . If ∪ L M L is dense in the space of measures, we have:

μ L , ψ -→ L→+∞ μ, ψ ∀ψ.
So we get the finitedimensional approximate state formulation of H(∂ t ) :

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ t ψ(t, ξ l ) = γ(ξ l ) ψ(t, ξ l ) + u(t), l = 1 : L ψ(0, ξ l ) = 0, l = 1 : L (H(∂ t )u) (t) L l=1 μ l L ψ(t, ξ l ). ( 15
)
More details can be found in [START_REF] Montseny | Représentation Diffusive[END_REF].

STATE FORMULATION OF THE POROUS MEDIUM MODEL

Consider the convolution operators H 1 (∂ t ) and H 2 (∂ t ) with respective symbols

H 1 (p) = 1/(p ρ eff (p)) and H 2 (p) = 1/(p χ eff (p)).
These functions are decreasing at infinity and analytic in C R -. So, by considering the contour γ defined by γ(ξ) = -|ξ| , we can show that operators H 1 (∂ t ) and H 2 (∂ t ) satisfy the hypothesis of theorem 3. After computations as described in section 2, we obtain the following measures μ 1 and μ 2 respectively associated to H 1 (∂ t ) and H 2 (∂ t ):

μ 1 (ξ) = a π ρ √ b ξ -1 ξ 2 + a ξ 2 -a 2 1 ξ>2a + k 1 δ ξ 1 (ξ), μ 2 (ξ) = a c π χ √ b ξ -1 ξ 2 (1 -c) 2 + a 2 ξ -a 2 1 ξ>2a + 1 χ δ 0 (ξ) + k 2 δ ξ 2 (ξ),
where:

ξ 1 = a( √ 17-1) 4 > 0, ξ 2 = a ( √ 1+16(1-c) 2 -1) 4(1-c) 2 > 0, k 1 = √ 17-1 ρ √ 17 > 0, and k 2 = c √ 1+16(1-c) 2 -1 χ(1-c) √ 1+16(1-c) 2 > 0.
Let us express the system (1) under the form:

u = -H 1 (∂ t ) ∂ z P/e P = -H 2 (∂ t ) ∂ z u/e.
From results of section 2, we get the following diffusive formulations of operators ∂ z P → u and ∂ z u → P :

∂ t ψ 1 = -ξ ψ 1 - 1 e ∂ z P u := μ 1 , ψ 1 , and ∂ t ψ 2 = -ξ ψ 2 - 1 e ∂ z u P := μ 2 , ψ 2 .
Then we get the following augmented model, defined on (t, z, x, ξ) ∈ R * + ×]0, 1[×Γ × R + and input-output equivalent to (1,3,4):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t ψ 1 = -ξ ψ 1 - 1 e μ 2 , ∂ z ψ 2 ∂ t ψ 2 = -ξ ψ 2 - 1 e μ 1 , ∂ z ψ 1 μ 1 , ψ 1 (t, 1, x, .) = 0 μ 2 , ψ 2 (t, 0, x, .) = w(t, x) y(t, x) = μ 1 , ψ 1 (t, 0, x, .) . ( 16
)
As μ 1 and μ 2 are positive, the functional

(ψ 1 , ψ 2 ) → (ψ 1 , ψ 2 ) μ,x = [ (μ 1 |ψ 1 | 2 + μ 2 |ψ 2 | 2 )dξdz] 1 2
is a semi-norm. We have:

Proposition 5. The functional E ψ x = 1 2 (ψ 1 , ψ 2 ) 2 μ,
x verifies on any solution (ψ 1 , ψ 2 ) of ( 16) :

dE ψ x dt = - 1 0 ∞ 0 ξ μ 1 |ψ 1 | 2 dξ dz - 1 0 ∞ 0 ξ μ 2 |ψ 2 | 2 dξ dz + 1 e w y. ( 17
)
Proof.

dE ψ x /dt = - i 1 0 ∞ 0 ξ μ i |ψ i | 2 dξ dz - 1 e 1 0 P ∂ z u dz - 1 e 1 0 ∂ z P u dz = - i 1 0 ∞ 0 ξ μ i |ψ i | 2 dξ dz - 1 e P (t, 1, x) u(t, 1, x) + 1 e P (t, 0, x) u(t, 0, x).
The two first terms of the right member of ( 17) are negative: they express the instantaneous dissipation of ( 16) for any given value of x, whereas the last term expresses the instantaneous exchanged power. Moreover, thanks to the property ψ(t, 0) = 0, we deduce from (17) the positiveness of the quadratic form w → Q T (w) := T 0 wy dt for any T > 0 and any x: the passive feature of the absorbent wall is restored by model ( 16).

COUPLING WITH AN AEROACOUSTIC WAVE MODEL

As a consequence of ( 17), the coupling of ( 16) with another passive dynamic model leads to a globally dissipative system provided that the energy functional is chosen such that the transfer between the two subsystems is balanced. This is highlighted in the sequel.

We consider the aeroacoustic model studied in [START_REF] Mazet | Control of Aero-acoustic Propagations with Wall Impedance Boundary Conditions: Application to a Porous Material Model[END_REF]. For simplicity we mainly consider the non convective case with rectangular domain ]0, Z[×]0, X[ and boundary control at z = Z:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t u ρ + c 0 0 ∇ ∇ T 0 u ρ = f u x (t, z, 0) = u x (t, z, X) = u z (t, 0, x) = 0 u z (t, Z, x) = φ y(t, x) w(t, x) = c 0 ρ 0 ρ(t, Z, x), ( 18
)
where u and ρ are the velocity and the density, f is the perturbation source, c 0 is the velocity of wave fronts inside the fluid medium and φ is the porosity coefficient of the porous material. The acoustic energy of (u, ρ) T is classically given by the functional:

E m = 1 2 (u, ρ) T 2 L 2 .
We can define a global energy for (16,18) as:

E = X 0 e(x)E ψ x dx + ρ 0 φ E m .
We have: Theorem 6. The system (16,18) is dissipative in the sense:

dE dt = - X 0 e(x) 1 0 ∞ 0 ξ μ 1 |ψ 1 | 2 dξ dz dx - X 0 e(x) 1 0 ∞ 0 ξ μ 2 |ψ 2 | 2 dξ dz dx 0
on any solution (u, ρ, ψ) T of (16, 18).

Proof. From ( 18) and:

dE m /dt = - X 0 1 0 c 0 0 ∇ ∇ T 0 u ρ • u ρ dz dx = -c 0 X 0 1 0 (∂ x ρu x + ∂ z ρu z + ∂ x u x ρ + ∂ z u z ρ) dz dx = -c 0 ( 1 0 [ρu x ] X x=0 dz + X 0 [ρu z ] Z z=0 dx) = -c 0 X 0 ρ(t, Z, x) u z (t, Z, x) dx.
Note that this global energy dissipation due to wave absorption at the boundary, is not accessible under the frequential or convolution forms (1) and ( 5). In this sense, the diffusive model ( 16) is conform to the physical interpretation of such media.

In the convective case [START_REF] Mazet | Control of Aero-acoustic Propagations with Wall Impedance Boundary Conditions: Application to a Porous Material Model[END_REF]), similar results can be obtained up to suitable technical adaptations.

DIFFUSIVE REALIZATION OF THE IMPEDANCE OPERATOR

In order to build state realizations of the boundary control defined by (1,3,4), it can be judicious to start from the symbol of the input-output operator w → y, denoted Q(∂ t ) in the sequel. From the numerical point of view, the realization of this operator is more accurate and cheaper than the simulation of the wave propagation inside the medium and so is more adapted for example to control purposes. Let us now compute the symbol of Q(∂ t ).

The equations of the porous medium (1,3) can be written:

∂ z X = -e(x) iω AX X(0, ω) = [v, 0] T , ( 19
)
with:

X = P û , A = 0 ρ eff (iω) χ eff (iω) 0 .
Diagonalization of A leads to:

A = M √ χ eff ρ eff 0 0 - √ χ eff ρ eff D M -1 ,
where M is the transition matrix; using the change of unknowns Y = M -1 X, we get the diagonal system of ODEs: ∂ z Y = -e iω DY with limit conditions deduced from (19). After simple computations, we obtain: (iω) .

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ P (ω, z) = P (ω, 0) A(ω)/B(ω, z) + B(ω, z) A(ω) + 1 û(ω, z) = P (ω, 0) χ eff (iω) ρ eff (iω) A(ω)/B(ω, z) -B(ω, z) A(ω) + 1 , (20) where A(ω) = e 2eiω √ χ eff (iω) ρ eff (iω) , B(ω, z) = e -eiωz √ χ eff (iω) ρ eff
Finally, the input-output correspondence w → y = Q(iω) w (which summarizes the frequency behavior of the porous medium) can be deduced by taking z = 0 and replacing iω by p ∈ C in (20):

Q(p) = χ eff (p) ρ eff (p) tanh e(x) p χ eff (p) ρ eff (p) .
It can be shown that the symbol

Q(p) = Q(p) p
satisfies the hypothesis of theorem 3. From section 2, we then get the following state-representation of Q(∂ t ):

∂ t ψ = -γψ + w y = ν, ∂ t ψ = -γ ν ψ dξ + ν dξ w, ( 21
)
where ν is the measure associated to Q(∂ t ) and γ is a suitable contour.

The presence of complex poles of Q(p) (see fig. 1) does not allow us to use the same contour γ as in section 3 (note in particular that

Q(p) ∼ ∞ K1 p tanh(K 2 p)).
For the realization of Q(∂ t ), we so consider a contour γ of the form:

γ(ξ) = |ξ| e i sign(ξ)( π 2 +α) , with α ∈]0, π
2 ] (α will be taken equal to 5π 180 for the numerical results of section 6).

SOME NUMERICAL RESULTS

As described in section 2, we can perform converging approximations of ( 16) and ( 21) by using a discretization (ξ l ) l=1,L of ξ and standard quadratures. We then get approximate dynamic realizations of the form:

Ẋ = A X + Bw y = CX with X(t) ∈ R L ,
such that y y in a suitable sense [START_REF] Montseny | Représentation Diffusive[END_REF]).

The frequency responses of the approximations of Note that by increasing the number of ξ l , the oscillations of the frequency response of Q(∂ t ) will be better approximated.

H 1 (∂ t ), H 2 (∂ t ) and Q(∂ t )
For illustration, the evolution of P obtained from simulation of ( 16) is given in figure 5. The explicit scheme used for this simulation and whose stability is proved in Casenave [2008], is written:

(∂ t ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ψ n+1 1 (x, ξ l ) = a l1 ψ n-1 1 (x, ξ l ) -b l1 G 21 j b j2 ψ n 2 (x, ξ j ) ψ n+1 2 (x, ξ l ) = a l2 ψ n-1 2 (x, ξ l ) -b l2 G 12 j b j1 ψ n 1 (x, ξ j ) u n+1 (x) = l b l1 ψ n 1 (x, ξ l ) P n+1 (x) = l b l2 ψ n 2 (x, ξ l ), (22) with u n+1 (x), P n+1 (x) ∈ R K and ψ n+1 1 (x, ξ l ), ψ n+1 2 (x, ξ l ) ∈ C K ,
where K = 300 is the number of points of discretization in z, and:

a li = e -ξ l 2Δt , b jk = c jk e -ξ j 2Δt -1 -ξ j , G 21 = G 12 = 1 2Δz ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 -1 0 1 . . . . . . . . . -1 0 1 -1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ,
where c li are some coefficients computed by simple quadrature of μ i (ξ l )Λ l (ξ)dξ, Λ l being the classical interpolation function:

Λ l (ξ) = ξ -ξ l-1 ξ l -ξ l-1 1 [ξ l-1 ,ξ l ] (ξ) + ξ l+1 -ξ ξ l+1 -ξ l 1 ]ξ l ,ξ l+1 ] (ξ).
The boundary conditions are: We clearly observe the dissipative and dispersive nature of this propagative model, due to the convolution operators H 1 (∂ t ) and H 2 (∂ t ).

In figure 6 we can see at a fixed time t = 0.09 ms, the functions ψ 1 involved in the synthesis of u. Fig. 6. Functions ψ 1 (t, ., x, ξ l ) l = 1 : L at time t = 0.09 ms

We compare in figure 7 the outputs y := u |z=0 obtained by simulation of ( 16) and ( 21) with the input w given by (24). Note that the difference between the two results is mainly due to the error approximation of z-discretization.

As the simulation of (21) doesn't need any z-discretization, we can consider a high number of ξ l to approximate the operator Q. Then the simulation of ( 21) is more accurate but remains cheaper than the scheme (22) (whose global dimension is 10500).

CONCLUSION

On the basis of the above energy dissipation analysis, we can conclude that the model (1) can be viewed as a dissipative absorbing controller for aeroacoustic waves. Furthermore the quality of the numerical results of the simulation of the boundary operator w → y confirms that the state-realizations ( 16) and ( 21) can be used for more deepened investigations on absorbing boundary control.

Fig. 1 .

 1 Fig. 1. Poles of Q(p) (e = 1)

  respectively obtained with L = 15, 20 and 150 are given in figures 2, 3, 4. In the case of H 1 (∂ t ) and H 2 (∂ t ) (respectively Q(∂ t )), 4 decades from 10 2 to 10 6 (respectively 3 decades from 10 2 to 10 5 ) are covered by the ξ-discretization. The numerical values of parameters are: Λ = Λ = 0.1 10 -3 m, ρ 0 = 1.2 kg.m -3 , P 0 = 10 5 Pa μ = 1.8 10 -5 kg.m -1 .s -1 , γ = 1.4, α ∞ = 1.3, e = 2 10 -2 m.

Fig. 5 .Fig. 2 .

 52 Fig. 5. Evolution of P = l b l2 ψ 2 (ξ l ) (N.B: the unit of length for the z-axis is 10 -2 m)

Fig. 3 .

 3 Fig. 3. Frequency response of operator H 2 (∂ t )

Fig. 4 .

 4 Fig. 4. Frequency response of operator Q(∂ t ) u(t, 1) = 0, (23)P (t, 0) = 1cos(2π × 1.3 10 4 t) 1 [0,T ] (t) = w. (24)

Fig. 7 .

 7 Fig. 7. Outputs y obtained by simulation of (16) (--) and (21) (-).