Keywords: Fluid-Structure Interaction, Partitioned procedure, Wave propagation, Aeroelasticity, Shock mitigation

Arrays of obstacles have proved to be an efficient way of attenuating shock waves generated by large scale explosions. The present study intends to take into account fluid-structure interactions that may occur when elastic obstacles are used. A tractable coupling tool based on a partitioned procedure is exposed, validated on the supersonic flutter of a panel and applied to a configuration composed of square section cylinders. Several numerical difficulties related to staggering are emphasized and workarounds discussed. A methodical procedure involving one and two-way coupled simulations highlights the influence of material properties as well as the acceleration of waves in the fluid when an initial motion is prescribed to the obstacles. Finally, it is shown that uncoupled simulations may be relevant to investigate shock mitigation in some given cases.

Introduction

The prediction of damage caused on structures or on the human body by blast and shock waves, generated by large scale explosions either from terrorist attacks or industrial hazards, has been investigated for many years. Ngo et al. [START_REF] Ngo | Blast loading and blast effects on structures -an overview[END_REF] provided an overview of the effects of explosions on structures and underlined that blast resistant building design requires a detailed understanding of both blast phenomena and dynamic structural responses. Concerning damage on the human body, DePalma et al. [START_REF] Depalma | Blast injuries[END_REF] outlined various mechanisms of blast injuries owing to their type: primary blast injuries are due to the effects of over and underpressure, secondary blast injuries result however from projectiles.

In terms of shock mitigation, several methods (metallic foams, aqueous foams, textiles, sandwich panels...) [START_REF] Igra | Review of methods to attenuate shock/blast waves[END_REF] can be used to attenuate shock and blast waves in order to reduce damage during an explosion. Among them, the use of arrays of solid obstacles stands out by its efficiency and ease of implementation and sparks the interest of many researchers [START_REF] Berger | Experimental investigation on the shockwave load attenuation by geometrical means[END_REF][START_REF] Chaudhuri | Computational study of shock-wave interaction with solid obstacles using immersed boundary methods[END_REF][START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF]. In particular, an in-depth investigation of the role played by several parameters such as the shape and the arrangement of the solid obstacles on shock wave mitigation has been conducted by Chaudhuri et al. [START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF]. Nevertheless, in these studies, only rigid obstacles are considered, which does not allow to account for fluid-structure interaction (FSI) effects.

In fact, FSI simulations involving structures submitted to blast or shock wave impact are not rare. For instance, Gong [START_REF] Gong | Mutual interactions between shock waves and structures[END_REF][START_REF] Gong | Coupled fluid-structure solver: The case of shock wave impact on monolithic and composite material plates[END_REF] reported several research works concerning the structural response of plates. Deiterding et al. [START_REF] Deiterding | Parallel adaptive fluid-structure interaction simulation of explosions impacting on building structures[END_REF] performed simulations of explosions under a bridge and in a multistory building with massive structural failure. Cirak et al. [START_REF] Cirak | Large-scale fluid-structure interaction simulation of viscoplastic and fracturing thin-shells subjected to shocks and detonations[END_REF] carried out 3D venting and rupture simulations of a thin walled aluminum tube due to the passage of internal detonation. Even if some studies such as [START_REF] Giordano | Shock wave impacts on deforming panel, an application of fluidstructure interaction[END_REF] investigate the response of the fluid, most of the time, they concentrate on the structure alone.

The aim of the present study is to investigate shock/flexible obstacles interactions with the same emphasis on the wave propagation in both media (i.e. fluid and structure). More specifically, we focus on three key-points: i) understanding the role of material properties in wave propagation in the structures, ii) evaluating the effects of the structure motion on shock speed and intensity in the fluid and iii) estimating the degree of interaction between both media for a given configuration. To achieve this goal, we will restrict the study to elastic structures.

The paper is organized as follows. The computational methods, on which the following FSI simulations rely on, are presented in Section 2. This section also provides a validation of the developed coupling tool. The impact of a shock wave on a column of elastic structures is explored in Section 3 and Section 4. Section 3 constitutes a preparatory work for Section 4: the complexity of the pressure loading exerted by the fluid on the fluid-structure interface is first addressed. Then, one-way coupled simulations examine the influence of the loadings on the structure response as well as the need for an appropriate structure solver. Section 4 benefits from this groundwork to investigate the influence of the elastic deformations of the structures on the propagation of several shocks: the first part of the section tackles the influence of a prescribed structure motion on the fluid dynamics by means of a one-way coupled simulation, while its last part concerns two-way coupled simulations. Finally, the conclusion and future research directions are drawn in Section 5.

Nomenclature

The following list is restricted to the non-standard notations employed in the paper. Notations that are commonly used in structural dynamics and fluid mechanics are not reported. Length of the diagonal of the obstacles (square section cylinders) F s

Eigenmode frequency (Table 2)

H s
Length of the side of the obstacles (square section cylinders) h s min Minimum element size of the structure mesh I Pressure impulse (Eq. ( 13))

m p , m s |(v f shock .t)|/c s p , |(v f shock .t)|/c s s ratios n
Outward normal unit vector T s = (F s ) -1 Eigenmode period t

Tangent unit vector at the interface t 0 Time of arrival of the considered wave (Eq. ( 13)) t c

Impulse integration time (Eq. ( 13))

Coupling strategy and validation

This section presents the adopted coupling strategy. The governing equations as well as the partitioned procedure that has been selected to perform the coupling are detailed and validated thanks to a widely-used test case.

Governing Equations

FSI problems involve a fluid domain Ω f and a structural domain Ω s that interact at the interface Σ = Ω f ∩ Ω s during a time interval [0, T ]. The present study is restricted to compressible inviscid flows. The coupling conditions on Σ consist in the continuity of normal velocity components v(x, t).n and the equilibrium of stress vectors σ(x, t).n: ∀(x, t) ∈ Σ× [0, T ]. v f .n f + v s .n s = 0 and σ f .n f + σ s .n s = 0

where x = (x, y), n denotes the outward normal unit vector and f and s refer to the fluid and the structure, respectively. The evolution of the flow is governed by the Euler equations which, in a conservative form, read:

