Christophe Rouzaud

Fabrice Gatuingt

Guillaume Hervé

Olivier Dorival

A new multi-frequency approach based on Padé approximants for the treatment of transient dynamics problems with the variational theory of complex rays

Keywords: model reduction, frequency resolution, Padé approximants, VTCR

come

A new multi-frequency approach based on Padé approximants for the treatment of transient dynamics problems with the variational theory of complex rays

B

vector form of L X amplitudes vector of the wave functions T truncation order χ coefficient vector N number of degrees of freedom, or the number of wave functions used L, M truncated orders P L , Q M Taylor series truncated respectively to order L and M p, q Taylor series coefficients x coefficients for the Taylor serie X A matrix of Padé coefficients

Introduction

When a projectile impacts a structure, several wave regimes occur that depend on the observation time and the structure properties and geometry. Classically these wave regimes can be divided into two types: tension-compression waves and flexural waves which propagate longitudinally. The first type is non-dispersive and mainly affects the damages in the impacted area, while the second one is dispersive and has a global effect on the structure. The vibratory response of elastic structures in the low-frequency range is not a problem even for complex structures and for the high-frequency range, it exists efficient numerical tools mainly based on energetic approaches.

In this work we are interested in the effect of an aircraft impact on large civil engineering industrial structures. This type of load covers a large frequency range, including the medium frequency range [START_REF] Hervé | On numerical implementation of a coupled rate dependent damage-plasticity constitutive model for concrete in application to high-rate dynamics[END_REF][START_REF] Hervé | Optimizing the analysis of airplane crash induced spectra by means of generic airplane methodology[END_REF]. The cut-off frequency for this type of load is typically within the 50-100 Hz range, which will be referred to as the medium frequency range for civil engineering structures [START_REF] Wolf | Response of equipment subjected to aircraft impact[END_REF]. Due to the complexity of the problem, it is impossible to use simplified approaches, as in [START_REF] Jacquelin | A modelling of an impacted structure based on constraint modes[END_REF], based on spring-mass more or less complex models. The finite elements method with a global model seems to be the most appropriate choice in our case. In order to properly describe the response of the impact area, the finite element method requires that space and time discretizations are sufficiently refined. Ten linear elements per wavelength if often used as a rule of thumb [START_REF] Babuška | A generalized finite element method for solving the helmholtz equation in two dimensions with minimal pollution[END_REF][START_REF] Marburg | Six boundary elements per wavelength: is that enough?[END_REF]. Time step is limited by the Courant condition for explicit time integration scheme used for this kind of simulations, so that the duration for a wave crossing a discrete spatial grid is less than the time for the wave to travel to adjacent grid points. For this reason, computation time quickly becomes prohibitive when entire civil engineering structures need to be calculated under impact loadings. As a consequence of the mesh refinement limitation, the medium frequency range is often ignored in this type of simulation, which alters the loading that the structure undergoes. This may result in wrong predictions from the simulations.

We can find in the literature some methods that can deal with complex structures. They can be based on structural partitioning as in [START_REF] Kassem | Structural partitioning of complex structures in the mediumfrequency range. an application to an automotive vehicle[END_REF], on methods developed for some types of structures [START_REF] Geng | Reconstruction of transient vibration and sound radiation of an impacted plate using time domain plane wave superposition method[END_REF], on specific finite-element developments [START_REF] Zhu | A hybrid finite element-energy finite element method for mid-frequency vibrations of built-up structures under multi-distributed loadings[END_REF]. In a previous work [START_REF] Rouzaud | A new way for the simulation of the impact on reinforced concrete structures[END_REF] we proposed a new way of calculating the shaking of civil engineering structures, that allows to account for medium frequency range. This methodology is based on the use of the Variational Theory of Complex Rays (VTCR) [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF][START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF][START_REF] Rouch | The variational theory of complex rays: a predictive tool for mediumfrequency vibrations[END_REF][START_REF] Rouch | Vibrational analysis of structures with stochastic interfaces in the mediumfrequency range: Experimental validation on a touch screen[END_REF], which is a wave-based computational approach dedicated to the resolution of forced vibration problems at a given frequency. The VTCR uses a weak formulation of the boundary conditions that is projected on a space of admissible shape functions that satisfy both the governing equation and the constitutive relation. Such shape functions are found to show a strong vibrational content as they describe classical waves in infinite or semi-infinite domains. In order to address transient dynamic problems, the proposed strategy is to use a shift from the time domain to the frequency domain that is achieved by a Fast Fourier Transform (FFT) technique [START_REF] Ladevèze | A new computational method for transient dynamics including the low-and the medium-frequency ranges[END_REF]. After solving the problem in the frequency domain with the VTCR, a time recombination is performed by an inverse Fast Fourier Transform (IFFT). The approach used in [START_REF] Rouzaud | A new way for the simulation of the impact on reinforced concrete structures[END_REF] can be very efficient because at each given frequency in the loading spectrum, the VTCR performs a fast resolution of the problem. However, when the loading spectrum contains a lot of frequencies, performances of the frequency-by-frequency treatment are decreased dramatically. The same kind of problem is found when the solution of the Helmholtz equation is required over a broad frequency range that can be found in various engineering applications.

The following article proposes an alternate multi-frequency approach based on Padé approximants for improving the reconstruction of the Frequency Response Functions (FRF). Padé approximants technique consists in the approximation of a function by a rational function of given order. The tech-nique was developed around 1890 by Henri Padé [START_REF] Padé | Sur la représentation approchée d'une fonction par des fractions rationnelles[END_REF]. The model reduction by Padé 's approximants we have developed is based on its use for acoustic problems given in [START_REF] Rumpler | A finite element approach combining a reduced-order system, padé approximants, and an adaptive frequency windowing for fast multi-frequency solution of poro-acoustic problems[END_REF] Using reduced-order models build via Padé 's approximations are now well established for systems exhibiting polynomial frequency dependency of second order kind and for independent frequency excitations. This paper deals with a more complex situation where a more complicated wave number dependency is encountered due to the frequency dependent loading.

The paper is structured as follows: Section 2 presents the reference problem and its solution using the VTCR; Section 3 gives a general presentation of the reconstruction method by Padé approximants; Section 4 illustrates an application of the strategy on a simple example and on an industrial structure. Finally, conclusions and perspectives are drawn in Section 5.

Description of the Variational Theory of Complex Rays for transient dynamics

In this section, a brief reminder of the time-frequency strategy (see Figure 1) and of the VTCR methodology are given. The Variational Theory of Complex Rays (VTCR) [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF][START_REF] Ladevèze | A new computational method for transient dynamics including the low-and the medium-frequency ranges[END_REF] allows to handle low and medium frequency ranges. High frequencies, which are known to vanish quickly and/or which do not contain a lot of energy, are neglected [START_REF] Vanmaele | An efficient wave based prediction technique for plate bending vibrations[END_REF].

The reference problem at a given frequency for an assembly of two substructures

Let us consider the case of homogeneous Kirchhoff-Love's thin plates Ω i which vibrate at a fixed frequency f = ω 2π . The thickness is denoted by h i and the density ρ i . Under the Kirchhoff-Love assumptions, the displacement U i of the average surface becomes:

U i (x, y, z) = u i (x, y) + w i (x, y) e 3i + z θ i θ i (x, y) = -grad w i (x, y) (1)
where u i is the displacement of the average surface and w i is the out-of-plane displacement. The average surface of the plate can be defined by two independent parameters α i and β i . The location of a point on the medium surface is defined by the position vector r i (α i , β i ) (see Figure 2). The tangent local plane to the plate at point r i (α i , β i ) is defined by two vectors, ∂r i ∂α i = A i e αi and ∂r i

∂β i = B i e β i (the vectors e αi and e β i are unit standards). The vector e 3i is defined by e 3i = e αi ∧ e β i .

The base e αi , e β i , e 3i is then orthogonal. Let us consider a domain Ω composed of two plates Ω i and Ω j with a common border Γ ij . The actions of the environment are modeled on Ω i/j by displacement prescribed on ∂ U Ω i/j and prescribed stresses on ∂ F Ω i/j .

Figure 3 illustrates these actions on an assembly of plates Ω i and Ω j .

Frequency (Hz) In elasticity , K CP i are Hooke's operators related to each area, under plane stress assumptions; they are defined as:

K CP i = (1 + iη i )K 0i CP = (1 + iη i )     E αi 1-ν αi ν βi ν αi E αi 1-ν αi ν βi 0 ν βi E βi 1-ν αi ν βi E βi 1-ν αi ν βi 0 0 0 √ E αi E βi 2(1+ √ ν αi ν βi )     (eαi,eβ i ,e 3i) (2)
where i is the imaginary unit, E α,βi is the Young's modulus, ν α,βi are the Poisson's ratio in e αi and e β i directions, h i is the plate thickness, ρ i is the density and η i is the structural damping coefficient of each sub-domain i.

The variational formulation associated with the VTCR

The first ingredient of the VTCR is a global weak formulation of the boundary conditions in terms of both displacements and forces. The weak formulation can be expressed as follows: find (U i , σ i ) ∈ S ad,i such as:

A s i s j , δs i δs j = L δs i δs j (3)
with the following general forms:

A s i s j , δs i δs j =Re -iω ∂ U Ω i δσ i .n i .U * i dS + ∂ F Ω i σ i .n 1 .δU * i dS + ∂ U Ω j δσ j .n j .U * j dS + ∂ F Ω j σ j .n j .δU * j dS + 1 2 Γ ij (δσ i .n i + σ j .n j ) .(U i + U j ) * +(σ i .n i -σ j .n j ) .(δU i -δU j ) * dS (4) L δs i δs j =Re -iω ∂ U Ω i δσ i .n i .U * id dS + ∂ F Ω i F id .δU * i dS + ∂ U Ω j δσ j .n j .U * jd dS + ∂ F Ω j F jd .δU * j dS (5)
where:

• the integral part on ∂ U Ω i/j satisfies the prescribed displacements on Ω i/j ,

• the integral part on ∂ F Ω i/j satisfies the prescribed stresses on Ω i/j ,

• the integral part on Γ ij satisfies the transmission conditions on the boundary Γ ij .

Re is the real part of a complexe quantity and * the conjugate part. Spaces S ad,i/j are the admissible fields associated with homogeneous conditions on the substructure Ω i/j .

The previous weak equation is projected on an approximation basis subspace S ad,i/j that is composed of shape functions that satisfies exactly both the constitutive relation and the dynamic equilibrium equation.

Derivations of two-scale shape functions and the discretized problem

The VTCR can be viewed as a Trefftz's method which particularity is the use of Herglotz wave functions to represent the vibrational field. Those functions involve an integral combination of plane waves in all the directions. The solution is assumed to be properly described locally as the superposition of an infinite number of these functions which can be written in the following manner:

U i (Z i , P i ) = P i ∈C i U ni (P i ) .e P i .Z i dP i on Ω i σ i (Z i , P i ) = P i ∈C i C ni (P i ) .e P i .Z i dP i on Ω i (6)
where Z i denotes the position vector. U ni and C ni contain the amplitudes of the associated local vibration waves. P i characterizes the direction of local waves and is defined on circle C i .

In order for these local waves (U i , σ i ) to be admissible, they must be in S ad,i and satisfy the constitutive and dynamic equilibrium equations. Combining these two equations leads to some remarkable shape functions (named complex rays) characterized by various expressions of P i . For example in the case of a flexural plate the vibration waves P are defined as :