∀(x, t) ∈ Ω f ×]0, T [ ∂ t ρ f + ∇.(ρ f v f ) = 0 (2) ∂ t (ρ f v f ) + ∇.(ρ f v f ⊗ v f ) + ∇P = 0 (3) ∂ t (ρ f E ) + ∇.[(ρ f E + P )v f ] = 0 ( 4
)
where ρ f , P and E are the density, the pressure and the total energy per unit mass. ∇ and ⊗ refer to the nabla operator and the tensor product. The operator ∂ t stands for time derivative. The following equation of state for an ideal gas closes the system:

P = (γ -1)(ρ f E -ρ f v f • v f /2) (5)
The simulations presented in the next sections have been carried out using air as a working gas with a ratio of specific heats γ = 1.4.

As far as elastic structures are concerned, the basic equilibrium equations of linearized elastodynamics in absence of body force are expressed in terms of the displacement field u(x, t) as follows:

∀(x, t) ∈ Ω s ×]0, T [ ρ s ∂ 2 tt u s -∇.σ s = 0 (6)
where ρ s is the density. The second order Cauchy stress tensor σ s is linked to the second order infinitesimal strain tensor ε and so to u s by the following equations:

σ s = C : ε (7) ε = 1 2 ∇u s + (∇u s ) T r (8)
where C denotes the fourth-order elasticity tensor, ":" the double dot product, and () T r the transpose. Lastly, prescribing boundary and initial conditions completes the definitions of the coupled flow evolution and elastic wave propagation problems. Note that, in the previous equations and in the following, the superscripts f and s referring to the fluid and the structure have been omitted for nonconfusing variables. Vectors and tensors are denoted using bold letters.

Partitioned Procedure

The numerical methods used to tackle FSI problems can be broadly classified in two categories [START_REF] Hou | Numerical Methods for Fluid-Structure Interaction -A Review[END_REF]: monolithic approaches and partitioned approaches. In a monolithic approach [START_REF] Hubner | A monolithic approach to fluidstructure interaction using space-time finite elements[END_REF], the fluid and structure subsystems are treated in the same mathematical framework: a unified algorithm solves the single system of equations that governs the entire problem. Interfacial conditions (Eq. ( 1)) are taken into account in an "implicit" way. This approach can achieve high accuracy, but the resulting numerical solver may be hard to develop and to maintain.

On the opposite, in a partitioned procedure [START_REF] Farhat | Provably second-order timeaccurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity[END_REF], the fluid and structure subsystems are solved separately with their respective numerical algorithms and mesh discretization. In this case, the interface conditions are taken into account in a more "explicit" way by exchanging numerical data between the fluid and the structure. They benefit from numerous advantages such as their ease of implementation and software modularity.

The present study falls within the framework of a challenging project dedicated to the simulation and the modeling of wave propagation in both fluid and structure. The underlying physics involves intricate multiphysics and multiscale phenomena in both media such as compressible turbulence and wave propagation in anisotropic heterogeneous materials. From that perspective, relying on sophisticated and constantly enhanced codes to effectively capture the wide range of relevant physical phenomena, while retaining a high level of flexibility, is crucial. Partitioned procedures meet these specifications and have been selected.

Two in-house codes based on well-grounded numerical methods are the key-ingredients of the coupling tool. The fluid is handled by the high-order structured finite differences hydrocode named Choc-Waves (Compressible High-Order Code using Weno AdaptiVE Stencils) developed at CORIA laboratory, France [START_REF] Chaudhuri | Numerical study of compressible mixing layers using high-order WENO schemes[END_REF][START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF]. Large-eddy simulations of shock/boundary layer interaction in transient nozzle flows [START_REF] Chaudhuri | Numerical investigations of transient nozzle flow separation[END_REF] are one example of application of this threedimensional compressible Navier-Stokes solver. A simplified version of the code is employed to solve the Euler equations [START_REF] Depalma | Blast injuries[END_REF][START_REF] Igra | Review of methods to attenuate shock/blast waves[END_REF][START_REF] Berger | Experimental investigation on the shockwave load attenuation by geometrical means[END_REF][START_REF] Chaudhuri | Computational study of shock-wave interaction with solid obstacles using immersed boundary methods[END_REF]. Specifically, a fifth-order WENO (Weighted Essentially Non-Oscillatory) scheme [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] for convective fluxes, combined with a third-order TVD (total variation diminishing) Runge-Kutta method for time advancement, makes it suited for simulations involving steep pressure gradients. The treatment of the fluid mesh is another major classification criterion for FSI solution procedures (the reader may refer to [START_REF] Hou | Numerical Methods for Fluid-Structure Interaction -A Review[END_REF] for a review). This version relies on a body-fitted grid method: the fluid mesh conforms to the structure mesh at the interface and remeshing is needed when addressing moving/deforming boundaries. This is performed by means of a dynamic mesh algorithm [START_REF] Batina | Unsteady euler airfoil solutions using unstructured dynamic meshes[END_REF][START_REF] Piperno | Numerical simulation of fluid-structure interaction phenomena (simulation numérique de phénomènes d'interaction fluide-structure)[END_REF][START_REF] Mouronval | Une méthodologie générale pour l'étude numérique du couplage fluide-structure par une approche multicode[END_REF]. As far as the structure is concerned, the object-oriented finite elements code OOFE (Object Oriented Finite Element program) developed at MSSMat laboratory, France, allows for the simulation of wave propagation in complex structures (composite materials such as honeycomb core sandwich panels [START_REF] Grédé | Numerical modeling of pyrotechnic shock wave propagation in ariane5's structures (Modélisation des chocs d'origine pyrotechnique dans les structures d'ariane5 : développement de modèles de propagation et d'outils de modélisation)[END_REF]). One of the strengthes of this code for linear/nonlinear static/transient structural mechanics lies in its adaptative mesh capacity [START_REF] Tie | Adaptive time discontinuous Galerkin method for numerical modelling of wave propagation in shell and 3d structures[END_REF] and in the diversity of implemented solvers. More precisely, one may solve Eqs. (6-8) using a time discontinuous space-time Galerkin method [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF][START_REF] Johnson | Discontinuous galerkin finite element methods for second order hyperbolic problems[END_REF][START_REF] Li | Implementation and adaptivity of space-time finite element method for structural dynamics[END_REF][START_REF] Tie | Adaptive time discontinuous Galerkin method for numerical modelling of wave propagation in shell and 3d structures[END_REF], a space discontinuous Galerkin method [START_REF] Tie | A unified variational framework for the space discontinuous galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF][START_REF] Käser | An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes i: The two-dimensional isotropic case with external source terms[END_REF][START_REF] Dumbser | An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes ii: The three-dimensional isotropic case[END_REF] or the classical Newmark method.

Staggered Algorithm

The fluid and structure subsystems are time-integrated independently using their own scheme. The choice of a coupling scheme is subjected to the specificities of the studied FSI problem. For instance, the interaction of an incompressible fluid with a thin flexible structure may lead to an "add-mass effect" [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF] and thereby requires special attention regarding the selection of the algorithm.

The most popular coupling algorithm for partitioned procedures is the socalled "CSS" (Conventional Serial Staggered procedure) [START_REF] Piperno | Design of efficient partitioned procedures for the transient solution of aeroelastic problems[END_REF]. This procedure is depicted in Fig. 1 and consists in four steps:

1. Transfer the motion of the structural boundary to the fluid system and update the position of the moving fluid mesh accordingly;

2. Advance the fluid system;

3. Transfer the computed pressure on the interface to the structure;

4. Advance the structural system under this flow induced load.

The CSS algorithm may present deficiencies in terms of temporal accuracy and stability under certain conditions since the interface conditions are not exactly verified at each time step, which induces a lack of momentum and energy conservation at the interface [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems -part II: energy transfer analysis and threedimensional applications[END_REF]. In particular, sub-cycling (i.e. advancing the fluid and the structure with different time steps) can deteriorate both the accuracy and the stability of CSS. This is the reason why many variants, among which generalized CSS [START_REF] Piperno | Design of efficient partitioned procedures for the transient solution of aeroelastic problems[END_REF] and Improved Serial Staggered [START_REF] Lesoinne | Higher-order subiteration-free staggered algorithms for non-linear transient aeroelastic problems[END_REF], have been proposed. The abilities and popularity of these schemes are reflected by the number of research works [START_REF] Cirak | Large-scale fluid-structure interaction simulation of viscoplastic and fracturing thin-shells subjected to shocks and detonations[END_REF] that rely on them. The CSS procedure remains a good candidate, for a class of FSI, when used with caution as the following validation case shows. For this reason, the simulations presented hereafter have been performed using this procedure.

Validation: Supersonic Flutter of a Panel

Supersonic flow M ∞ , P ∞ , ρ f ∞

Fixed wall

Simply supported panel Simulating the aeroelastic response of a flat panel under a supersonic flow is a widely-used test case when validating coupling tools in FSI. In fact, this test case provides a means of estimating both the stability and the accuracy of the underlying coupling scheme.

Still air

Figure 2 is a schematic of the test configuration. The panel has an infinite aspect ratio (i.e. with infinite dimension along the orthogonal direction to the flow). The supersonic flow is characterized by its Mach number M ∞ , density ρ f ∞ =0.4 kg/m 3 and pressure P ∞ =13000 Pa. The panel has a length L s =0.5 m, a uniform thickness b s =1.35 mm, a Young modulus E=77 GPa, a density ρ s =2710 kg/m 3 and a Poisson coefficient ν p =0 and is simply supported at both edges (in terms of plate kinematics, this implies that, on the corresponding boundaries, the displacement degrees are set to zero and that the rotation degrees are not prescribed). Its upper side is submitted to the pressure of the flow and its lower side to the pressure P ∞ of still air. An interaction between the supersonic flow and the panel is initiated by a perturbation of the latter.

The resulting instability problem may be investigated analytically using a shallow shell theory for the structure and the linearized formulation of the fluid problem. Such an analysis reveals that a certain value of the freestream Mach number M ∞ triggers a dynamic instability: the two first modes of the panel coalesce leading to an exponential growth of its displacements. For the present configuration, this critical Mach number is M th cr = 2.11. In order to assess the quality of the coupled tool, a series of numerical simulations is performed by varying M ∞ and the agreement between the obtained critical value M num cr and its analytical counterpart is evaluated. For this essentially two-dimensional problem, the panel is modeled as a beam. Its discretization consists in 50 beam elements of Bernoulli type. The implicit Newmark solver (average acceleration method) is selected in OOFE. The first mode of the structure is used as an initial condition. As far as the fluid is concerned, time advancement is performed using the explicit Runge-Kutta solver with a CFL number of 0.8. The fluid mesh consists in 100×126 cells and coincides with the structural mesh at the fluid-structure (FS) interface. Slip boundary conditions are prescribed to the fluid along the fixed wall and the flexible structure alike. Supersonic conditions are used at the fluid inlet and freestream conditions are prescribed at the oulet. The coupling is performed using the CSS procedure described in Section 2.3. The fluid time step is employed in the structure solver so that no sub-cycling is used. Note that more details concerning the numerical set-up can be found in [START_REF] Mouronval | Une méthodologie générale pour l'étude numérique du couplage fluide-structure par une approche multicode[END_REF]: the present simulation only differs by the code used for the structure.

Figure 3 demonstrates the time evolution of the dimensionless vertical displacement u s y /b s at the middle of the beam for four values of M ∞ . An increase followed by a slow decrease is observed for M ∞ = 2.17 and 2.18, whereas a continuous growth is seen for M ∞ = 2.19 and 2.20. Moreover, the growth coefficient is higher for M ∞ = 2.20 than 2.19. This indicates that M num cr ≈ 2.19 with an error of less than 4% with respect to M th cr .

The results of this classical test problem demonstrate the good capacities of the coupling tool and its correct implementation. CSS with no sub-cycling can be adequate to predict physical instabilities in FSI. Nevertheless, this test case does not demonstrate that the coupling tool will be appropriate for any class of FSI. For example, numerical simulations of shock/elastic obstacles implicate some numerical difficulties that differ significantly. Indeed, unlike in the panel flutter simulation for which the pressure evolution along the FS interface is smooth, shock/elastic obstacles simulations involve large gradients of the loading at the FS interface (because of discontinuities in pressure profiles at the FS interface) which require special attention, as the next section shows.

0 0.1 0.2 0.3 t(s) -0.2 -0.1 0 0.1 0.2 u y s /b s M ∞ =2.17 M ∞ =2.18 M ∞ =2.19 M ∞ =2.20

Preliminary Analysis of the Impact of a Shock on a Column of Elastic Obstacles

We now address the issue of a planar shock wave impacting a column of infinite square section cylinders as sketched in two dimensions in Fig. 4 (a). The present section lays the groundwork for the more complex simulations presented in the next section: we review in detail the different kinds of waves that propagate in the fluid and in the structures, when the deformations of the latters are not taken into account by the fluid solver. Furthemore, we highlight some numerical difficulties and workarounds. These preliminary considerations will be very valuable for the simulations reported in Section 4, in which the flow is submitted to some imposed or induced displacements and deformations of the structures.

Studied Configuration and Loading Key-Parameters

Although only three obstacles are depicted in Fig. 4 (a), the number of obstacles in the y direction is assumed to be infinite. The length of the cylinders is chosen large enough compared to the side length of their section so that threedimensional phenomena can be neglected allowing two-dimensional simulations to be carried out.

A rectangular section of size The structure is clamped around its corners.

L f x × L f y ≈ 62 × 5.
Points S 4 , S 5 , F 1 and F 2 are used for postprocessing.

Point O is the axes origin. The compressible fluid is air and the flow is considered to be governed by the Euler equations. At t=0 s, the planar shock wave is located at the leading edge of the column and propagates at a Mach number M shock equal to 1.4 (the corresponding speed, ||v f shock ||, is 485 m/s). This abrupt change in all flow quantities (pressure, density, velocity and temperature) separates the stagnant gas (referred to as "right state" and denoted by the subscript r) from the shocked gas ("left state", denoted by the subscript l hereafter). Providing M shock and the definition of the right state, the Rankine-Hugoniot relation [START_REF] Anderson | Modern Compressible Flow with Historical Perspective[END_REF] enables to fully specify the left state.

The small displacement and deformation theory and the plain strains assumption are retained for the flexible cylinders. At this point, it should be underlined that the speeds (c s p and c s s ) of the pressure waves ("P-waves") and shear waves ("S-waves") in a structure depend solely on its material properties (ρ s , E and ν p ) and may be computed according to:

c s p = E(1 -ν p ) ρ s (1 + ν p )(1 -2ν p ) , c s s = E 2ρ s (1 + ν p ) (9)
This shows that the stiffer the structure is, the faster P and S-waves will propagate. As mentioned before, we intend to illustrate the prominent influence of the material properties on the shock-induced wave fronts that develop in the elastic structures. Hence, two significantly different materials were selected and the simulations were performed using aluminum and a softer material (denoted by mat2 in the following) whose properties and associated wave speeds are reported in Table 1.

Aluminium Mat2 E (MPa) 70×10 iii) m p < m s < 1: subsonic loading (the loading moves along the structure boundary slower than P and S-waves).

The obtained values of these ratios (see Table 1) clearly indicate that the incident shock wave acts as a supersonic loading for the cylinders made of mat2 and a subsonic loading for those made of aluminium. The numerical scheme selected in Choc-Waves is unchanged. As far as OOFE is concerned, two schemes will be tested: the implicit space-time Galerkin method (tDG) and the explicit space discontinuous Galerkin solver (sDG). In sDG, in contrast to continuous finite element methods, the approximated solution is discontinuous from one element to another and a numerical flux is introduced in order to "link" adjacent elements together. tDG relies on a subdivision of the space-time domain Ω s ×]0, T [ into a series of space-time slabs Ω s ×]t n , t n+1 [, with t n+1 = t n + ∆t s and ∆t s the time step of the structure. In each space-time slab, a variational formulation for the displacement and velocity fields is written simultaneously in space and time. Within each space-time slab, continuous finite elements are used, but between two successive slabs (i.e. at t n and t n+1 ), both displacement and velocity fields are discontinuous. The average acceleration Newmark method employed in Section 2.4 is not selected since its lack of damping makes it unsuitable for loadings with strong gradients.

Simulations were carried out according to a specific procedure in order to analyze the different kinds of waves in the fluid and in the structures and to tackle numerical difficulties at the same time. First, we investigated the flow patterns triggered by the impact of the shock on rigid structures (Section 3.2). Next, we performed simulations based on "one-way coupling" where the deformations of the structures are not taken into account by the fluid solver (Section 3.3). These simulations are carried out using structures made of aluminum and the softer material mat2 to put into evidence the influence of the material properties on the pressure and shear wave fronts that propagate in the structures.

Shock/Rigid Structures Interactions

Understanding the complex flow patterns generated by the impact of a planar shock wave on rigid structures is an essential step.

Such a flow involves multiple shock waves and discontinuities, requiring the use of a low-dissipation, high-order shock-capturing scheme. WENO (Weighted Essentially non-oscillatory) schemes are a good choice to avoid numerical damping and to prevent spurious numerical oscillations near discontinuities. Indeed, these schemes employ a superposition of several sub-stencils with adaptive coefficients to construct a high-order approximation of the solution. This procedure avoids interpolation across discontinuities and preserves uniformly high-order accuracy at all points where the solution is smooth. In the present study, we use a finite-difference, flux-based, fifth-order WENO scheme [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] with the Lax-Friedrichs splitting for calculating the numerical fluxes at cell interfaces. Three sub-stencils, each containing three grid points, are combined to achieve fifthorder accuracy. The smoothness indicators that measure the smoothness of the flux function on each candidate sub-stencil are computed owing to [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF]. A more detailed description of the method can be found in [START_REF] Chaudhuri | Computational study of shock-wave interaction with solid obstacles using immersed boundary methods[END_REF]. The reader may refer to [START_REF] Mouronval | Numerical Study of the Starting Process in a Supersonic Nozzle[END_REF] for a validation and an illustration of the capabilities of the adopted solver on the starting process in a supersonic nozzle.

335

In the following, the cylinders are supposed rigid and no rigid body motion is allowed. Simulations rely solely on Choc-Waves. The duration of the simulation is equal to 60 µs. The computational domain is composed of 400 × 100 cells and is refined around the structure (see Appendix B for a mesh refinement study). The CFL number is set to 0.7.

Chaudhuri et al. [START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF] recently conducted similar simulations to clarify the role of some geometrical parameters (such as the shape of the structures) in terms of shock attenuation. These simulations constitute a foundation on which to analyze the complex flow patterns encountered hereafter. Figures 5 (a) -5 (c) present numerical schlieren pictures obtained from the density gradient field [START_REF] Hadjadj | Computation and flow visualization in highspeed aerodynamics[END_REF]. The symmetries of the computational domain were exploited during the post-processing step in order to partially reproduce the column and to facilitate understanding the flow dynamics. At t=0 s, the incident shock wave (IS) is located at the leading edge of the column. Its propagation along the structure wall gives rise to a regular reflection (RS). Next, part of the IS (that is to say the transmitted shock TS) passes through the open passage between the two structures and generates vortices (V) [START_REF] Sun | Vorticity production in shock diffraction[END_REF] with shocklets and others reflected waves (RS2). As Chaudhuri et al. [START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF] pointed out, the interaction between transverse waves and vortices produces a secondary shock (SS).

Note that the profiles of the transmitted and reflected shocks are less sharp as they propagate since the mesh is progressively coarsened away from the structures. Figures 5 (d) -5 (f) depict the overpressure ratio profiles (P -P r )/P r recorded at the structure wall (boundary S 1 S 2 at -4.4 × 10 -3 < x < 0 and S 2 S 3 at 0 < x < 4.4 × 10 -3 ). From these profiles, four points demonstrate the complexity of the loading on the structure:

1. The spatial evolution of the pressure loading exerted on the structures is characterized by strong gradients. This is especially true in the very early stages of the impact, when the incident shock moves along S 1 S 2 ; 2. The reflected wave (RS) induces a clear rise in pressure along S 1 S 2 : the overpressure ratio behind the incident shock increases from (P l -P r )/P r =

1.1 at t=0 s to more than 3 at t=5.8 µs (Fig. 5 (d)); 3. The induced waves move either downstream (IS, TS...) or upstream (RS, RS2...); 4. The "overpressure" ratio (P -P r )/P r can become negative close to S 2 .

Comparing the present results to those of the viscous simulations by Chaudhuri et al. [START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF] for a single column of structures indicates that the dominant wave pattern is adequately taken into account by non-viscous computations in these early stages of interaction.

One-way Coupling: Fluid → Structure

In the following, the influence of the fluid loading on the structure is investigated using a one-way coupling as depicted in Fig. 6. As alluded earlier, the induced displacements of the structure are not taken into account at this point. This implies that the pressure field in the fluid during the whole simulation will be the same as the one obtained for a rigid structure (Section 3.2). The goals pursued are twofold: comparing supersonic and subsonic loadings by changing the material properties of the structure and selecting an appropriate solver for the structure in presence of strong gradients. About 11000 Lagrange elements of first order (T3 type) are used to discretize the structure. Both sides of the square contain 100 elements, giving a characteristic size of the mesh h s close to 62 µm. The computational domain of the fluid is unchanged (a mesh refinement study is presented in Appendix B). The overpressure P -P r at the FS interface is considered when computing the loadings on the interface.

Structure Made of Mat2

When using the soft material mat2, the loading is supersonic (m p = 2.1 and m s = 3.1): P and S-waves propagate slower in the structure than the incident shock moving along the FS interface. As a result, their fronts form Mach cones whose angles α s p and α s s are given analytically by:

α s p = arctan(1/ m 2 p -1) , α s s = arctan(1/ m 2 s -1) (10)
Using the initial amplitude of v f shock , we get α s p = 27 • and α s s = 19 • . The speed of the incident shock varies slowly in the early stages of the impact and those angles will remain quasi-constant during this period of time.

Figures 7 (a) and 7 (b) display the pressure field in the fluid as well as the S-waves in the structure. These waves are visualized using the rotational of displacements ||rot(u s )||, while P-waves are revealed by their divergence, owing to the Helmholtz decomposition. P-waves are not shown here because they propagate in a similar way as S-waves. These results were obtained with the sDG solver. No sub-cycling is used (both codes rely on the same time step and the loading is transferred at the end of each step, Fig. 6 (a)). Further details regarding the numerical implementation will be given and discussed below. Note that the structure is clamped in the vicinity of its corners (see Fig. 4 (b)) which causes gradients of displacements in this region. At t=8 µs, the incident shock (IS) has almost reached S 2 . Its position at the interface coincides with the tip of the Mach cone (MC) in the structure. The cone is characterized by a 19 • angle, which is in perfect agreement with the analytical solution. At t=23 µs, the transmitted shock (TS) and its associated Mach cone are located close to S 3 . The intricate wave patterns in the structure are directly linked not only to S-waves reflections inside the structure but also to shock waves reflections in the fluid. For instance, one may observe that the reflected shock wave (RS2) is responsible for another Mach cone (MC2) that moves towards S 1 . ). These points start moving at t ≈4.5 µs and t ≈17 µs when the incident (or transmitted) shock wave reaches them. One may also notice that the maximum of the induced displacements is larger for point S 4 than S 5 since the incident shock wave is stronger than the transmitted one (see Figs. 5 (d) -5 (e)). However, these displacements remain small compared to the edge size H s of the structure: ||u s || max /H s is equal to 0.33% at point S 4 and 0.24% at point S 5 .