(P .P) 4 = 12(1 -ν 2 )ρω 2 (1 + iη)Eh 2 (P .P) 2 (7)
Note that the term e P i .Z i describes the short scale variation, which is known in a closed form, while U ni and C ni are the long scale variations, which remain the unknown to be determined. Let us note that admissible space S ad,i is of infinite dimension since all the directions of propagation P i are taken into account. To end up with a finite dimension problem that can be solved numerically, one needs to discretize S ad,i into a finite dimension space S h ad,i . The displacement of any point of the substructure is generated by a basis of admissible complex waves. The unknowns are the generalized amplitudes U h i (P i ) of the basis on Ω i . Accounting for all the directions ϕ i in C i leads to an integral over C i . This integral takes the following form:

U h i (z i ) = ϕ i ∈C i U n h i (ϕ i )e P i (ϕ i ).z i dϕ i (8)
σ h i (z i ) = ϕ i ∈C i C n h i (ϕ i )e P i (ϕ i ).z i dϕ i (9)
The integrals in Equation ( 8) and Equation ( 9) can be discretized as Dirac's functions and one can consider the approximate amplitude U h (P i (ϕ i )). The angular distributions of the plane waves for all points in the substructure are assumed to be well-described by this discontinuous angular distribution. The choice of the angular discretization and therefore the number of waves required for solving this kind of problem is related to the number of waves in the structure and on each edge [START_REF] Kovalesky | Optimized approximation space for trefftz-discontinuous galerkin methods, application to the variational theory of complex rays[END_REF].

The number of wavelength n i depends on the celerity of the waves: for pressure waves,

c p = E α/βi ρ i (1-ν αi ν βi ) , for shear waves, c sh = √ E αi E βi 2ρ i (1+ √ ν αi ν βi ) and for bending waves c sv = √ ω 4 E α/βi h 3 i 12ρ i h i (1-ν αi ν βi )
. Finally at a given frequency ω, the VTCR equation results in a linear problem of the following form:

A(ω)X(ω) = B(ω) (10)
where A and B come from the discretization of operators A (Eq. 4) and L (Eq. 5). X contains the unknown slow scale amplitude of the shape functions. Due to the strong vibrational sense of the shape functions, the VTCR offers a substantial reduction of the size of the problem to solve compared to finite element methods.

A new approach for broadband analysis by the VTCR

The treatment of an aircraft crash on a civil engineering structure by the strategy proposed in [START_REF] Rouzaud | A new way for the simulation of the impact on reinforced concrete structures[END_REF] results in a frequency-by-frequency resolution by the VTCR before recombining the solution in the time domain. This may degrade the strategy performances when lots of frequencies need to be accounted for. This article provides an alternate approach based on Padé approximants for the rapid reconstruction of frequency response functions. The aim is to reduce the number of frequencies at which a complete problem must be solved, and then to interpolate the solution between these frequencies. To this aim, a reconstruction based on Padé approximants allows a good complement, involving the calculation of the complete solution and its derivatives for a limited number of frequencies, and a reconstruction around these points. This approach is helped by the expression of the frequency dependency in the form of scalar functions in the overall VTCR equation.

Padé approximants

All VTCR problems, with harmonic vibrations f = ω 2π , have the following form:

A(ω)X(ω) = B(ω) (11)
In a scheme of direct solution, the calculation of the complete solution of Equation ( 11) on only few frequency points, coupled to a reconstruction of an approximated solution between them, was found to be an effective approach [START_REF] Avery | Fast frequency sweep computations using a multi-point padébased reconstruction method and an efficient iterative solver[END_REF]. Thus, given a function X(ω) of ω, the Taylor series development around ω 0 gives the following approximate solution:

X(ω 0 + ∆ω) ≈ X T (∆ω) = T i=0 χ i (∆ω) i (12)
T is the truncation order of the Taylor series expansion. The coefficient vector χ i is given by:

χ i = X (i) (ω 0 ) i! (13)
where

X (i) (ω 0 ) is the i th derivative of X over ω, taking ω = ω 0 and considering that X (0) (ω 0 ) = X(ω 0 ).
The determination of the first T derivatives from X is then required to estimate an approximate solution around ω 0 . However, in the context of the series expansion, it has been shown that, for a function with poles as it is the case for structural FRFs, an expansion at a ratio of two power series generally performs a better representation (greater convergence range and better convergence rate) compared to the Taylor series for a given truncation order [START_REF] Orszag | Advanced mathematical methods for scientists and engineers[END_REF][START_REF] Baker | Padé approximants[END_REF][START_REF] Cochelin | Méthode asymptotique numérique: une technique de résolution des équations non linéaires[END_REF]. In the next part, the approximation at a given frequency ω 0 is presented under the methodology for determining Padé approximants of Equation [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]. For the sake of clarity, let us consider the j th component of X(ω 0 ), which is sought around ω 0 as :

X j (ω 0 + ∆ω) ≈ P Lj (∆ω) Q M j (∆ω) (14)
where P Lj (∆ω) and Q M j (∆ω) are two Taylor series, truncated respectively to order L and M and defined as:

P Lj (∆ω) = L i=0 p i j (∆ω) i Q M j (∆ω) = M i=0 q i j (∆ω) i (15)
The approximation of X j (ω 0 + ∆ω) is given by the single determination of the coefficients p i j and q i j . An efficient approach [START_REF] Baker | Padé approximants[END_REF] to determine these coefficients is to identify the coefficients of a Taylor series expansion X L+M j (∆ω) to the order L + M (Equation 12), and so solving the following linear system:

P Lj (∆ω) -X L+M j (∆ω)Q M j (∆ω) = 0 (16)
where the coefficients of order equal to ∆ω are extracted to form a set of (L + M + 1) equations. With the form of Equation ( 16), the linear system of equations allows zero for Q M j (∆ω), which makes it more suitable to account for a resonance frequency in the original answer. However, the system of underdetermined Equations ( 16) gives a solution of (L + M + 2) coefficients up to a multiplicative constant. Therefore, they are usually normalized so that the denominator coefficient of the zero order q 0 j is taken equal to 1 [START_REF] Baker | Padé approximants[END_REF]. In the following, all the equations can then be put as:

p 0 j = x 0j p 1 j -x 0j q 1 j = x 1j . . . p L j -x L-1 j q 1 j -• • • -x 0j q L j = x Lj -x Lj q 1 j -x L-1 j q 2 j -• • • -x L-M +1 j q M j = x L+1 j . . . -x L+M -1 j q 1 j -x L+M -2 j q 2 j -• • • -x Lj q M j = x L+M j                      (17)
with

x ij = 0 if i < 0 q i j = 0 if i > M (18)
Equation ( 13), gives the expression of the coefficients x ij . Equations ( 17) can be explicitly written by taking into account the (L + M ) first derivatives of X j (ω) to ω 0 :

i!p i j - i k=1 i k k!X (i-k) j (ω 0 )q k j = X (i) j (ω 0 ), for i = 0, ..., L + M (19)
where

p i j = 0 if i > L q k j = 0 if k > M (20)
and where the binomial coefficient are given by:

i k = i! k!(i -k)! (21)
This system can be solved efficiently in the following matrix form for each DoF j involving derivatives up to L + M order of the solution vector :

[A] j          p 0 . . . p L q 1 . . . q M          j =          X(ω 0 ) . . . X (L) (ω 0 ) X (L+1) (ω 0 ) . . . X (L+M ) (ω 0 )          j , for j = 1, ..., N (22)
Based on the solution in a given frequency, X(ω 0 ), this method is defined in the literature as the Padé approximation in a single point. It can easily be extended to several points, as shown for example in [START_REF] Avery | Fast frequency sweep computations using a multi-point padébased reconstruction method and an efficient iterative solver[END_REF]. The goal is to estimate the combination of this approach coupled with VTCR. The use of a multipoint Padé approximation is considered as a natural extension of this work. Similarly, an approximation of the solution using the Taylor series expansion, solving Equations ( 22) requires the determination of the (L + M ) th derivative from the solution vector X(ω) in ω 0 . This is discussed in the next section.

Successive derivatives of the solution vector

The (L + M ) successive derivatives of X to the pulsation ω can be obtained by Equation [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]. The successive derivatives can be computed by a recursive formula, leading to the solution of a linear system with multiple second members, which is very efficient from a numerical point of view : only one decomposition of the coefficient matrix is to be done. To the order i, the expression becomes:

i k=0 i k A (i-k) (ω 0 )X (k) (ω 0 ) = B (i) (ω 0 ), for i = 1, ..., (L + M ) (23)
where the zero order derivatives correspond to the undifferentiated functions. Extracting the term of highest order of the sum in Equation ( 23) leads to the following recursive expression for the i th order derivative of X in ω 0 :

A(ω 0 )X (i) (ω 0 ) = B (i) (ω 0 ) - i-1 k=0 i k A (i-k) (ω 0 )X (k) (ω 0 ), for i = 1, ..., (L + M ) (24)
This implies that the successive derivatives of X required for the determination of the Padé approximation can be computed efficiently as the solution of a N dimension set of equations. Indeed, the solution can be performed by direct methods using a decomposition of the system matrix. Factorising is then the most numerical consuming step but is processed only one time. Regarding the multiple secondary vectors, they are built from derivatives of the system matrix and derived from lower order vector solution.

Procedure for the reconstruction of the frequency response using the Padé approximants

This section summarizes the calculation steps to approximate the solution over a frequency range using Padé approximants. They can be separated into two parts:

• first, at a given frequency ω 0 the solution and its derivatives to successive orders (L + M ) are calculated according to Equation [START_REF] Djellouli | A fast method for solving acoustic scattering problems in frequency bands[END_REF]. This implies the solution of (L + M + 1) problems of size N for each frequency at which the solution must be approached;

• then the solution around ω 0 is reconstructed, involving solving N problems of (L + M + 1) dimensions (see Equation ( 22)) to determine the Padé approximants for each degree of freedom. Following this, N evaluations of the rational fraction (Equation ( 14)) are made for each ∆ω, where the solution evaluated is approached.

Given a discretization of the frequency domain into intervals associated with their center frequency ω 0 , the detailed procedure corresponding to a frequency range is presented in the algorithm given in Table 1.

Table 1: Reconstruction algorithm of the response over a frequency range using Padé approximants.

Algorithm of multi-frequency solution using Padé approximants 1. Calculate the matrix and the second member in ω 0 from Equation (11) 2. Solve Equation [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF] for the solution X(ω 0 ) 3. For i = 1 to L + M 4.

Solve Equation ( 24) for X (i) (ω 0 ) 5. end 6. For j = 1 to N 7.

Solve Equation ( 22) for the Padé approximants 8.

Evaluate the multi-frequency solution approximated around ω 0 (Equation ( 14)) 9. end Three points can be considered in view of the given procedure. The first point is to select the secondary frequencies where the solution will be evaluated, or an appropriate discretization of the frequency range. Secondly, for each coarse frequency, the order of the polynomial expansions and frequency range around ω 0 should be selected consequently, which raises the question of values of L and M . Third, one must evaluate the efficiency conditions for this approximation scheme.

Although the second and third points are discussed in details in [START_REF] Avery | Fast frequency sweep computations using a multi-point padébased reconstruction method and an efficient iterative solver[END_REF] for structural or acoustic applications, the choice of secondary frequencies has not received much attention in the state of art. A suggestion is made in [START_REF] Rumpler | A finite element approach combining a reduced-order system, padé approximants, and an adaptive frequency windowing for fast multi-frequency solution of poro-acoustic problems[END_REF]. Regarding the order of polynomial expansions to consider for both the numerator (L) and the denominator (M ) in Equation ( 14), the following constraint proposed in [START_REF] Avery | Fast frequency sweep computations using a multi-point padébased reconstruction method and an efficient iterative solver[END_REF] reduces the range of possibilities giving satisfying approximations:

M = L + 1 (25)
In addition, there is an upper limit for the maximum order which can be defined, due to ill conditionning in the resulting matrix system from Equations [START_REF] Baker | Padé approximants[END_REF]. Therefore, depending on the frequency discretization chosen for a scan of the frequency range, the Padé approximation should be more effective if it could converge on several frequency increments around ω max . The time allocated to solve several different linear systems of size N leads us to solve (L+M +1) systems with the second member of N equations (1 decomposition, and (L + M + 1) solutions) and N systems (L + M + 1) equations, with (L + M + 1)

N . This methodology has been extremely effective when applied to important structural and acoustic problems with single field and a a priori discretization of the frequency range [START_REF] Djellouli | A fast method for solving acoustic scattering problems in frequency bands[END_REF][START_REF] Avery | Fast frequency sweep computations using a multi-point padébased reconstruction method and an efficient iterative solver[END_REF].

Application of multi-frequency analysis

Multi-frequency analysis on a simple case study

In this study to compare the multi-frequency analysis via Padé approximants with a classical frequencyby-frequency solution, we set ourselves in the case of a bending problem applied to a single Kirchhoff-Love plate (see Figure 5). The case presented is a plate clamped on one of its edges, free on two others and loaded by a out-ofplane line force on its last edge. The plate is square and has the following dimensions: 1m × 1m × 0.7mm. The loaded edge is thus subjected to a unit line force F d with a direction perpendicular to the mean surface of the plate. The frequency range is [2π × 0 Hz ; 2π × 350 Hz]. Materials and structural parameters are given in Table 2. The comparison between a finite element solutions calculated with Cast3m [START_REF]Cast3m[END_REF] and the VTCR solution at ω = 2π × 150 Hz is presented in Table 3. To study the efficiency of multi-frequency solution using Padé approximants, the flexural plate described is studied for a 40 Hz medium frequency range, from 290 to 330 Hz. In this first study, the goal is to investigate the efficiency of Padé approximants applied to the VTCR problem in terms of both results and calculation time compared to a classical solution at each frequency. This solution repeated each 1 Hz will be considered as the reference solution. Figure 6 shows the convergence of the frequency response function (FRF) for several sets of parameters L and M . Here a single analysis is calculated with 310 Hz as a center frequency. One can observe that increasing the order L and M of Padé approximants quickly leads to the convergence to the reference solution. Application of Padé approximants to VTCR is found to be efficient when the approximants order is high enough to describe the solution over the frequency range. Nevertheless, if this analysis does not require one inversion of the bilinear matrix for the central frequency, it requests to calculate at each frequency the successive derivatives of this matrix and the second member. One can easily imagine that these derivatives may involve significant computational cost if they are not written recursively. This cost also increases with orders L and M as the order of derivation is equal to L + M . Figure 7 compares the computation time for multi-frequency analysis with two sets of parameters and the classical analysis. Note that the approach using Padé approximant is particularly interesting with an important frequency range. Indeed, in the case of a frequency-by-frequency solution, the overall calculation time for 40 Hz is 95 s, and includes assembly and solution times required for each frequency. For the multi-frequency analysis, the computation time is 47 s for orders L = 1 and M = 2, and 65 s for orders L = 5 and M = 6. In these two cases, a high proportion of the calculation time is dedicated to the construction of the successive derivatives.

Indeed, in this study, the matrix to invert is rather of small size as we treat only one subdomain. Thus, the matrix inversion in the case of a problem with a few plate structures is numerically inexpensive. However, for problems involving a larger number of subdomains, the inversion time will be much higher, and so the analysis via Padé approximants can be expected to be more efficient. From this first example one can therefore conclude that the choice of the discretization analysis is dependent on the frequency content that one seeks to represent. To solve this, an analysis where the frequency ranges are define by dichotomy from an error criterion can be used.

Numerical application of the method on a civil engineering structure

In this section the VTCR coupled with Padé's multi-frequency approach is used to calculate the frequency response of a structure subjected to a time versine loading. Experimental shaking or vibration studies on large civil engineering structures are lacking in actual literature. For that reason, the comparisons between the VTCR-Padé strategy and the classical explicit time finite element approach is here based on international rules [START_REF]Office for Nuclear Regulation, Safety assessment principles (saps)[END_REF] for the design of nuclear civil engineering structures. The comparison between simple VTCR (frequency-by-frequency analysis) and FEM calculation for this structure is given in [START_REF] Rouzaud | A new way for the simulation of the impact on reinforced concrete structures[END_REF].

The load applied on the civil engineering structure sketched in Figure 10 may be seen as an impulsive load. First the discrete Fourier transform of the load is calculated. Then the frequency response at two chosen points (P 1) and (P 2) of the structure (see in Figure 10) is given by the VTCR.

A concrete structure with a thickness of 1.0 m is considered where the mechanical properties of concrete are calculated according to the Eurocode 2 rules:

• concrete B30 = 30 MPa, • Young's modulus = 34 GPa, • Poisson's ratio = 0.2,

• mass density = 2500 kg/m 3 ,

• damping coefficient = 0.04 (a hysteretic damping is used),

The structure is subjected to a time loading applied at the center of a side wall (point P 0 in Figure 10). This loading and the boundary conditions of the structure are modeled and described in Figure 10. The versine force F P 0 (t) with an amplitude of 1 MN is defined in Figure 11. Frequency (Hz)

1.E+00

1.E-01

1.E-02

1.E-03

1.E-04

1.E-05 Two hundred waves are sufficient to properly represent the frequency response [START_REF] Rouzaud | Ébranlement de structures en béton armé soumises à un phénomène transitoire (shaking of reinforced concrete structures submitted to transient dynamic phenomenon[END_REF]. Table 4 shows the VTCR frequency solution obtained at four different frequencies. A verification of the boundary conditions is done to check the quality of the solution. 1.E-01

1.E-02

1.E-03

1.E-04

1.E-05

1.E-06

1.E-07

1.E-08

Figure 13: FRF over an adapted wide frequency bandwidth, with a standard criterion of 10 %, for point P 1.

Figure 14 shows the results for the point P 2 with the same previous bandwidth discretization. We can observe that the error can be higher than 10 % on some frequencies. These errors are mainly due to the fact that the modal content of the floor at the point P 2 is slightly different than that of the area of the point P 1. Nevertheless, The error remains below 20 % which is better than finite element calculations which one not refined enough. Frequency area with an error ≥10%

1.E-02

1.E-03

1.E-04

1.E-05

1.E-06

1.E-07

1.E-08

1.E-09

Figure 14: FRF over the previous adapted wide frequency band for point P 2.

To go further in the analysis, it is possible to return to the time domain in order to analyze the impact of these errors on the response of the structure. Figure 15 shows the displacement at point P 2 of the floor. One can thus see that despite the 20 % maximum error at some frequencies, the time displacement response is only slightly affected. The error in this example is less than 1 mm. The displacement response allows us to observe easily and graphically the low frequency content. It is recommended to use the acceleration response to focus on errors in medium frequencies ranges. So Figure 16 allows us to qualify errors with multi-frequency analysis in medium frequencies. One can see that the approximation error remains low (≤ 4 g in our case). To conclude, in this application a multi-frequency analysis with Padé's approximants provides good results due to the fact that the modal content of the studied structure is known. In this case, the frequency bands could be adapted to the modal content. Nevertheless, it is interesting to note that in a first approach, the proposed adapted multi-frequency analysis allows to describe an approximate response of the structure with a computational cost reduced by approximately 30 %. In the case studied, the computational cost for a frequency-by-frequency resolution is around 105 min.

If the modal content is not known, the same analysis can be done by choosing first arbitrary frequences. In a second step the bandwidth can be divided by two and an error between the two solutions computed. If this error seems not acceptable the bandwidth can be divided by two again and the news error computed. This approach can be done until the differences between two solutions is assumed acceptable.

Conclusions

In this paper, a new way has been presented for the treatment of impact problems on civil engineering structures, as well as the determination of the shaking induced on industrial structures. The approach is based on a FFT/VTCR/IFFT methodology coupled with a multi-frequency resolution based on use of Padé approximants.

The previous FFT/VTCR/IFFT methodology, based on a frequency-by-frequency resolution, has already been shown to be an accurate way for solving transient dynamic impact problems over a low and medium frequency range. However, it could exhibit a loss of efficiency in the case of large spectrum loadings. The multi-frequency resolution developed in this paper, based on a Padé expansion of the VTCR problem in the frequency domain, can be seen as a large band analysis methodology. The main conclusions about this approach are:

• the Padé expansion gives accurate results on plate assemblies, but requires either the modal content of the structure to be a priori known, or to proceed by dichotomy with an error criterion as proposed in this article;

• to overcome this issue, a multi-point approach based on a bandwidth adaptation related to a convergence criterion can be used. This methodology provides good results with insignificant errors on FRFs and on time displacement and acceleration responses. However, the results depend on the eigenmodes excited by the input load, and more particularly by specific modes such as floor modes or wall modes which may vary between the studied area and the area used to define the frequency bandwidth;

• application of the Padé expansion is efficient and accurate on plate assemblies compared to the previous frequency-by-frequency approach. In the case studied in this article, the time saving was estimated to be nearly 30 % (32 min) compared to a classic frequency-by-frequency resolution;

• however, application of the proposed strategy to structures involving shell substructures seems to be more problematic. Indeed the method requires the calculation of successive derivatives of the VTCR terms with respect to the frequency what is numerically a tough task in the implementation for shell substructures VTCR shape functions. Further developments should be carried out in order to adapt the Padé approximants approach to structures including shell parts.

2. 1

 1 Time-to-frequency-to-time domain strategy

Figure 1 :

 1 Figure 1: The time-frequency strategy.

FrequencyFigure 2 :

 2 Figure 2: Geometry of a plate Ω i .

Figure 3 :

 3 Figure 3: The reference problem.

Figure 4 :

 4 Figure 4: The discretized amplitudes.

Figure 5 :

 5 Figure 5: Padé approximants example: description of the geometry and boundary conditions.

Finite

  elements/wavelength) 4 × 5 edge waves out-of-plan displacement (m)

Figure 6 :

 6 Figure 6: FRF for several sets of parameter L and M over a frequency range of 40 Hz.

Figure 7 :

 7 Figure 7: Comparison of computation time for both types of analysis.

Figures 8 and 9

 9 Figures 8 and 9 define the FRF calculated over the full frequency range from 0.1 to 350 Hz. Here Padé approximants orders of L = 4 and M = 5 were chosen. This parameter set presents the best error/computation time ratio on a small frequency range. A study on the discretization analysis is presented here. This discretization can be observed through the center frequencies chosen to represent the FRF solution. First, in Figure8, the choice of bandwidth is fixed at 20 Hz and therefore defined a priori.

Figure 8 : 20 Figure 9 Figure 9 :

 82099 Figure 8: FRF over a wide frequency range with a bandwidth of 20

Figure 10 :

 10 Figure 10: Geometry of civil engineering structure studied.

Figure 11 :

 11 Figure 11: Time loading function.

Figure 12

 12 Figure 12 describes the frequency content of the loading function applied.

Figure 12 :

 12 Figure 12: Frequency content of F P 0 versine loading function.