Considering the second goal of these simulations, namely to select a robust and accurate solver in OOFE to deal with large gradients of the loading, computations were carried out using two candidates: sDG and tDG. Before getting into the issue itself, let us remind the time integrators that are combined with the WENO (Choc-Waves), sDG and tDG (OOFE) schemes, respectively. tDG is implicit and unconditionally stable. Classical linear elements are used in space and one linear element is used in time, the scheme is third-order accurate in time and second-order in space for displacements [START_REF] Johnson | Discontinuous galerkin finite element methods for second order hyperbolic problems[END_REF]. Explicit third-order Runge-Kutta methods are employed for sDG and WENO. They are subjected to a Courant-Friedrichs-Lewy (CFL) stability condition. For the fluid, this con-dition reads CF L f < 1. In contrast, except in some specific cases, there is no analytical formula to predict the stability limit of conditionally stable Runge-Kutta methods when used with sDG [START_REF] Toulorge | CFL conditions for Runge-Kutta discontinuous Galerkin methods on triangular grids[END_REF]. As a matter of fact, the maximum value of the CFL depends on the shape of elements (T3, Q4...), the degree of the polynomials used for space discretization, and the geometry (one or multiple spatial dimensions) and is often determined empirically for a given computation. In the present case, numerical tests demonstrate that this maximum value is close to 0.23.

When no coupling is performed, the time step of the fluid ∆t f is prescribed in Choc-Waves owing to:

∆t f = CF L f max cells (v f x + c f )/∆ f x + (v f y + c f )/∆ f y ( 11
)
where c f is the speed of sound in the fluid, ∆ f x and ∆ f y are the sizes of the cells in x and y directions, respectively. ∆t f varies during the simulation. In practical, we perform simulations using CF L f = 0.7.

On the contrary, the time step used to advance the structure (∆t s ) in absence of coupling is chosen constant. For sDG, ∆t s is set according to:

∆t s = CF L s h s min /c s max (12)
where c s max is the maximum speed of the waves in the structure (here c s max = c s p ) and h s min is the minimum element size of the mesh. For tDG, there is no time step restriction in terms of stability. However, ∆t s must be small enough so as not to damp the high frequencies of interest. In practice, we use ∆t s = 2 h s min /c s max . An analogy with Eq. ( 12) would give an equivalent CF L s of 2.

A comparison of the ratios ∆t s /∆t f where ∆t s are computed from Eq. ( 12) with CF L s = 0.23 for sDG and CF L s = 2 for tDG and ∆t f is recorded during the simulation of the shock/rigid structures impact (Section 3.2) indicates that, most of the time, ∆t s /∆t f is close to 50 for tDG and 7 for sDG. Such values are to be related with the material properties of the structure: the softer the material is, the smaller c s p is. When faced with selecting a staggered scheme to couple two codes, beginners in the FSI field are often drawn to sub-cycling free algorithms because of their ease of implementation. In this instance, one requires the fluid and the structure solvers to advance with the same time step (i.e. ∆t = min(∆t s , ∆t f )). When the forces exerted by the fluid on the FS interface exhibit strong gradients, it is relevant to examine whether enforcing ∆t s = ∆t f will affect the ability of the structure solver to handle such gradients. Figures 8 (a) and 8 (b) display velocity fields, showing the Mach cones associated with both P and S-waves, and velocity profiles for two simulations using sDG and tDG without sub-cycling: OOFE advances using ∆t f and receives data from Choc-Waves at each time step (Fig. 6 (a)). The spurious oscillations that can be observed for tDG are evidence of an improper use of this solver. tDG may handle strong gradients when an appropriate amount of damping is introduced. Specifically, damping in tDG is due to a dissipation of energy that is related with the jump terms in time. These jumps terms are reduced when ∆t s is decreased. As far as sDG is concerned, the velocity field is well-described. This may be explained not only by the fact that ∆t f is not as smaller as the usual range of ∆t s for sDG as for tDG but also by the fact that damping in sDG is related to the discontinuous nature of the approximated solution between elements.

At this point, taking advantage of sub-cycling to use tDG in a more appropriate way seems a suitable alternative. Simulations using tDG with ∆t s = 2h s min /c s max have been performed. Figure 8 (c) presents the results. As expected, no spurious oscillations appear. However, some localized variations in speed, which at first sight look like some low frequency oscillations, are visible beyond the tip of the Mach cones along the FS interface. This is due to the fact that, in the early stages of the simulation, the shock crosses approximately 5 elements during ∆t s . Updating the pressure loadings every ∆t s does not provide the structure with a correct description of the shock motion: the temporal resolution of the boundary condition becomes too coarse when sub-cycling is used, which results in an unsurprising deterioration of the quality of the solution in the structure. It is worth noting that the pressure loadings that are sent to the structure correspond to values at the end of the fluid time steps required to reach ∆t s . Another choice may be to average the pressure field over the fluid time steps [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems -part II: energy transfer analysis and threedimensional applications[END_REF] but this may lead to a smoothing of the gradients induced by the shock displacement.

In the following, we will use sDG with no sub-cycling even if this is not an optimal choice in terms of computational efficiency. For this basic two-dimensional configuration, efficiency is not crucial and designing or testing more complex partitioned procedures is not in the scope of the present work. The previous results clearly demonstrate some of the difficulties one may encounter when using a staggered procedure and the need to continue research studies on this topic.

They also illustrate how convenient a partitioned procedure is compared to a monolithic one. Relying on dedicated codes that are well established in their respective domains and are the fruits of a long-term numerical development efforts makes overcoming numerical difficulties easier. This is largely due to the suitability of the specialized schemes they include.

Structure Made of Aluminium

For this case, the loading is subsonic (m p = 0.12 and m s = 0.19): P and S-waves propagate faster in the structure than the incident shock moving along the FS interface.

Plotting the ratios ∆t s /∆t f where ∆t s are computed from Eq. ( 12) with CF L s = 0.23 for sDG and ∆t f is recorded during the simulation of the shock/rigid structures impact (Section 3.2) shows that ∆t s is always smaller than ∆t f which is not usual in FSI simulations. Indeed, most of the time, ∆t s /∆t f ≈ 0.35. Using the coupling algorithm depicted in Fig. 6 without sub-cycling, requires advancing the fluid using the time step of the structure. In this case, the corresponding CF L f value for the fluid is close to 0. reducing the time step of the fluid does not affect the results of the simulation by comparing the overpressure ratio profiles (P -P r )/P r obtained using CF L f = 0.7 and CF L f ≈ 0.25 along the bottom boundary of the fluid domain (y = 0) at the end of the simulation.

Comparing Fig. 7 with Fig. 9 demonstrates the difference between a supersonic and a subsonic loading. For the structure made of aluminium, no Mach cone is observed since P and S-waves propagate faster than the shocks in the fluid. Point S 4 starts moving at t ≈1.5 µs and S 5 at t ≈2.1 µs. Indeed, the corners of the structure are clamped along each side on a distance equal to H s /10 so the incident shock triggers a displacement of the structure after t = 0.1H s /|(v f shock .t)| ≈0.9 µs. The resulting P-waves require approximately 0.5 µs and 1.2 µs to reach S 4 and S 5 . Interestingly, point S 5 starts moving approximately 16 µs before the transmitted shock wave reaches it. In the same way, Fig. 9 (a) shows S-waves around S 3 whereas the incident shock is close to S 2 . In theory, this infers that in the case of a two-way coupling, the fluid region located close to the edge S 2 S 3 will be informed of the impact of the structure before the arrival of the transmitted shock wave, via the edge movement. However, these displacements are very small compared to H s (||u s || max /H s is equal to 3×10 -4 % for point S 4 and 2×10 -4 % for point S 5 ) and are unlikely to affect the fluid in case of a two-way coupling, as Section 4.2 will confirm. The results reported in this section provide a sound basis for carrying out more complex simulations. Indeed, various wave patterns have been observed and explained in both the fluid and the structures in the early stages of the impact. The evolution of the flow (reflected and transmitted shocks, vortices...) is in good agreement with the results of Chaudhuri et al. [START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF], which confirms the abilities of the WENO solver. The angles of the Mach cones for both P and S-waves in the case of mat2 perfectly match their analytical values and sDG has been shown to be well-suited for such simulations.

Influence of the Elastic Obstacles Deformations on the Shock Propagation at Early Stages of the Impact

The present section relies on the previous preparatory study and aims at enriching the investigation of shock mitigation by solid obstacles by taking into account the displacements and the deformations of the structures in the fluid.

One-way Coupling: Structure → Fluid

In this first part, we investigate the influence of a prescribed motion of the structure on the reflected and transmitted shock waves (RS and TS) in the fluid. We perform a one-way coupling where the structure transfers the position of the FS interface to the fluid but does not receive any pressure loading from it (Fig. 10).

Structure Fluid

Without sub-cycling An initial velocity is prescribed to the structure in such a way that the induced motion of the structure corresponds to a free eigenmode vibration. Its shape and maximum amplitude (||u s || max /H s = 1.28%) are the same for both materials, however the corresponding eigenfrequencies F s and maximum velocities (||v s || max = 2πF s ||u s || max ) differ (see Table 2).

Deformed shape

Al. Mat2 The incident shock wave requires T f shock ≈23 µs to travel, in the fluid, from point S 1 to point S 3 . Computing the ratio between this duration and the period of the eigenmode T s = (F s ) -1 indicates that the structure will vibrate slowly compared to the shock motion for the case mat2 and fastly for aluminium. This is confirmed by preliminary simulations that highlight the fact that the prescribed motion of the structure made of mat2 induces very small pressure variations in the fluid, contrary to the structure made of aluminium. So, for the sake of brievety, we will focus on the structure made of aluminium.

We are interested in evaluating the influence of this prescribed motion on the intensity and velocity of the shock waves (RS and TS) in the fluid in the vicinity of the vibrating structure only. Therefore, we decided to stop the simulation when the transmitted shock TS has reached the point located at (x, y) ≈(11.0,0) mm = (1.25D s , 0) (Fig. 14 and Fig. 15), which corresponds to t=31.2 µs. The length of the fluid computational domain was shortened (L f

x ≈ 32 mm) accordingly. A mesh refinement study was performed to make sure to capture accurately the pressure variations induced by the structure motion (Appendix B). As a result, the number of fluid was doubled in each direction, leading to a total number of cells of N f

x × N f y = 800 × 200. The structure mesh is also refined to be coincident (200 elements along each edge side are used giving four times as many elements as in Section 3.3). These mesh refinements induce a decrease of both ∆t s and ∆t f . ∆t s is equal to 1.27×10 -9 s and ∆t s /∆t f = 0.5 approximately. The coupling algorithm is not sub-cycled (both codes are advanced using ∆t = min(∆t f , ∆t s )). Evaluating the influence of the structure motion on RS and TS requires comparisons with the case of a rigid structure (Section 3.2). Assessment is made using the same meshes and time steps (some of the previous simulations were rerun). The eigenmode prescribed to the structure is displayed in Table 2. At t = 0 s, the corresponding velocity field is applied to the structure and the initial displacement field is set to zero. Each edge (S 1 S 2 and S 2 S 3 ) adopts a periodic "W"-like shape motion. Figure 11 (a) illustrates the normal displacement along S 1 S 2 at time t = T s /4, when the amplitude of the displacement along this interface reaches its maximum value. The normal displacement along S 2 S 3 is identical and thus not presented. Note that the corners of the structure are clamped along each side on a distance equal to H s /10. A curve fitting has been performed, so that one may carry out the same simulation without computing the eigenmode. The obtained function u s,f it (s, t).n s is given in Appendix A. As can be seen, the maximum normal displacement at point S 4 is about 8×10 -5 m which corresponds to 1.28%H s . In comparison, the maximum value obtained at point S 4 due to the shock impact (see Section 3.3.2) is more than 4000 times smaller (see Fig. 9). Figure 11 (b) presents the temporal evolution of the displacement at point S 4 at the end of the simulation. One can check that the initial motion has almost not been damped by the numerical scheme in the structure. This 'W"-like shape motion of the edges of the structure triggers expansion and compression waves in the fluid. Figures 12 (a) and (b) display the temporal pressure evolution at points S 4 and S 5 . Three cases are compared: i) the present case ("S→F"), where the incident shock impacts the vibrating structure;

ii) the rigid one, where the shock impacts a rigid structure;

iii) a third one ("S→F without shock"), where the vibrating structure is placed in a flow at rest (this corresponds to the configuration depicted on Fig. 4 with the left state set equal to the right state).

In the "S→F without shock" case, the pressure evolution is the same for both points, as expected. The intensity of the expansion and compression waves is a function, among others, of the local speed of sound in the fluid, whose value is modified by the waves themselves. The strong amplitude of the compression waves, which can reach values up to 6P r (P r being the ambiant pressure), is to be related with the high speed of the structure motion (||v s || max =465 m/s).