6000Figure 13

 13 Figure13presents the results for the point P 1. Bandwidth is refined in this case studied with the 10 % criterion based on the error between the reference VTCR frequency-by-frequency solution and the Padé approximants solution. In the case where the reference solution is not known, bandwidth refinement by dichotomy, as stated in paragraph 4.1, could be used. Here this bandwidth begins at 10 Hz and is initially refined to 4 Hz. In this case, only 13 matrix inversions are enough to present good results on a frequency bandwidth of 100 Hz.

Figure 15 :Figure 16 :

 1516 Figure 15: Time displacement response over the previous adapted wide frequency band for the point P 2.

Table 2 :

 2 Material and structural parameters of the studied steel plate.

	Notation Physical quantity Value and unit
	E	Young's modulus	210 GPa
	ν	Poisson's ratio	0.3
	ρ	mass density	7800 kg/m 3
	ξ	damping factor	0.01

Table 3 :

 3 FEM Cast3m solution (left), VTCR solution with 40 DoF (right).

Table 4 :

 4 VTCR solutions of civil engineering example at four different frequencies.

Acknowledgements

This work has been funded by AREVA. The authors would like to thank N. Moussallam for hosting this research and his valuable comments.


G Allaire

C Dapogny

A Faure

G Michailidis

SHAPE OPTIMIZATION OF A LAYER BY LAYER MECHANICAL CONSTRAINT FOR ADDITIVE MANUFACTURING

The purpose of this article is to introduce a new functional of the domain, to be used in shape optimization problems as a means to enforce the constructibility of shapes by additive manufacturing processes. This functional aggregates the self-weights of all the intermediate structures of the shape appearing in the course of its layer by layer assembly. Its mathematical analysis is performed and an algorithm is proposed to accelerate the significant computational effort entailed by the implementation of these ideas. Eventually, a numerical validation and a concrete example are discussed.

Résumé. Nous introduisons dans cet article une nouvelle fonctionnelle dépendant du domaine qui, utilisée comme contrainte dans un problème d'optimisation de forme, impose la constructibilité par les procédés de fabrication additive. Cette fonctionnelle agrège les poids propres de toutes les structures intermédiaires de la forme mises en jeu au cours du processus d'assemblage par strates. Après son analyse mathématique, nous proposons un algorithme pour accélérer significativement les calculs coûteux entraînés par l'implémentation de ces idées. Une validation numérique ainsi qu'un exemple concret sont enfin présentés.

Introduction

The additive manufacturing technologies have demonstrated a unique potential in realizing structures with a high degree of complexity, thereby allowing to process almost directly the designs predicted by shape and topology optimization algorithms [START_REF] Gibson | Additive manufacturing technology: rapid prototyping to direct digital manufacturing[END_REF]. These breakthroughs come along with new opportunities, and with new challenges; see [START_REF] Clausen | Topology optimization for additive manufacturing[END_REF] for a detailed and comprehensive overview of both, and the references therein.

To summarize its main features in a few words, additive manufacturing is a common label for quite different methodologies, which share the fact that the construction process starts with a slicing procedure: the Computer Aided Design (CAD) model for the input shape (which is often supplied by means of a mesh under the popular STL format [START_REF]Stereolithography Interface Specification[END_REF]) is converted into a series of two-dimensional layers (see Figure 1). Thence, these layers are assembled individually, one above the other, according to the selected technology. As far as these technologies are concerned, two important categories are the following:

• Material extrusion methods, such as Fused Filament Fabrication (FFF), act by selectively extruding the molten material through a nozzle; such methods are typically used to process plastic (ABS). • Powder bed fusion methods (such as Selective Laser Sintering, or Electron Beam Melting) are on the contrary used to process metals; at the beginning of the construction of each layer, metallic powder is spread within the build chamber and a laser (or an electron beam) is used to bind the grains together. These technologies show competing features in terms of cost, speed, accuracy, ... and also in terms of the restrictions they impose on the manufacturing process and the constructed shape. Beyond their differences, one challenge is faced by all additive manufacturing technologies, that of building shapes showing large overhangs, i.e. regions hanging over void (or powder) without sufficient support from the lower structure.

• In the case of material extrusion methods, parts of the boundary showing large overhangs cannot be produced as is, since this demands depositing material on void. • In the case of powder bed fusion methods, the rapid melting then solidification of the material induces large thermal variations in the structure; this creates residual stresses, and causes the structure to warp. This phenomenon is all the more likely to occur in regions which are unanchored to the lower structure (in particular, overhanging regions); see [START_REF] Mumtaz | A method to eliminate anchors/supports from directly laser melted metal powder bed processes[END_REF]. Another source of difficulties in the assembly of overhanging regions lies in that the fused material may drip between the unfused powder of the lower structure, thus leaving the processed boundary with rough patches [START_REF] Calignano | Design optimization of supports for overhanging structures in aluminum and titanium alloys by selective laser melting[END_REF]. One way to cope with the presence of overhangs is to erect scaffolds (or supports) at the same time that the shape is constructed, with the purpose of anchoring the overhanging regions [START_REF] Dumas | Bridging the Gap: Automated Steady Scaffoldings for 3D Printing[END_REF]. This scaffold structure has to be removed manually at the end of the process, which is cumbersome and time-consuming. Another way is to constrain the presence of overhangs in the formulation of the shape optimization problem guiding the design of the shape. Hitherto, ad hoc criteria, based on a minimum angle between the structural boundary and the horizontal directions, have been used to tackle this issue [START_REF] Gaynor | Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design[END_REF][START_REF] Langelaar | Topology optimization of 3D self-supporting structures for additive manufacturing[END_REF][START_REF] Mirzendehdel | Support structure constrained topology optimization for additive manufacturing[END_REF].

The present article is devoted to the modeling and the mathematical analysis of a new mechanical constraint for the optimization of shapes which are processed by an additive manufacturing method. Several variants of this contraint, and various numerical examples and discussions of engineering applications will be presented in a forthcoming article [START_REF] Allaire | Structural optimization under overhang constraints imposed by additive manufacturing technologies[END_REF]. Under the simplifying assumption that the components of one single layer of material are built simultaneously during the manufacturing process, we introduce a new constraint functional for shape optimization problems, which appraises the constructibility of shapes at each stage of their assembly. In particular, overhang constraints are naturally addressed by this formulation which appeals to their mechanical origin. To achieve our purpose, in the setting of the optimization problem, we distinguish the mechanical situation where the final shape Ω is utilized, on which the optimization criterion is based, and that where Ω (and all the successive, intermediate structures) is under construction, which guides the definition of our constraint functional. Our first main result is to provide a shape derivative for this new constraint functional (see Theorem 3.1), which is not a completely standard matter since the upper boundaries of the intermediate shapes are not subject to optimization, being dictated by the additive manufacturing process.

Our second main result is an acceleration method for the computation of our new constraint functional and for that of its shape derivative. Indeed, their expressions involve mechanical problems posed on all the intermediate structures of the considered shape, the number of which is precisely the number of layers in the additive manufacturing process (typically of the order of a few hundreds). Therefore, these evaluations are quite costly in numerical practice. Our idea is to interpolate, with piecewise affine functions of the height, the values of the functional and of its derivative, thus relying on the derivatives of the mechanical solutions with respect to the height of the intermediate structures.

This article is organized as follows. In Section 2, we introduce our shape optimization problem. In Section 3, we describe the mechanical context in which shapes are constructed, we formulate our new manufacturing constraint and we prove our first main mathematical result, Theorem 3.1, concerning its shape derivative. As we have mentionned, the resulting functional of the domain and its shape derivative are costly to evaluate in practice. Thus, in Section 4 we propose an interpolation method for accelerating significantly these calculations. For this purpose we introduce a variant of the Hadamard method of shape deformations where only the upper horizontal boundary of the intermediate structures is allowed to vary, while the rest of the boundary is fixed. In this setting, first-order Taylor approximations of the mechanical performances of the intermediate shapes can be computed in terms of the height. Eventually, a numerical validation of our acceleration process and an optimization example are provided in Section 5.

Presentation of the shape optimization problem

A shape is a bounded, regular domain Ω ⊂ R d , d = 2, 3, filled with a linear elastic material with Hooke's law A. In the context of its final utilization, Ω is clamped on a subset Γ D ⊂ ∂Ω, and it is submitted to surface loads f ∈ L 2 (Γ N ) d applied on a region Γ N of ∂Ω disjoint from Γ D ; the remaining part of the boundary Γ :

= ∂Ω \ (Γ D ∪ Γ N ) is traction-free. The elastic displacement u m Ω is the unique solution in H 1 Γ D (Ω) d := u ∈ H 1 (Ω) d
, u = 0 on Γ D to the mechanical system:

(2.1)

       -div(Ae(u m Ω )) = 0 in Ω, u m Ω = 0 on Γ D , Ae(u m Ω )n = 0 on Γ, Ae(u m Ω )n = f on Γ N .
For simplicity, the objective J(Ω) driving the optimization problem is the compliance:

(2.2) J(Ω) = Ω Ae(u m Ω ) : e(u m Ω ) dx = Γ N f • u m Ω ds.
Our optimization problem then reads:

(2.3) min

Ω∈U ad J(Ω), such that P (Ω) ≤ α.
In (2.3), U ad is a set of admissible shapes Ω, which are assumed to be of class C ∞ for simplicity, and whose boundaries enclose the non optimizable regions Γ D , Γ N and Γ 0 (the latter is defined in Section 3 below), i.e.

U ad = Ω ⊂ R d is open, bounded, and of class C ∞ , Γ D ∪ Γ N ∪ Γ 0 ⊂ ∂Ω , P (
Ω) is our new constraint functional, whose definition and properties are discussed in the next sections, and α is a tolerance threshold. In practice, so that (2.3) be physically relevant, other constraints (e.g. on the volume Vol(Ω) of shapes) may be added, which we omit for the mathematical analysis; see Section 5.

Most popular optimization algorithms for the numerical resolution of (2.3) rely on the derivatives of J(Ω) and P (Ω) with respect to the domain; these are understood in the framework of Hadamard's method (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Soko Lowski | Introduction to shape optimization: shape sensitivity analysis[END_REF]): variations of a shape Ω are considered under the form:

Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1.
A generic function F (Ω) of the domain is then shape differentiable if the underlying mapping θ → F (Ω θ ), from W 1,∞ (R d , R d ) into R, is Fréchet differentiable at θ = 0; the corresponding derivative is denoted by F (Ω)(θ). In practice, the deformations θ featured in this definition are restrained to a subset of W 1,∞ (R d , R d ); in the following, we shall consider the sets

Θ k = θ ∈ C k,∞ (R d , R d ), θ = 0 on Γ D ∪ Γ N ∪ Γ 0 ,
where k ≥ 1, and C k,∞ (R d , R d ) is the set of k times continuously differentiable functions from R d into itself, whose derivatives up to order k are uniformly bounded.

For instance, if f is smooth, it is well-known (see [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]) that the objective (2.2) is shape differentiable when deformations are restrained to Θ k , k ≥ 1, and that its shape derivative reads:

∀θ ∈ Θ k , J (Ω)(θ) = - Γ Ae(u m Ω ) : e(u m Ω ) θ • n ds.

Description and analysis of the mechanical constraint

In this section, we introduce and analyze mathematically our new mechanical constraint functional P (Ω) describing the manufacturing process of shapes.

Formulation of the constraint functional P (Ω).

The constraint P (Ω) relies on the mechanical situation of Ω in the course of the manufacturing process: assuming a vertical build direction e d (the d th vector of the canonical basis (e 1 , ..., e d ) of R d ), Ω is enclosed in a box D = S × (0, H), where S ⊂ R d-1 . In practice, D represents the build chamber. For h ∈ (0, H),