In the "S→F" case, the incident shock reflects on the vibrating structure and reaches point S 4 around t=5 µs (RS), leading to an increase in pressure. As far as point S 5 is concerned, a smaller increase appears at t=15 µs, slightly before the arrival time of the transmitted shock (TS) in the rigid case, indicating that the shock has been accelerated by the vibration. Next, the pressure at S 5 decreases compared to "S→F without shock", due to complex interactions between several shocks, some low pressure zones with vortices (see Figs. Comparing more precisely these pressure evolutions between the "S→F" case and its rigid counterpart is achieved by computing the relative difference, that is, 100 × (P (S→F)-P (rigid))/P (rigid), at both points. Figure 13 presents the results. After the shock arrival (i.e. around t=5 µs and t=15 µs for S 4 and S 5 , respectively), the pressure increase is, on average, of the order of 150%. The high values obtained for point S 5 , compared to point S 4 , between t=5 µs and t=15 µs, are due to the low values of P (rigid) at point S 5 during this time (P (rigid) = P r ). In order to have a comprehensive view of the complex interactions between the prescribed vibration and the different shocks in the fluid, the pressure gradient fields at the end of the simulation have been computed for both cases (Fig. 14). The compression waves due to the vibration of the structure are clearly visible, as is the small shock wave (PS) associated with the onset of the prescribed vibration. Several shocks whose extremities were propagating along the interface in the rigid case are now detached from it. Both the reflected and the transmitted shocks have been accelerated by the vibration. Moreover, the shape of the transmitted shock has been significantly modified. This accelaration is an interesting result: even if in the context of shock attenuation one desires to reduce the speed and the strength of the reflected and the transmitted shock waves, this suggests that one may actively control them by prescribing a motion to the structure. Analyzing the pressure profiles along the lower and the upper boundaries of the fluid domain (lines y = 0 and line y = y max (x), see Fig. 4 (b)), for both the "S→F" case and its rigid counterpart, clarifies the reasons for this change in shape of the transmitted shock. It is recalled that the fluid mesh conforms to the structure mesh at the FS interface: y = y max (x) coincides with this interface for -4.4×10 -3 < x < 4.4×10 -3 . Figure 15 (a) presents the results at the end of the simulation for y = 0. Both RS and TS have been accelerated by the vibration and the pressure oscillations due to the prescribed oscillations of the structure can be observed. As far as the upper boundary, y = y max (x), is concerned, Fig. 15 (b) puts into evidence the strong pressure variations at the FS interface, as already revealed by the pressure evolutions recorded for points S 4 and S 5 , as well as a stronger increase in speed for TS as compared to y = 0, which explains its change of shape. Computing the corresponding relative differences, 100×(P (S→F)-P (rigid))/P (rigid), confirms these conclusions (figure not shown for the sake of brievety).

P (Pa) Temporal evolution of pressure P at points F 1 and F 2 (a) for the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F") and in the case of a rigid structure placed in the same flow ("rigid"). Temporal evolution of the relative differences 100 × (P (S→F)-P (rigid))/P (rigid) at points F 1 and F 2 , in the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F")(b).

Assessing more accurately the amplification of the transmitted shock can be achieved by computing pressure impulses for the present case and the rigid one. This requires recording pressure histories at some probes situated in the fluid initially at rest, downstream of the obstacles. For this purpose, two probes located close to the obstacles at (x, y) = (0.5D s , 0) and (x, y) = (0.75D s , 0) and denoted by F 1 and F 2 on Fig. 4 (b) and Fig. 15 (a) have been used. Figure 16 (a) depicts the obtained pressure signals. The corresponding pressure impulses have been computed owing to [START_REF] Berger | Experimental investigation on the shockwave load attenuation by geometrical means[END_REF][START_REF] Chaudhuri | Numerical study of shock-wave mitigation through matrices of solid obstacles[END_REF]:

I (t c ) = tc t0 P dt (13)
with t 0 the time at which the shock reaches the probe, determined by an increase of 0.1% in pressure (i.e. P > 1.001P r with P r the ambiant pressure). Note that the shock due to the onset of the vibration (PS) is taken into account when computing I (t c ) since it generates an overpressure that may damage buildings in the same way as the transmitted shock does. The impulse integration time t c is chosen to be the final time (i.e. t c =31.2 µs). The values of t 0 and I (t c ) for the two probes are reported in Table 3: the impulse rises by 46% for point F 1 and 65% for point F 2 , due to the prescribed vibration of the obstacles. The relavite differences between the present case and the rigid one, 100 × (P (S→F)-P (rigid))/P (rigid) are plotted on Fig. 16 (b) and put into evidence the fact that the wave associated with the onset of the structure perturbation (PS) results in a non-negligible increase in pressure (up to 50%), which contributes significantly to the increase of the pressure impulse at both points. To conclude, let us mention that we have chosen to focus on pressure impulses downstream of the obstacles where the buildings and the population are supposed to be located, for the sake of concision. The same type of analysis can be conducted for the reflected shock (RS). The previous sections highlight the influence of the fluid on the elastic structure (F→S simulations) and reciproquely (S→F simulations). In this section, the fully-coupled behavior of the system is investigated. More precisely, two-way coupled simulations based on the algorithm depicted in Fig. 1, without subcycling, are performed to estimate the degree of coupling of the FS interaction for the case of the structure made of aluminium without any initial perturbation of the structure and next with the same initial eigenmode perturbation as in Section 4.1.

In order to make comparisons with one-way coupled simulations, the following notations have been adopted: u s refers to the present two-way coupling, with or without initial perturbation owing to the case, whereas u s (S→F) and u s (F→S) refer respectively to the solution obtained for the S→F one waycoupling (Section 4.1) and the F→S one way-coupling (Section 3.3). The same notations are used for the pressure field (P , P (S→F) and P (F→S)). Comparisons are made using the same meshes, duration and time steps as in Section 4.1.

Two-way Coupling Without Initial Perturbation of the Structure

In this case, the incident shock impacts the structure at rest and triggers displacements that may give rise to a feedback effect by modifying the pressure field in the fluid. Comparing the obtained results with those of the F→S simulation makes it possible to evaluate the degree of interaction between the fluid and the structure.

The temporal pressure evolutions at points S 4 and F 1 (Fig. 17 (a) and Fig. 17 (c)) give a first indication of the degree of interaction: the curves related to both types of coupling coincide (the results for points S 5 and F 2 are similar and have been omitted to facilitate the reading of the curves). The corresponding relative differences, 100 × (P -P (F→S))/P (F→S), do not exceed 0.06 (Fig. rise to a feedback. No additional loadings act on the structure when the system is fully-coupled and its displacements remain unchanged too (Fig. 18 (a)). Even if this difference increases with time, it remains lower than 0.008, which confirms the absence of any feedback phenomenon. The conclusion of this first two-way coupled simulation is clear: ignoring the influence of structural displacements on the fluid flow is relevant; a F→S simulation describes both the fluid and structure evolutions in an accurate manner when the structure has not undergone any initial perturbation.

Two-way

Coupling With an Initial Eigenmode Perturbation of the Structure Before addressing in detail the results of this second two-way coupled simulation, we have found helpful to review the underlying ideas of the preceding one-way coupled simulations: i) F→S coupling neglects the effects of structure displacements on the fluid.

This implies that the displacements u s (F→S) result exclusively from the pressure loading associated with the shock impact. The pressure field P (F→S) is the same as the one obtained when the shock impacts a rigid structure (Section 3.2).

ii) S→F coupling neglects the effects of the fluid on the structure. This implies that the displacements u s (S→F) result exclusively from the initial perturbation of the structure. It is recalled that they are about 4000 times larger than u s (F→S) for the chosen perturbation. The pressure field P (S→F) differs from P (F→S) solely by the significant pressure variations induced by the eigenmode vibration of the structure.

In this second F↔S simulation, the structure is perturbed using the same initial velocity as in Section 4.1 and is subjected to the pressure loadings exerted by the fluid, whose evolution is affected by the induced displacements of the structure. Our goal is to evaluate the structure displacements that arise only from the interactions with the fluid, that is, u su s (S→ F). Three types of pressure loadings contribute to these displacements: i) P l impact loadings related to the pressure field P (F→S) associated with the shock impact (Sections 3.2 and 3.3, see also P ("rigid") on Fig. 12);

ii) P l eigenmode loadings related to the pressure variations P (S→F) -P (F→S) induced by the eigenmode motion of the structure (Section 4.1);

iii) P l f eedback loadings related to the pressure variations P -P (S→F) induced by the displacements due to P l impact and P l eigenmode (i.e. feedback effect). Figure 19 demonstrates that the pressure variations associated with P l f eedback are negligible with respect to P (S→F): P -P (S→F) remains lower than 2%P (S→F) for point S 4 and point F 1 (once again the results for points S 5 and F 2 are similar and have been omitted to facilitate the reading of the curves). These pressure variations associated with the feedback are higher than those previously obtained when the structure was not perturbed (Fig. 17) but are still very low.

Moreover, they can also be neglected with respect to P (S→F) -P (F→S), as Fig. 20 shows, which implies P l eigenmode P l f eedback . At this point, it is especially suited to compare three structure displacements: ı) ||u su s (S→ F)||, that is, the displacements triggered by P l impact + P l eigenmode + P l f eedback ; ii) ||u s -u s (S→ F)-u s (F→S)||, the displacements generated by P l eigenmode + P l f eedback ; iii) ||u s (F→S)||, the displacements associated with P l impact .

Given that P l eigenmode +P l f eedback ≈ P l eigenmode , ||u s -u s (S→ F)-u s (F→S)|| is an accurate estimate of the displacements generated by the pressure variations associated with the eigenmode motion of the structure. Figure 21 ( In summary, there is not any noticeable feedback effect for this configuration either. Both the pressure field in the fluid and the displacements in the initially perturbed structure can be described with a good accuracy by conducting a S→F simulation.

It can be concluded from the previous fully-coupled simulations (without and with an initial perturbation of the structure), that, in absence of a prescribed motion of the structure, the degree of coupling for this configuration and this material is low: in this case, carrying out simulations with structures considered as rigid to evaluate shock wave attenuation is relevant for the present configuration. When an initial displacement or velocity is prescribed to the structure, a S→F one-way coupling is required to take into account the resulting decrease of the mitigation process. In neither case is a fully-coupled simulation required.

Conclusion and Future Work

In this paper, a tractable coupling tool based on a partitioned procedure has been presented and applied to investigate shock/elastic obstacles interactions. A step-by-step computational methodology has allowed for identifying the role of each component (fluid, structure and coupled system). The complex flow patterns and the associated pressure loadings have been analyzed. A one-way F→S coupling has put into evidence two very different types of loading (supersonic and subsonic) owing to the material properties of the structures. Furthermore, some numerical difficulties related to the strong gradients of the loadings and the discrepancies between the time steps of the fluid and the structure have been highlighted and have led us to select a space discontinuous Galerkin method for the structure and to avoid using sub-cycled procedures. Next, the acceleration of both the reflected and transmitted waves in the fluid by a prescribed motion of the structures has been exhibited. Pressure impulses downstream of the obstacles have been evaluated to quantify the amplification of the transmitted waves. Finally, two-way coupled simulations have demonstrated that uncoupled simulations using rigid structures may be appropriate to evaluate shock mitigation for the studied configuration when no motion is prescribed to the structures. When the structure is initially perturbed, a S→F one-way coupling is needed to provide a correct description of the mitigation process.

In fact, prescribing a vibration to the structures can significantly modify the flow field. Increasing the number of columns of structures and varying the eigenmodes constitute one of the research topics that we intend to investigate in the future, by means of advanced numerical methods for the fluid, such as Immersed Boundary Methods combined with adaptive multiresolution [START_REF] Soni | On the accuracy and efficiency of pointvalue multiresolution algorithms for solving scalar wave and Euler equations[END_REF].

As far as Fluid → Structure simulations (Section 3.3) are concerned, the following criteria were used to assess whether a grid combination was well-suited or not: -flow patterns (shocks, vortices, slip lines...) must be captured since we need to understand the flow dynamics; -pressure evolutions at the interface must be converged enough so that the induced displacements are not a function of the mesh combination. The pressure evolution and the induced displacements at points S 4 and S 5 , as well as numerical schlierens of the flow field, have been used to compare the results for M 1 , M 2 and M 5 (figures not shown for the sake of brievety). M 2 meets the criteria mentioned above and was used for the simulations of Section 3.

Structure → Fluid simulations (Section 4.1) with structures made of aluminium aim at investigating the influence of a free eigenmode vibration of the structures on the flow and more precisely on the transmitted shock in the vicinity of the obstacles. The amplitude of this free vibration should remain almost constant during the simulation and the induced pressure variations at the FS interface should be well-described. The fluid mesh must be fine enough to reach convergence for the reflected and transmitted shock positions (RS and TS). The evolution of the pressure impulses at F 1 and F 2 with the grid must also be taken into account when selecting a mesh combination. The influence of the six mesh combinations reported in Table B.4 on these quantities has been evaluated by plotting the displacements at point S 4 , the pressure evolution at this point and along y = 0 at the end of the simulation and the temporal pressure evolution at points F 1 and F 2 (figures not shown) and by computing the pressure impulses at F 1 and F 2 (Table B.5). These six mesh combinations provide coherent results: even if the number of cells and elements and the associated time step greatly vary, the obtained results display the same overall behaviour (acceleration of TS and RS, pressure evolution at the FS interface...). When using M s 1 (M 1 combination) and its associated time step, the free vibration of the structure is significantly damped, therefore M 1 can not be selected for Structure → Fluid simulations with structures made of aluminium. M s 2 (which is part of M 2 , M 3 and M 4 ) gives an acceptable prediction of the displacements at S 4 whereas M s 3 (which is part of M 5 and M 6 ) provides a very good estimation. Increasing the number of cells in the y direction, N f y , has an influence on the pressure results and N f y = 200 appears to be a good choice (M 3 and M 4 give similar results). From M 3 , increasing the number of cells in the x direction, N f

x (see combinations M 3 , M 5 and M 6 ), in the obstacles zone and outside, leads to a sufficient convergence in terms of pressure evolution at the FS interface, shock positions along y = 0 and pressure impulses at F 1 and F 2 . For example, the difference between the pressure impulses obtained with M 6 differ by less than 1% from those obtained with M 5 (Table B.4). M 5 could be a relevant choice for Structure → Fluid simulations with structures made of aluminium.

However, the mesh combination that consists of the finest fluid and structure meshes, M 6 , provides a slightly better resolution of the reflected and transmitted shocks along y = 0 at the end of the simulation and numerical schlierens of better quality outside the obstacles zone. These are the reasons why M 6 was selected for Structure → Fluid simulations. The same tests were performed for the two-way coupled simulations of Section 4.2 and M 6 was selected for these simulations too.

  Left state in the fluid (Fig. 4) r Right state (ambiant state in the fluid, Fig. 4) ∞ Freestream value Latin Letters b s Panel thickness (Section 2.4) c s p , c s s Speed of pressure and shear waves, respectively (Eq. (9)) c s max = c s p Maximum speed of the waves propagating in the structure D s

Figure 1 :

 1 Figure 1: Conventional Serial Staggered algorithm. The dashed lines symbolize data transfer between the fluid and structure solvers.

Figure 2 :

 2 Figure 2: Flat panel with infinite aspect ratio in a supersonic air flow.

Figure 3 :

 3 Figure 3: Time evolution of the dimensionless vertical displacement u s y /b s of the point located in the middle of the beam for M∞ = 2.17, 2.18, 2.19, and 2.20.

8 mm 2 1 S 3 S 2 S 4 S 5 F 1

 2132451 is chosen for the computational domain of the fluid (Fig.4 (b)). The lengths of the diagonal and the side of the square are D s =8.8 mm and H s ≈ 6.2 mm, respectively. This results in an effective flow area (open passage) between the obstacles ε = 1.0-0.5D s /L f y Schematic of the studied configuration at initial time S Computational domains for fluid and structure.

Figure 4 :

 4 Figure 4: Two-dimensional schematic of a planar shock wave impacting a column of infinite square section cylinders and computational domains.

  The pressure loading induced by the propagation of the incident shock wave along the FS interface is equivalent to a moving load. Initially, the incident shock wave propagates along the FS interface at speed|(v f shock .t)|,where t is the tangent vector at the interface. Depending on the ratios m p = |(v f shock .t)|/c s p and m s = |(v f shock .t)|/c s s , three cases can be identified: i) m s > m p > 1: supersonic loading (the loading moves along the structure boundary faster than P and S-waves); ii) m s > 1 > m p : transonic loading (the loading moves along the structure boundary faster than S-waves but slower than P-waves);

Figure 5 :

 5 Figure 5: Shock/rigid structures simulations. Numerical Schlierens (a)-(c) and overpressure ratio profiles (P -Pr)/Pr at the structure boundaries S 1 S 2 and S 2 S 3 (d)-(f), at t=5.8 µs (pictures (a) and (d)), t=16 µs (pictures (b) and (e)) and t=25.6 µs (pictures (c) and (f)). S 1 , S 2 and S 3 are located at x=-4.4×10 -3 m, 0 and 4.4×10 -3 m, respectively. IS: incident shock, RS and RS2: reflected shocks, TS: transmitted shock, V: vortex, SS: secondary shock.

aFigure 6 :

 6 Figure6: One-way coupling scheme (F→S). The fluid transfers the pressure loading at the FS interface to the structure. The structure is advanced using these external forces but the induced displacements are not taken into account by the fluid.

5 Figure 7 :

 57 Figure 7: One-way coupling: F→S (mat2). Pressure field in the fluid (min=39 kPa, max=473 kPa) and S-waves in the structure made of mat2 at (a) t=8 µs and (b) t=23 µs. MC and MC2: Mach cones. Temporal evolution of the displacement ||u s || at points S 4 and S 5 (c). See Fig. 4 for points location.