(3.1) Ω h := Ω ∩ x = (x 1 , ..., x d ) ∈ R d , 0 < x d < h
is the intermediate shape describing the stage where the final shape Ω is assembled up to height h. The boundary ∂Ω h is decomposed in a different fashion from that of Section 2:

∂Ω h = Γ 0 ∪Γ u h ∪Γ l h , where
• Γ 0 = {x ∈ ∂Ω h , x d = 0}
is the contact region between Ω and the build table,

• Γ u h = {x ∈ ∂Ω h , x d = h} is the upper side of the intermediate structure,

• Γ l h = ∂Ω h \ (Γ 0 ∪ Γ u h
) is the lateral surface. Eventually, we define h := {x ∈ ∂Ω, x d = h}, the part of the boundary ∂Ω lying at height h (typically a curve in three space dimensions); see Figure 2 about these notations.

Each intermediate shape Ω h is clamped on Γ 0 , and is only subjected to gravity effects, accounted for by a body force g

∈ L 2 (R d ) d . Its elastic displacement u c Ω h ∈ H 1 Γ0 (Ω h ) d satisfies the system: (3.2)    -div(Ae(u c Ω h )) = g in Ω h , u c Ω h = 0 on Γ 0 , Ae(u c Ω h )n = 0 on Γ l h ∪ Γ u h , The compliance c Ω h of Ω h then reads: (3.3) c Ω h = Ω h Ae(u c Ω h ) : e(u c Ω h ) dx = Ω h g • u c Ω h dx.
Our constraint P (Ω) of the final structure Ω aggregates the compliances of all the intermediate shapes:

(3.4) P (Ω) = H 0 j(c Ω h ) dh,
where j : R → R is a given, smooth function. Note that, as pointed out in the introduction, P (Ω) only involves the intermediate stages Ω h of the construction of Ω where the successive layers are completed (and not all the stages where these layers are themselves under construction, and partially assembled). 3.2. Differentiability of P (Ω) with respect to the domain.

Throughout this section, we consider a fixed shape Ω ∈ U ad . The rigorous exposition of the shape differentiability analysis of P (Ω) requires that we introduce two open sets O 1 O 2 in R d and a smooth function χ : R d → R such that:

x ∈ ∂Ω \ Γ 0 , n(x) • e d = ±1 ⊂ O 1 , 0 ≤ χ ≤ 1, χ ≡ 0 on O 1 , and χ ≡ 1 on R d \ O 2 .
In other words, O 1 is an open neighborhood of the 'flat horizontal regions' of ∂Ω \ Γ 0 , and χ is a cutoff function whereby these regions will be ignored. Using these notations, the relevant sets for deformations of Ω are the Banach spaces

(3.5) X k = θ = χ θ, θ ∈ Θ k , equipped with the quotient norm ||θ|| X k = inf || θ|| C k,∞ (R d ,R d ) , θ = χ θ .
Among other things, vector fields θ ∈ X k vanish near the points of ∂Ω \ Γ 0 where the normal vector n is parallel to e d .

Let us explain the roles of the cutoff function χ and of the space X k of shape perturbations:

• Deformations Ω θ of Ω, with θ ∈ X k , can be equivalently described by 'horizontal' perturbations (see Proposition 3.2), which are the only ones for which the shape derivative of the compliances c Ω h can be rigorously calculated owing to Lemma 3.1.

• Deformations θ ∈ X k vanish around any point x ∈ ∂Ω where n(x) is 'vertical'. This is because at such points x (where the tangent plane is 'horizontal') the intermediate structure Ω h at height h = x d shows a turning point at x, and therefore may be not even Lipschitz regular around x. The poor regularity of u c Ω h in this region (see [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF] §3.3 for a related study) would be an obstruction to our mathematical analysis. The main result of this section is the following: Theorem 3.1. The functional P (Ω) given by (3.4) is shape differentiable at Ω, in the sense that the mapping

θ → P (Ω θ ), from X k into R is differentiable for k ≥ 1.
Its derivative is:

(3.6) ∀θ ∈ X k , P (Ω)(θ) = ∂Ω\Γ0 D Ω θ • n ds,
where the integrand factor D Ω is defined, for a.e. x ∈ ∂Ω \ Γ 0 , by:

(3.7) D Ω (x) = H x d j (c Ω h ) 2g • u c Ω h -Ae(u c Ω h ) : e(u c Ω h ) (x) dh.
The shape sensitivity P (Ω) expressed in Theorem 3.1 does not result so easily from standard arguments, since P (Ω) involves all the intermediate structures Ω h of Ω, which are only Lipschitz regular (in particular, they show angles at the tip of the upper boundary). We perform the proof in several steps.

(1) In Section 3.2.1, we start by proving that θ → P (Ω θ ) is differentiable if θ is restricted to the subset

X k H := θ ∈ X k , θ • e d = 0 of 'horizontal' perturbations in X k ,
and we show that Formulae (3.6), (3.7) hold in this case.

(2) We prove in Section 3.2.2 that for a given θ ∈ X k , there exists a horizontal deformation ξ ≡ ξ(θ) ∈ X k H accounting for the same perturbed shape: Ω θ = Ω ξ(θ) . We also prove that the mapping θ → ξ(θ) is differentiable, and we calculate its derivative.

(3) Theorem 3.1 arises in Section 3.2.3 as a consequence of chain rule and of the two previous points. Remark 3.2. Formulae (3.6) and (3.7) have an intuitive structure: the shape gradient of P (Ω) at a point x = (x 1 , ..., x d ) ∈ ∂Ω \ Γ 0 involves the elastic energy in x for all the intermediate structures Ω h , h > x d .

3.2.1.

Step 1: Shape differentiability of θ → P (Ω θ ) when θ ∈ X k H .

Let us start with the following lemma:

Lemma 3.1. The compliance θ → c (Ω θ ) h at level h, defined by (3.3), is Fréchet differentiable over X k H for k ≥ 1,
and the corresponding derivative reads:

(3.8) ∀θ ∈ X k H , c Ω h (θ) = Γ l h 2g • u c Ω h -Ae(u c Ω h ) : e(u c Ω h ) θ • n ds.
Proof. The key feature of horizontal deformations lies in the following relation, which holds for θ ∈ X k H small enough:

(

Ω θ ) h = {x ∈ Ω θ , 0 < x d < h} ≡ (Ω h ) θ = (Id + θ)(Ω h ).
Hence, because it only involves deformations θ ∈ X k H , Lemma 3.1 merely boils down to the differentiation of the compliance Ω → c Ω defined by (3.3) at Ω = Ω h . The only difference with the usual setting, as in e.g. [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], is that the domain Ω h is not smooth. More precisely, Ω h is piecewise smooth and exhibits corners in two dimensions, ridge edges in three dimensions, at points x ∈ h (see again Figure 2).

However, Ω h is necessarily locally convex around the sharp features formed by the points x ∈ h where the normal vector n(x) is not parallel to e d (in two dimensions, the angles corresponding to these corners are in (0, π)). As a consequence, the theory of elliptic equations in polygonal domains implies that the solution [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF], Remark 3.2.4.6). Since deformations θ ∈ X k (in particular θ ∈ X k H ) identically vanish on O 1 , which contains the 'bad points' where n(x) is parallel to e d (see (3.5)), the classical arguments (see [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], §5.3) leading to the expression (3.8) for the shape derivative of Ω → c Ω at Ω h , involving the calculation of the Lagrangian and Eulerian derivatives of the mapping Ω → u c Ω at Ω h , and elementary (but tedious) calculations based on the Green formula, can be worked out in our particular situation. Similar calculations are performed in a slightly more general context in Proposition 4.2 below, and we do not replicate the argument.