Figure 7 (

 7 Figure 7 (c) presents the displacements ||u s || at the points S 4 and S 5 located in the middle of S 1 S 2 and S 2 S 3 (see Fig. 4 (b)). These points start moving at t ≈4.5 µs and t ≈17 µs when the incident (or transmitted) shock wave reaches them. One may also notice that the maximum of the induced displacements is larger for point S 4 than S 5 since the incident shock wave is stronger than the transmitted one (see Figs. 5 (d) -5 (e)). However, these displacements remain

Figure 8 :

 8 Figure 8: One-way coupling: F→S (mat2). Velocity field ||v s || (min=0 m/s, max=1.07 m/s) at t=9 µs for (a) sDG without sub-cycling, (b) tDG without sub-cycling and (c) tDG with sub-cycling.

5 Figure 9 :

 59 Figure 9: One-way coupling: F→S (aluminium). Pressure field in the fluid (min=39 kPa, max=473 kPa) and S-waves in the structure made of aluminium at (a) t=8 µs and (b) t=23 µs. Temporal evolution of the displacement ||u s || at points S 4 and S 5 (c).

Figure 10 :

 10 Figure 10: One-way coupling scheme (S→F). The structure transfers the displacements of the FS interface to the fluid. The fluid is advanced using the new position of the FS interface. Pressure loads are not taken into account by the structure.

FTable 2 :

 2 s (kHz)T s (µs) T f shock /T s ||v s ||max (m/sOne-way coupling: S→F. Deformed shape (scaled by a factor of 3) of the structure at a given time with amplitude of displacements ||u s ||, eigenfrequency (F s ) and eigenmode period (T s = (F s ) -1 ) of the selected mode, ratio between the time required by the incident shock wave to travel in the fluid from point S 1 to point S 3 (T f shock ≈23 µs) and the eigenmode period (T s ), and maximum velocity ||v s ||max at time t = 0 s.

4 Figure 11 :

 411 Figure 11: One-way coupling: S (aluminium)→F. Normal displacement u s (s, t).n s (with s the curvilinear abscissa) along S 1 S 2 , recorded during the simulation at t = T s /4 and curve corresponding to the fitting function given in Appendix A at the same time (a). Temporal evolution of the norm of the displacement ||u s || at the end of the simulation at point S 4 (b).

Figure 12 :

 12 Figure12: One-way coupling: S (aluminium)→F. Temporal evolution of pressure P at points S 4 (a) and S 5 (b) for the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F") or in a flow at rest ("S→F without shock") and in the case of a rigid structure placed in the flow containing the incident shock ("rigid"). The ambiant pressure, Pr, is reported on the y-axis for comparison.

  5 (b) and (c) in the rigid case) and the vibrating structure.

5 Figure 13 :

 513 Figure13: One-way coupling: S (aluminium)→F. Temporal evolution of the relative differences 100 × (P (S→F)-P (rigid))/P (rigid) at points S 4 and S 5 , in the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F").

Figure 14 :

 14 Figure 14: One-way coupling: S (aluminium)→F. Pressure gradient field (represented in a schlieren-like manner) in the fluid and displacement field (||u s ||) in the structure for the "S→F" case (top) and the rigid one (bottom), at t=31.2 µs. PS: shock associated with the onset of the prescribed vibration.

Figure 15 :

 15 Figure15: One-way coupling: S (aluminium)→F. Temporal evolution of pressure P along y = 0 (a) and y = ymax(x) (b) for the perturbed structure ("S→F") and for a rigid structure, at t=31.2 µs. The structure is located between x=-4.4×10 -3 m and 4.4×10 -3 m. PS: shock associated with the onset of the prescribed vibration. "F 1 " and "F 2 " on the x-axis of (a) refer to the positions of the two probes mentioned on Fig.16.

Figure 16 :

 16 Figure16: One-way coupling: S (aluminium)→F. Temporal evolution of pressure P at points F 1 and F 2 (a) for the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F") and in the case of a rigid structure placed in the same flow ("rigid"). Temporal evolution of the relative differences 100 × (P (S→F)-P (rigid))/P (rigid) at points F 1 and F 2 , in the case of the perturbed structure vibrating in the flow containing the incident shock ("S→F")(b).

17

 17

  (b) and Fig.17

1 TSFigure 17 : 4 Figure 18 :

 117418 Figure 17: Two-way coupling without initial perturbation of the structure: F↔S (aluminium). Temporal evolution of pressure P compared to P (F→S) at points S 4 (a) and F 1 (c). Relative difference 100 × (P -P (F→S))/P (F→S) at points S 4 (b) and F 1 (d).

Figure 18 (

 18 Figure 18 (b) displays the relative difference at point S 4 , in terms of the norm of the displacement associated with feedback effects, namely u su s (F→ S).Note that because the displacement amplitude at a given point of the interface is null before the arrival of the incident shock, the maximum value at point S 4 , ||u s (F→ S)|| max , has been used to compute the relative difference, which therefore takes the following form: 100 × ||u su s (F→ S)||/||u s (F→ S)|| max .

1 Figure 19 :

 119 Figure 19: Two-way coupling with an initial eigenmode perturbation of the structure: F↔S (aluminium). Temporal evolution of pressure P compared to P (S→F) at points S 4 (a) and F 1 (c). Relative difference 100 × (P -P (S→F))/P (S→F) at points S 4 (b) and F 1 (d).

- 3 S 4 Figure 20 :

 3420 Figure 20: Two-way coupling with an initial eigenmode perturbation of the structure: F↔S (aluminium). Ratio (P -P (S→F))/|P (S→F)-P (F→S)|max at point S 4 (note that the maximum value of |P (S→F)-P (F→S)| is used to compute this ratio, instead of |P (S→F)-P (F→S)|, in order to avoid divisions by zero).

  a) shows that ||u su s (S→ F)|| is higher than ||u s (F→S)|| and very close to ||u su s (S→ F)-u s (F→S)||: the displacements generated by the pressure variations associated with the eigenmode motion of the structure constitute the main contribution to ||u su s (S→ F)||. Although ||u su s (S→ F)|| increases during the simulation, Fig.21 (b) indicates that its remains very small compared to the eigenmode vibration ||u s (S→F)|| itself. This is the reason why the present fully-coupled and the S→F simulations lead to the same pressure field and consequently to the same pressure impulses and shock accelerations.

7 S 4

 74 ||u s -u s (S->F)|| S 4 ||u s -u s (S->F)-u s (F->S)|| S 4 ||u s (F->S)|| b u s (S->F)||/||u s (S->F)|| max x 100 ||u s -u s (S->F)-u s (F->S)||/||u s (S->F)|| max x 100

Figure 21 :

 21 Figure 21: Two-way coupling with an initial eigenmode perturbation of the structure: F↔S (aluminium). Temporal evolution of the displacement ||u su s (S→F)|| and ||u su s (S→F)-u s (F→S)|| compared to ||u s (F→S)|| (a). Relative differences 100 × ||u su s (S→F)||/||u s (S→ F)||max and 100×||u s -u s (S→F)-u s (F→S)||/||u s (S→ F)||max at point S 4 (b).

Table 1 :

 1 Material properties for aluminium and academic material mat2 with corresponding wave speeds.

	3	100

Table 3 :

 3 One-way coupling: S (aluminium)→F. Time of arrival of waves (TS or PS) t 0 and pressure impulses at probes F 1 and F 2 for the rigid and the "S→F" cases.

	4.2. Two-way Coupling: Fluid ↔ Structure

  Mesh I (t c ) (Pa.s) at point F 1 I (t c ) (Pa.s) at point F 2

	M 1	3.697	2.845
	M 2	3.576	2.768
	M 3	3.646	2.822
	M 4	3.630	2.825
	M 5	3.664	2.809
	M 6	3.664	2.801

Table B.5: One-way coupling: Structure → Fluid (aluminium). Pressure impulses at points F 1 and F 2 .

Acknowledgments

This work was granted by the ANR (Agence Nationale de la Recherche) -MAPIE project, ANR-13-MONU-002. Computations were performed using HPC resources from the computing center of CentraleSupelec and ENS Paris-Saclay.

Appendix A. Fitting Function Corresponding to the Prescribed Displacements along S 1 S 2

The fitting function corresponding to the normal displacements along S 1 S 2 prescribed to the structure reads as follows:

u s,f it (s, t).n s = -(a 0 + a 1 cos (sw) + b 1 sin (sw) + a 2 cos (2sw) + b 2 sin (2sw) + a 3 cos (3sw)+b 3 sin (3sw)+a 4 cos (4sw)+b 4 sin (4sw))×(H (s-H s /10)-H (s-9.0H s /10.0)) × sin (2πF s t) with s the curvilinear abscissa, H the heaviside function, a 0 = 0.1443, a Table B.4 summarizes the grid combinations that were tested to choose the meshes employed for the simulations presented in Sections 3 and 4. Table B.4: Mesh combinations used for the mesh refinement study. M f and M s : mesh of the fluid domain and mesh of the structure domain, respectively; L f x : size of the fluid domain in the x direction; N f x , N f y : number of cells of the fluid mesh in the x and y directions; N f x,cyl : number of cells of the fluid mesh in the x direction in the obstacle zone between S 1 and S 3 ; N s elt =N f x,cyl /2: number of elements for each side (S 1 S 2 and S 2 S 3 ) of the obstacles.

The fluid and the structure solvers are advanced with the same time step, ∆t = min(∆t s , ∆t f ). One may also note that Structure → Fluid simulations and two-way coupled simulations (Section 4) focus on the early stages of the interaction between the shocks and the obstacles and can be performed on a shorter fluid domain than Fluid → Structure simulations.

The choice of the six combinations presented in Table B.4 was motivated by some case-dependent criteria, as explained below. For example, M 5 and M 6 differ by the number of cells in the fluid domain outside of the obstacle zone in the x direction (N f x -N f x,cyl ): the main aim of this modification was to increase the accuracy of the pressure prediction for the probe F 2 which is located in this outside region.