u c Ω h to (3.2) enjoys H 2 regularity in Ω h \ O 1 (see
We are now in good shape for proving Theorem 3.1 in the special case where θ ∈ X k H . Proposition 3.1. The mapping θ → P (Ω θ ) defined by (3.4), from X k H into R (for k ≥ 1), is Fréchet differentiable at θ = 0. Its shape derivative reads:

∀θ ∈ X k H , P (Ω)(θ) = ∂Ω\Γ0 D Ω θ • n ds
, where D Ω is given by (3.7).

Proof. Let us first discuss the shape differentiability of P (Ω); denoting by m(θ, h) = j(c (Ω h ) θ ), it follows from Lemma 3.1 (and arguments similar to those involved in its proof) that:

• The mapping (θ, h) → m(θ, h) is continuous on X × (0, H), where X is a neighborhood of 0 in X k H . • For any h ∈ (0, H), the mapping θ → m(θ, h) is Fréchet differentiable on X . • The (partial) Fréchet derivative (θ, h) → ∂m ∂θ (θ, h) of θ → m(θ, h) is continuous from X × (0, H) into the dual space of X k
H . Then, it follows from the Lebesgue dominated convergence theorem and Lemma 3.1 that θ → P (Ω θ ) is Fréchet differentiable at θ = 0, and that its derivative reads:

∀θ ∈ X k H , P (Ω)(θ) = H 0 j (c Ω h ) c Ω h (θ) dh.
Recalling the expression (3.8) of c Ω h (θ) and using the shorthand

I(x, h) ≡ (2g • u c Ω h -Ae(u c Ω h ) : e(u c Ω h ))(x), this rewrites: (3.9) P (Ω)(θ) = H 0 j (c Ω h ) Γ l h I(x, h)(θ • n)(x) ds(x) dh.
The rest of the proof relies on repeated applications of the Fubini theorem, which in particular entails:

Γ l h ρ(x) dx = h 0 z ρ(x) d (x) dz,
for a (smooth) function ρ. From (3.9), we obtain successively:

P (Ω)(θ) = H 0 h 0 z j (c Ω h )I(x, h)(θ • n)(x) d (x) dz dh = H 0 H z z j (c Ω h )I(x, h)(θ • n)(x) d (x) dh dz = H 0 z H x d j (c Ω h )I(x, h) dh (θ • n)(x) d (x) dz = ∂Ω\Γ0 H x d j (c Ω h )I(x, h) dh (θ • n)(x) ds(x),
which is the desired conclusion.

3.2.2.

Step 2: Parameterization by horizontal perturbations.

In this section, we prove that 'horizontal' perturbations θ ∈ X k H are the only needed ingredient to describe variations of Ω of the form Ω θ , at least when θ ∈ X k (here again, the cutoff function χ plays a key role). For this purpose, it is convenient to introduce an alternative way to express identities of the form Ω θ = Ω ξ , for θ, ξ ∈ X k . Let us consider an implicit representation for Ω, that is, a smooth function φ : R d → R such that:

(3.10)    φ(x) < 0 if x ∈ Ω, φ(x) = 0 if x ∈ ∂Ω, φ(x) > 0 if x ∈ R d \ Ω. For arbitrary vector fields θ, ξ ∈ W 1,∞ (R d , R d ) with ||θ|| W 1,∞ (R d ,R d ) < 1, ||ξ|| W 1,∞ (R d ,R d ) < 1, we define F(θ, ξ) = φ • (Id + θ) -1 • (Id + ξ)
as an element in the set C(∂Ω) of continuous functions on ∂Ω.

Lemma 3.2. Let θ, ξ ∈ W 1,∞ (R d , R d ) be such that (3.11) ||θ|| W 1,∞ (R d ,R d ) < 1, and ||ξ|| W 1,∞ (R d ,R d ) < 1 .
Then the domains Ω θ and Ω ξ coincide if and only if F(θ, ξ) = 0 in C(∂Ω).

Proof. We only need to discuss the 'if' part of the above statement: assume that F(θ, ξ) = 0. Because of (3.11), Ω θ and Ω ξ are bounded, Lipschitz regular domains which are Lipschitz diffeomorphic, see [START_REF] Evans | Measure theory and fine properties of functions[END_REF]. For the same reason, ∂Ω θ and ∂Ω ξ are compact, Lipschitz submanifolds of R d which are also Lipschitz diffeomorphic. Moreover, since the function φ • (Id + θ) -1 implicitly describes Ω θ (in the sense that (3.10) holds with Ω θ instead of Ω), and since

F(θ, ξ) = 0 implies that φ • (Id + θ) -1 vanishes on (Id + ξ)(∂Ω) = ∂Ω ξ , if follows that ∂Ω ξ ⊂ ∂Ω θ .
Let us now write ∂Ω ξ = N i=1 B i as the disjoint reunion of its connected components; each B i is a compact, Lipschitz submanifold in R d , and so it is a connected component of ∂Ω θ too. Hence, if the inclusion ∂Ω ξ ⊂ ∂Ω θ were strict, there would exist one Lipschitz submanifold B ⊂ R d , disjoint from ∂Ω ξ , such that ∂Ω θ = B ∪ ∂Ω ξ , in contradiction with the fact that ∂Ω θ and ∂Ω ξ are Lipschitz diffeomorphic.

Before proceeding, let us introduce an additional notation: when x ∈ ∂Ω, n H (x) := n(x) -(n(x) • e d )e d denotes the orthogonal projection of the normal vector to Ω on the 'horizontal space', spanned by e 1 , ..., e d-1 .

Recall (see [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], Prop. 5.4.14) that the normal vector field n can be extended from ∂Ω to R d as a whole, into a vector field of class C ∞ which has unit norm in a neighborhood of ∂Ω. In the following, with a small abuse of notations, we still denote by n (resp. n H ) this extended normal vector field (resp. its projection on the horizontal space).

The desired result is the following; see Figure 3 for an illustration.

Proposition 3.2. For every k ≥ 1, there exists a mapping θ → ξ(θ), from a neighborhood X of 0 in X k+1 into X k H such that, for θ ∈ X , Ω θ and Ω ξ(θ) coincide. In addition, θ → ξ(θ) is Fréchet differentiable on X and the restriction to ∂Ω of its derivative at θ = 0 reads:

(3.12) ∀θ ∈ X k+1 , ξ (0)(θ) ∂Ω = 1 |n H | 2 (θ • n) n H .
Proof. The proof is decomposed into three steps.

(i): Let us define the Banach spaces

F k = ζ ∈ C k (∂Ω), ζ = 0 on Γ 0 ∪ Γ D ∪ Γ N , and
(3.13) Y k = ζ = χ ζ, ζ ∈ F k , equipped with the norm ||ζ|| Y k = inf || ζ|| C k (∂Ω) , ζ ∈ F k s.t. ζ = χ ζ .
Introducing a sufficiently small neighborhood X (resp. Y) of 0 in X k+1 (resp. in Y k ), let us define:

(3.14) ∀θ ∈ X , ∀ζ ∈ Y, G(θ, ζ) = φ • (Id + θ) -1 • (Id + ζn H ) ∈ C k (∂Ω). Our first observation is that G actually maps X × Y into Y k . Indeed, if θ ∈ X and ζ ∈ Y, there exist θ ∈ Θ k and ζ ∈ F k such that θ = χ θ and ζ = χ ζ.
Then, for x ∈ ∂Ω, letting z = (Id + ζn H )(x), we first calculate:

(3.15) (Id + θ) -1 (z) -z = (Id + θ) -1 (z) -(Id + θ) -1 (z + θ(z)), = χ(z) 1 0 ∇((Id + θ) -1 )(z + tθ(z)) θ(z) dt.
Using Taylor's formula on χ(z) in the right-hand side of the above formula yields:

χ(z) = χ(x) + χ(x) ζ(x) 1 0 ∇χ(x + tζ(x)n H (x)) • n H (x) dt.
It follows that there exists r ∈ Θ k such that:

(Id + θ) -1 • (Id + ζn H ) = Id + χr.
Then, for arbitrary x ∈ ∂Ω, we obtain:

G(θ, ζ)(x) = (φ • (Id + θ) -1 (Id + ζn H ))(x) = φ(x + χ(x)r(x)) -φ(x), = χ(x) 1 0 ∇φ(x + tχ(x)r(x)) • r(x) dt, which confirms that G(θ, ζ) ∈ Y k .
(ii): Using analogous arguments, it is easily seen that G is differentiable; its partial differential at (0, 0) with respect to the ζ variable reads:

(3.16) ∀ ζ ∈ Y k , d ζ G(0, 0)( ζ) = (∇φ • n H ) ζ = |n H | 2 ζ.
Since the partial differential d ζ G(0, 0) defined by (3.16) is invertible from Y k into itself, the Implicit Function theorem (see e.g. [START_REF] Lang | Fundamentals of differential geometry[END_REF], Chapter I, Theorem 5.9) allows to conclude that, after possibly taking smaller neighborhoods X and Y of 0 in X k+1 and Y k respectively, there exists a mapping

X θ → ζ(θ) ∈ Y of class C 1 such that: ∀θ ∈ X , ζ ∈ Y, G(θ, ζ) = 0 ⇔ ζ = ζ(θ).
(iii): At this point, the only remaining operation is to extend the scalar function ζ(θ) : ∂Ω → R into a vector field defined on R d as a whole. To this end, we denote by p ∂Ω (x) the closest point on ∂Ω to an arbitrary point x ∈ R d ; since ∂Ω is smooth, x → p ∂Ω (x) is well-defined and smooth on a tubular neighborhood V of ∂Ω; see again [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], Prop. 5.4.14. Let also γ : R d → R d be a smooth function such that γ ≡ 1 on a smaller neighborhood of ∂Ω and γ ≡ 0 on R d \ V. We now define, for θ ∈ X ,

ξ(θ) = (ζ(θ) • p ∂Ω ) γn H ∈ X k .
It readily follows from Lemma 3.2 that Ω θ = Ω ξ(θ) . Eventually, differentiating the relation G(θ, ζ(θ)) = 0 with respect to θ, then evaluating at θ = 0, we obtain the Fréchet derivative of θ → ζ(θ) at θ = 0:

∀θ ∈ X k+1 , ζ (0)(θ) = 1 |n H | 2 θ • n,
which readily leads to (3.12) and terminates the proof. Remark 3.3. Notice that the horizontal deformation ξ(θ) supplied by Proposition 3.2, giving rise to the same variation Ω θ of Ω than the argument θ, has one degree of regularity less than θ. This technical point is a side effect of our application of the Implicit Function theorem to the function G defined by (3.14) and of the underlying choice (3.13) of functional spaces (see in particular Formula (3.15) where we need one more derivative for θ). We do not know whether this result can be improved.

End of the proof of Theorem 3.1.

For k ≥ 1, let θ → ξ(θ) be the mapping from X k+1 into X k H supplied by Proposition 3.2. Then, as a consequence of definitions, it holds, for θ ∈ X k+1 small enough that Ω θ = Ω ξ(θ) and thus

P (Ω θ ) = P (Ω ξ(θ) ).
The combination of the chain rule with Proposition 3.1 allows to conclude. Remark 3.4. Formulae (3.6)-(3.7) for the shape derivative of J(Ω) can be retrieved in a formal way by using the interesting results in [START_REF] Santambrogio | A semi derivation lemma on BV functions[END_REF][START_REF] Vialard | Hamiltonian approach to geodesic image matching[END_REF], about the differentiation of functions of the form:

t → D (a • Φ t ) b dx,
where a and b are special functions with bounded variations (e.g. characteristic functions), and t → Φ t is the flow generated by a vector field. Remark 3.5. A careful investigation of the proof of Theorem 3.1 shows that it was not necessary to assume that the final shape Ω ∈ U ad is of class C ∞ . Rather a shape of class C k+1 is enough.

⌦ D 0 ⌦ ✓ ✓ ⇠(✓)

Practical calculation of the mechanical constraint and its derivative

The numerical evaluation of P (Ω) and P (Ω)(θ), or equivalently D Ω , by means of Formulae (3.4) and (3.6), (3.7) relies on a discretization of the height interval (0, H) with a sequence 0 = h 0 < h 1 < ... < h N = H. The intuitive, '0 th -order' method to calculate approximations P 0 N and D 0 N of P (Ω) and D Ω consists in replacing c Ω h and u c Ω h by piecewise constant quantities on each interval I i := (h i , h i+1 ) before applying (3.4) and (3.7):

(4.1) c Ω h ≈ c Ω h i+1 and u c Ω h ≈ u c Ω h i+1 on Ω h , for h ∈ I i = (h i , h i+1 ).
This procedure is costly since the piecewise constant approximation (4.1) is low-order: so that its accuracy is guaranteed, the subdivision {h i } i=0,...,N of (0, H) has to be quite fine, which brings about many numerical resolutions of the elasticity system (3.2) for the u c Ω h i . The efficiency of the '0 th -order' method can be improved by constructing a higher-order reconstruction of the mappings h → c Ω h and h → u c Ω h on each interval I i . This requires the calculation of the derivatives of these mappings in an adequate sense, which is the main purpose of this section.

A review of shape differentiation using diffeomorphisms.

This section takes place in a slightly different setting from that of Section 3 where shape derivatives were computed by Hadamard's method, as described in Section 2, i.e., by differentiating functionals of the type θ → F (Ω θ ). Here, we rather rely on shape variations described by a parameter-dependent diffeomorphism t → T t . This change in point of views will come in handy in Section 4.2 below: it will allow to describe intermediate shapes Ω h-t close to Ω h (see (3.1) for the definition of Ω h ) as variations of Ω h , i.e. Ω h-t = T t (Ω h ). This paves the way to a 'natural' notion of differentiation of quantities such as the compliance h → c Ω h defined by (3.3)

and the elastic behavior h → u c

Ω h given by (3.2) of the intermediate shapes Ω h with respect to the height parameter.

For the moment, we drop the index h and we consider a bounded domain Ω in R d , which is only assumed to be Lipschitz regular (as is Ω h ). Its boundary reads as the disjoint reunion ∂Ω = Γ 0 ∪ Γ, where Γ 0 is a non optimizable subset of ∂Ω of positive (d -1)-dimensional Hausdorff measure.

In this context, we denote by v Ω ∈ H 1 Γ0 (Ω) d the unique solution to the system:

   -div(Ae(v Ω )) = g in Ω, Ae(v Ω )n = 0 on Γ, v Ω = 0 on Γ 0 ,
where g is a given function in H 1 (R d ) d . As announced above, variations of Ω are performed by means of a mapping t → T t , defined on the interval (-t 0 , t 0 ) for some t 0 > 0, which satisfies the properties:

(4.2) For any t ∈ (-t 0 , t 0 ), T t is a diffeomorphism of R d such that T t (Γ 0 ) = Γ 0 . (4.3) For any t ∈ (-t 0 , t 0 ), (T t -Id) ∈ W 1,∞ (R d , R d ),
and

the mapping t → (T t -Id) is of class C 1 from (-t 0 , t 0 ) into W 1,∞ (R d , R d ).
We define

V (x) := dTt(x) dt | t=0 ∈ W 1,∞ (R d , R d ).
We now calculate the Eulerian and Lagrangian derivatives of the mapping Ω → v Ω , with respect to variations of Ω driven by T t . Although the involved arguments are quite classical (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Novotny | Topological derivatives in shape optimization[END_REF]), the (tedious) proof of these formulae is not so easily found in the literature in the context of the linearized elasticity system. For the sake of convenience, we recall the main steps in the Appendix. Proposition 4.1. Let t → T t be a mapping satisfying (4.2) and (4.3). Then, (i) The mapping Ω → v Ω has a material derivative in the sense that the transported function t

→ v t := v Tt(Ω) • T t , from (-t 0 , t 0 ) into H 1 Γ0 (Ω) d
is Fréchet differentiable at t = 0; its derivative vΩ satisfies:

(4.4)    -div(Ae( vΩ )) = div(g ⊗ V + (divV )Ae(v Ω ) -AC(v Ω , V ) -Ae(v Ω )∇V T ) in Ω, Ae( vΩ )n = -(divV )Ae(v Ω )n + AC(v Ω , V )n + Ae(v Ω )∇V T n on Γ vΩ = 0 on Γ 0 ,
where we have used the shorthand:

(4.5) C(v, V ) = 1 2 (∇v∇V + ∇V T ∇v T ). (ii) Assuming that ∇v Ω V ∈ H 1 Γ0 (Ω) d , the mapping Ω → v Ω has a Eulerian derivative v Ω := vΩ -∇v Ω V in H 1 Γ0 (Ω) d
, where v Ω is the solution to the system:

(4.6)    -div(Ae(v Ω )) = 0 in Ω, Ae(v Ω )n = -∂ ∂n ((Ae(v Ω )n) (V • n) + Ae(v Ω )(∇ Γ (V • n)) on Γ, v Ω = 0 on Γ 0 , with ∇ Γ ζ = ∇ζ -(∇ζ • n)n, the tangential gradient of a (smooth enough) function ζ : Γ → R.
Remark 4.1. A word about notations. For (smooth) vector fields v, w : R d → R d , we denote by ∇v the d × d Jacobian matrix of v, that is, the matrix with entries: (∇v) ij = ∂vi ∂xj , i, j = 1, ..., d. Accordingly, ∇vw is the vector field with components (∇vw

) i = (w • ∇)v i := d j=1 ∂vi ∂xj w j .
As an easy consequence, we obtain the following result about the shape differentiation of the compliance

c Ω = Ω g • v Ω dx = Ω Ae(v Ω ) : e(v Ω ) dx.
Corollary 4.1. Let t → T t be a mapping satisfying (4.2) and (4.3). Then, t → c Tt(Ω) is differentiable at t = 0, and its derivative reads:

(4.7) d dt (c Tt(Ω) ) t=0 = Γ (2g • v Ω -Ae(v Ω ) : e(v Ω )) V • n ds.

Using shape variations to identify close layers.

In this section, we consider a fixed shape Ω ∈ U ad , and a height h ∈ (0, H) satisfying:

(4.8) for any x ∈ h , the normal vector n(x) is not parallel to e d .

Our purpose is to show that, for t > 0 small enough, the intermediate structures Ω h and Ω h-t can be expressed in terms of one another via a shape variation of the form (4.2), (4.3), namely that there exists a diffeomorphism T t of R d such that T t (Ω h ) = Ω h-t (see Figure 4). Then we shall use the material in Section 4.1 to differentiate the mappings h → c Ω h and h → u c Ω h ; see Section 4.3 below.

h h t • • h h t x T t (x) • V (x) • • • • • Figure 4. Example of one diffeomorphism T t of d mapping Ω h onto Ω h-t .
Lemma 4.1. Under the assumption (4.8), there t 0 > 0 and a mapping (-t 0 , t 0 ) t → T t satisfying (4.2) and (4.3), as well as the additional property:

(4.9) For t ∈ (-t 0 , t 0 ), T t is a diffeomorphism from Ω h onto Ω h-t .
For any such mapping, let

V ∈ W 1,∞ (R d , R d ) be defined by V (x) = dTt(x) dt t=0
. Then,

(i) For any point x ∈ Γ u h , V (x) • e d = -1. (ii) For any point x ∈ Γ 0 ∪ Γ l h , V (x) • n(x) = 0.
Proof. First observe that, under the hypothesis (4.8), the existence of a mapping t → T t satisfying (4.2), (4.3) and (4.9) follows from elementary considerations of differential geometry. Indeed, since (4.8) holds, there exist t 0 > 0 and a smooth diffeomorphism τ of R d which maps the slice

K := {x ∈ Ω, h -t 0 < x d < h + t 0 } onto the straight cylinder Γ u h × (h -t 0 , h + t 0 ) ⊂ R d-1 × R in such a way that: ∀x = (x 1 , ..., x d ) ∈ R d , π d (τ (x)) = x d ,
where π d : R d-1 × R → R is the standard projection onto the last coordinate: π d (y, s) = s for y ∈ R d-1 and y ∈ R; see Figure 5. Let us now introduce a smooth function l : (-t 0 , t 0 ) × R → R such that:

• For every t ∈ (-t 0 , t 0 ), l(t, •) : R → R is a strictly increasing, one-to-one function.

• For t ∈ (-t 0 , t 0 ), l(t, s) = s for all s ∈ (-∞, h -t 0 ] ∪ [h + t 0 , +∞), and l(t, h) = h -t.
We finally define:

∀(y, s) ∈ R d-1 × R, L t (y, s) = (y, l(t, s)).
Then the mapping T t = τ -1 • L t • τ has the desired properties (4.2), (4.3) and (4.9).

We now turn to the proof of (i). Let x 0 ∈ Γ u h be given, and let ε > 0 be so small that the ball B ε (x 0 ) with center x 0 and radius ε is compactly contained in Ω. Let ψ be an arbitrary function of class C ∞ with compact support in B ε (x 0 ). On the one hand, a change of variables produces:

p(t) := Ω h-t ψ(x) dx = Ω h |det∇T t | ψ • T t dx,
whence, differentiating at t = 0 and using Green's formula:

(4.10) p (0) = Ω h ((divV )ψ + ∇ψ • V ) dx = Γ u h ∩Bε(x0) ψ V • e d ds.
On the other hand, since ψ has compact support in B ε (x 0 ), one may alternatively perform the change of variables:

=

Ω h |det∇ T t | ψ • T t dx,
where T t is the diffeomorphism of R d defined by T t (x) = (x 1 , x 2 , ..., x d-1 , x d -t). Hence, (4.11) p (0) = Γ u h ∩Bε(x0) ψ V • e d ds,
where

V (x) = d Tt(x) dt t=0
= -1. Since both expressions (4.10) and (4.11) hold for arbitrary ψ ∈ C ∞ c (B ε (x 0 )), one infers in particular that V (x) • e d = -1, which is the desired result.

The proof of (ii) relies on similar arguments. For a given point x 0 ∈ Γ l h ∪ Γ 0 , take ε > 0 small enough so that the ball B ε (x 0 ) is compactly contained in the half-space x ∈ R d , x d < h . Let also ψ be an arbitrary function of class C ∞ with compact support in B ε (x 0 ). On the one hand, one has, for t > 0 small enough: (4.12) q(t) :=

Ω h-t ψ(x) dx = Ω h ψ(x) dx,
and so q (0) = 0. On the other hand, using the same change of variables as that leading to (4.10), we obtain:

(4.13) q (0) = Γ l h ∪Γ0
ψV • n ds.

Since both expressions (4.12) and (4.13) hold for arbitrary ψ, it follows that (V • n)(x 0 ) = 0, and the desired conclusion follows.

Derivatives of the mappings

h → c Ω h and h → u c Ω h .
Our first result is concerned with the derivative of the compliance h → c Ω h defined by (3.3). The key observation is that, for t > 0 small enough, c Ω h-t = c Tt(Ω h ) , for any mapping t → T t furnished by Lemma 4.1. Therefore, combining Corollary 4.1 with Lemma 4.1 straightforwardly yields: Proposition 4.2. Let Ω ∈ U ad , and h ∈ (0, H) be such that (4.8) holds; then the mapping h → c Ω h is differentiable at h and:

(4.14) d dh (c Ω h ) h = Γ u h (2g • u c Ω h -Ae(u c Ω h ) : e(u c Ω h )) ds.
Let us now turn to giving a suitable meaning to the derivative of h → u c Ω h . Roughly speaking, this derivative is defined as the Eulerian derivative (in the sense of Proposition 4.1) of t → u c Tt(Ω h ) , associated to any diffeomorphism t → T t mapping Ω h onto Ω h-t .

To make these considerations precise, let us summarize the results from Sections 4.1 and 4.2:

• There exist t 0 > 0 and a mapping (-t 0 , t 0 ) t → T t satisfying the properties:

(4.15) (i) For t ∈ (-t 0 , t 0 ), T t is a diffeomorphism of R d , mapping Ω h onto Ω h-t such that T t (Γ 0 ) = Γ 0 , (ii) The mapping (-t 0 , t 0 ) t → (T t -Id) ∈ W 1,∞ (R d , R d ) is of class C 1 and we define V (x) := dT t (x) dt t=0 ∈ W 1,∞ (R d , R d ). • The mapping t → u Ω c h-t • T t is differentiable from (-t 0 , t 0 ) into H 1 Γ0 (Ω h ) d . Its derivative Y Ω h at t = 0 may be interpreted as the Lagrangian derivative of h → u Ω c h . • The function U Ω h := Y Ω h -∇u Ω c h V is the solution in H 1 Γ0 (Ω h ) d to the system: (4.16)        -div(Ae(U Ω h )) = 0 in Ω h , U Ω h = 0 on Γ 0 , Ae(U Ω h )n = 0 on Γ l h , Ae(U Ω h )n = ∂ ∂n (Ae(u Ω c h )n on Γ u h .
It is natural to refer to U Ω h as the Eulerian derivative of the mapping h → u c Ω h . As is clear from (4.16) (and as expected), U Ω h is independent of the diffeomorphism t → T t used in its construction, as long as T t satisfies the intuitive properties (4.15). In particular, U Ω h does not depend on V since V (x) • e d = -1 for a.e. x ∈ Γ u h . Remark 4.2. Notice that, from a formal point of view, the Eulerian derivative U Ω h is the derivative of t → u c Ω h-t at t = 0, and not that of t → u c Ω h+t ; the reason for this seemingly unintuitive convention will find proper justification in Section 4.4 (see Formula (4.18)).

Practical algorithm.

The considerations of Section 4.3 suggest the following procedure for calculating first-order approximations P 1 N and D 1 N of P (Ω) and D Ω respectively. This allows for an accurate and computationally efficient calculation of these quantities, using a coarser subdivision {h i } i=1,...,N of (0, H) than in the calculation of the 0 th -order approximate values P 0 N and D 0 N , defined by (4.1). (1) For i = 0, ..., N calculate the compliances c Ω h i as (3.3) and the displacements u c Ω h i by solving (3.2).

(2) For i = 0, ..., N , calculate the derivative d dh (c Ω h ) h=hi of the compliance by using Proposition 4.2. (3) For i = 1, ..., N , calculate the Eulerian derivative U Ω h i at h i by using (4.16). (4) On each interval I i = (h i , h i+1 ), i = 0, ..., N -1, the compliance c Ω h is approximated by a cubic spline c i (h) which is uniquely determined by the data:

(4.17)

c i (h i ) = c Ω h i , c i (h i+1 ) = c Ω h i+1 , c i (h i ) = d dh (c Ω h ) hi , and
c i (h i+1 ) = d dh (c Ω h ) hi+1 .
(5) For i = 0, ..., N -1 and h ∈ I i = (h i , h i+1 ), u c Ω h is approximated by u h defined by: (4.18) and U Ω h i+1 are well-defined on Ω h ⊂ Ω hi+1 (see Remark 4.2).

u h (x) = u c Ω h i+1 (x) + (h i+1 -h) U Ω h i+1 (x

Numerical illustrations

Let us consider the 2d MBB Beam test case: the shapes Ω are contained in a rectangular domain D of size 6 × 1. Due to symmetry, only half of D is meshed by 300 × 100 Q 1 elements. In the context of their final utilization (described by the system (2.1)), the horizontal displacement of shapes is fixed on a small part of their lower-left side and both horizontal and vertical displacements are fixed on a small part of its lower-right side, and a unit vertical load f = (0, -1) is applied at the middle of their upper side. When it comes to their construction (modelled by (3.2)), shapes are built vertically from bottom to top, so that Γ 0 coincides with the lower side of D. The function j : R → R used in the definition (3.4) of P (Ω) is simply the identity: j(s) = s. The design Ω 0 on Figure 6 (top) is used for the numerical validation of our methods in Section 5.1 and it is the initial guess for the shape optimization of Section 5.2.

Validation of the approximations of Section 4.

At first, we calculate the functional P (Ω) and its shape derivative D Ω in the particular case where Ω = Ω 0 , by using a uniform subdivision of (0, H) made of 100 layers and the 0 th -order approximation scheme, i.e. we evaluate P 0 100 and D 0 100 , which serve as reference values for the comparisons in this section. We then calculate the 0 th -and 1 st -order approximations P i N and D i N , i = 0 and 1, associated to several subdivisions of (0, H) made of N intervals with equal length. We are interested in the behavior of the relative errors:

err(P, N, i) = |P i N -P 0 100 | P 0 100 and err(D, N, i) = ||D i N -D 0 100 || L 2 (∂Ω\Γ0) ||D 0 100 || L 2 (∂Ω\Γ0)
.

The results are displayed on Figure 6 (bottom): while the 1 st -order approximation method does not bring a lot of improvement when it comes to evaluating the constraint functional P (Ω), it allows for a substantial gain (i.e. a faster convergence with respect to the number N of subdivisions) in the evaluation of its derivative.

A numerical example.

Recalling that J(Ω) is the compliance, defined by (2.2), we now turn to the shape optimization problem:

(5.1) min

Ω∈U ad J(Ω) such that Vol(Ω) ≤ 0.2 Vol(D).
We first solve (5.1), starting from the initial design Ω 0 , by using an SLP-type algorithm in the spirit of that presented in [START_REF] Dunning | Introducing the sequential linear programming level-set method for topology optimization[END_REF], and the level set method on a fixed Cartesian mesh when it comes to tracking the deformation of shapes [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]. The optimized design Ω * is shown in Figure 7; in particular, several overhanging parts appear in Ω * . We now add our mechanical constraint P (Ω) to this problem, and solve:

(5.2) min

Ω∈U ad J(Ω) such that Vol(Ω) ≤ 0.2 Vol(D), P (Ω) ≤ 0.5 P (Ω * ).
The resulting optimized shapes, obtained by using 0 th -and 1 st -order approximations of P (Ω) and D Ω with different (uniform) subdivisions of (0, H) are represented on Figure 8. The computational effort is significantly different: about 237 h. are needed when the 0 th -order approximation process is used with N = 100 layers, whereas the total calculation takes 'only' 82 h. when using 1 st -order approximations and N = 25 layers. The values of the corresponding quantities of interest are collected in Table 1.

As is clear from the above computational times, the implementation of our algorithm has not been optimized at all. The reasons for such large CPU times are obvious, as are the possible remedies in a near future. Indeed, the optimization algorithm and the Finite Element analyses for the mechanical systems (2.1) and (3.2) are carried out in different softwares: we re-use a previous Scilab script [START_REF] Campbell | Modeling and simulation in Scilab/Scicos[END_REF] for the optimization algorithm, while we rely on FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF] for the mechanical analyses. Then, the communications between these two softwares is done through file exchanges, a notorious source of unefficiency. Also, we did not investigate the straightforward parallelization of the Finite Element resolutions of all the linear elasticity systems (3.2) posed on the intermediate structures Ω h associated to a common shape Ω which are independent one from another. Our point in giving these CPU data is only to emphasize the improvement in computational efficiency allowed by the use of the 1 st -order method instead of the 0 th -order one. Notice that, on the optimized designs of Figure 8, several overhangs placed at the lower part of the optimal shape Ω * without manufacturing constraint have vanished. Still, a few overhangs remain in the superior regions of the optimized shapes of Figure 8; this may be explained in two ways:

• The definition (3.4) of our mechanical constraint P (Ω) focuses on the performance (in terms of the self-weight (3.3)) of the lower intermediate shapes; indeed, high values of the self-weight c Ω h 0 of an intermediate structure Ω h0 generally cause high values of the self-weights c Ω h of some of the upper intermediate structures Ω h , h > h 0 , whereas the converse does not hold. • The constraint P (Ω) has been devised under the simplifying assumption that each layer of material is assembled at once, and does not bring into play the stages where these layers are themselves under construction. Hence, completely flat parts such as those observed in the designs of Figure 8 are not so 'bad' in terms of P (Ω) as long as they are anchored to the lower structure.

We refer to the companion article [START_REF] Allaire | Structural optimization under overhang constraints imposed by additive manufacturing technologies[END_REF] for further discussions about the practical use of the mechanical constraint P (Ω), and variations of it, built upon the same philosophy.

It is also remarkable that the value of the objective function is lower for the constrained problem, meaning that the constraint has the (surprising) effect of driving the algorithm in a lower local minimum (this may be due to the larger number of iterations in the latter case). Shape Ω J(Ω) Vol(Ω) P (Ω) Iterations Evaluations Figure 7 104.165 0.600 0.730 25 38 Figure 8 Inserting (A.4), (A.5), (A.6) and (A.7) into (A.3), we obtain successively:

(g ⊗ V + Ae(∇v Ω V ) + (divV )Ae(v Ω ) -AC(v Ω , V ) -Ae(v Ω )∇V T ) ij = d k=1 ∂ ∂x k (Ae(v Ω ) ij ) V k + ∂V k ∂x k Ae(v Ω ) ij -∂ ∂x k (Ae(v Ω ) ik )V j -Ae(v Ω ) ik ∂Vj ∂x k , = d k=1 ∂ ∂x k (Ae(v Ω ) ij V k ) -∂ ∂x k (Ae(v Ω ) ik V j ) ,
whence, as expected,

-div(Ae(v Ω )) i = d j,k=1 ∂ 2 ∂x j ∂x k (Ae(v Ω ) ij V k -Ae(v Ω ) ik V j )) = 0.
Let us now rearrange the boundary condition featured in (A.3). Using again (A.4), (A.6), (A.7), and the fact that Ae(v Ω )n = 0 on Γ, we obtain the following identity on Γ, for i = 1, ..., d:

(A.8) (Ae(v Ω )n) i = - d j,k=1 ∂ ∂x k (Ae(v Ω ) ij )V k n j + d j=1 Ae(v Ω ) ij (∇V T n) j , = d j,k=1 Ae(v Ω ) ij ∂V k ∂xj n k -∂ ∂x k (Ae(v Ω ) ij )V k n j .
Now, taking advantage of the fact that Ae(v Ω )n = 0 on Γ, we infer, by taking derivatives in the direction of a tangential vector field:

d j=1 ∇(Ae(v Ω ) ij ) • V Γ n j = d j=1 Ae(v Ω ) ij ∇n j • V Γ ,
where V Γ = V -(V • n)n is the tangential part of the vector field V . Hence (A.8) becomes:

(A.9) (Ae(v Ω )n) i = - Notice that, in passing from the first to the second line in (A.9), we have used the classical facts from tangential calculus (see [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], Chap. 5): ∇n = ∇n T , and ∇nn = ∇n T n = 0 on a neighborhood of ∂Ω.

Eventually, using once again these facts together with the boundary condition Ae(v Ω )n = 0 on Γ, we end up with:

Ae(v Ω )n = - ∂ ∂n ((Ae(v Ω )n) (V • n) + Ae(v Ω )(∇ Γ (V • n)) on Γ,
which is the announced result (4.6).

Figure 1 .

 1 Figure 1. Rough sketch of the slicing procedure, initiating any additive manufacturing process.

Figure 2 .Remark 3 . 1 .

 231 Figure 2. Intermediate shape Ω h at height h during the construction of the final structure Ω: the red zone is the lower boundary Γ 0 and the blue zone is the upper boundary Γ u h .

Figure 3 .

 3 Figure 3. Perturbation of a two-dimensional shape Ω by a vector field θ ∈ X k , and by the corresponding horizontal deformation ξ(θ) ∈ X k H .

Figure 5 .

 5 Figure 5. Illustration of one diffeomorphism τ of R d used in the construction of T t in the proof of Lemma 4.1.

0Figure 6 .

 6 Figure 6. (Top) Setting of the validation experiment and initial shape Ω 0 ; (bottom) relative errors of the 0 th -and 1 st -order approximations of P (Ω 0 ) and its derivative D Ω0 .

Figure 7 .

 7 Figure 7. Optimized design Ω * for the shape optimization problem (5.1).

Figure 8 .

 8 Figure 8. Optimized shapes for (5.2) using the respective approximations for P (Ω) and D Ω : (up) P 0 100 and D 0 100 ; (down) P 1 25 and D 1 25 .

  (v Ω ) ij )n j (V • n) + d j=1 Ae(v Ω ) ij d k=1 ∂V k ∂xj n k + ∇n j • V Γ , = -∂ ∂n ((Ae(v Ω )n) (V • n) + ∂n k ∂xj V k .

  ), a.e. x ∈ Ω h ; notice that the above relation does make sense for x ∈ Ω h regardless of the height h ∈ (h i , h i+1 ) since u c

	Ω h i+1

Table 1 .

 1 Values of the shape functionals and iteration numbers(A.6) (Ae(w)∇V T ) ij =

	(up)	98.484	0.599 0.343	127	143
	Figure 8 (down) 99.313	0.600 0.343	187	206

Acknowledgements. We thank Rafael Estevez for frutiful discussions during the course of this research. G. A. is a member of the DEFI project at INRIA Saclay Ile-de-France.

Now, relying on the identities, for any

and using test functions of the form w • T -1 t , w ∈ H 1 Γ0 (Ω) d in (A.1), we obtain the following variational formulation for v t :

At this point, a classical argument using the Implicit Function theorem (see e.g. [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF]) reveals that the mapping t → v t is differentiable, from (-t 0 , t 0 ) into H 1 Γ0 (Ω) d (up to decreasing the value of t 0 ). Its derivative vΩ at t = 0 is the solution to the following variational problem:

where C(v, V ) is defined by (4.5), as follows from a straightforward (yet tedious) calculation. The expression (A.2) can be rearranged owing to the following identities, valid for w ∈ H

It follows from (A.2) that:

Eventually, integrating by parts, we end up with the classical formulation (4.4) for the problem characterizing the material derivative vΩ ∈ H 1 Γ0 (Ω) d .

Proof of (ii): We now consider the Eulerian derivative v Ω of v Ω , defined from vΩ via the formula:

Using (4.4), v Ω ∈ H 1 Γ0 (Ω) d is characterized as the solution to the following problem: