Keywords: moyenne fréquence, coques, TVRC, chocs pyrotechnique, fusée Ariane, vibrations mid-frequency, shells, VTCR, pyrotechnic shocks, Ariane, vibrations

Titre : Simulation de la réponse en basse-et moyenne-fréquence de chocs par une approche en fréquence

Mots clés

Contents

Simulation of impact response using a frequency approach This demonstrate that the degradation effect highlighted in Sections 3.7.1 and 3.7.2 is solved. The problem is described in Section 4.1.2. . . . . . .

List of Figures

Simulation of impact response using a frequency approach 4.10 Performance analyzes over a frequency band [2000; 4000] Hz of the numerical example described in Section 4.1.2. VTCR and FEM are compared considering time and memory consumptions. DoFs required are reported for the sake of completeness, a direct comparison is unmeaningful due to theory differences. . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Changement of ki to address different wave speeds in orthotropic materials.100 4.12 three-step convergence process of VTCR with in-plane inertia of the numerical example described in Section 4.2 at 3700 Hz. This process is explained in general in Section 3. 5.1 HSS3+ structure before and after the pyrotechnic explosions, its cross section with sensor locations, and a particular of the explosion paths. . . . 5.2 HSS3+ structure before and after the pyrotechnic explosions. . . . . . . . 5.3 Figure 5.3a shows the rocket section during the pyrotechnic detachment.

Detonations ignite in A and B and end in C and D following the four explosion paths indicated by the red arrows. The explosions are supersonic. Figure 5.3b highlights their Mach cones. The PFL cylinder is detached from the rocket More details are provided in Section 5.2. . . . . . . . . . 5.4 Equivalent HSS3+ structure before and after the pyrotechnic explosions. . 11 Time-frequency plots of a sensor near the explosion and the satellite. . . . 5.12 Time and frequency domain responses of an accelerometer near the explosion on a small time window at the beginning of the signal. . . . . . . 5.13 Time and frequency domain responses of an accelerometer near the satellite on a small time window at the beginning of the signal. . . . . . . . .

Simulation of impact response using a frequency approach 5.14 Time and frequency domain responses of an accelerometer near the explosion on a small time window at the end of the signal. . . . . . . . . . . 5.15 Time and frequency domain responses of an accelerometer near the satellite on a small time window at the end of the signal. . . . . . . . . . . . . 5.16 Qualitative behavior of the input function in time and frequency domain determined studying the HSS3+ real data. . . . . . . . . . . . . . . . . . 5.17 Comparison of the frequency responses on different time windows of a sensor near the satellite to determine the damping coefficient η. . . . . . . 5.18 In-plane time-frequency plot comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. . . . . . . 5.19 Out-of-plane time-frequency plot comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. . . . 5.20 Frequency response comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. . . . . . . . . . . . 5.21 Frequency response comparisons in a semilog scale between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. . 5.22 Generic SDOF of the SRS illustrated in Section 5.4. . . . . . . . . . . . . 5.23 SRS comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. . . . . . . . . . . . . . . . . . . 5.24 SRS comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. A semilog scale is used to highlight magnitude differences. . 3.4 Ray number and type for initial guess VTCR and converged VTCR. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.3 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . 3.5 DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 3.4.3. . . . . . . . . 3.6 Quantities of interest of the example in Section 3.4.4. . . . . . . . . . . . 3.7 Ray number and type for initial guess VTCR and converged VTCR. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.4 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . 3.8 DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 3.4.4. . . . . . . . . 3.9 Quantities of interest of the example in Section 3.4.5. . . . . . . . . . . . 3.10 Ray number and type for initial guess VTCR and converged VTCR. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.4 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . 3.11 DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 3.4.5. . . . . . . . . 3.12 Quantities of interest of the example in Section 3.7.2. . . . . . . . . . . . 3.13 Ray number and type of the VTCR solution illustrated in Figure 3.31 of the problem exposed in Section 3.7.2. . . . . . . . . . . . . . . . . . . . 4.1 Ray number and type for initial guess VTCR and converged VTCR. converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.4 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . .

Simulation of impact response using a frequency approach 4.2 DoFs, time, and memory consumption comparisons between new and old VTCR and FEM of the numerical example described in Section 3.4.4 and reanalyzed in Section 4.1.1. . . . . . . . . . . . . . . . . . . . . . . . . 91 4.3 Quantities of interest of the supported cylinder section. . . . . . . . . . . 95 4.4 Ray number and type for initial guess VTCR and converged VTCR with in-plane inertia described in Section 3.1.1. converged VTCR is the result of the three-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Figure 4.6 where the error threshold is err V TCR ≤ 0.01. NS and SS stands for Normal and Shear Stresses respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.5 Orthotropic material properties and frequency examined of the numerical example described in Section 4.2. . . . . . . . . . . . . . . . . . . . . . 101 4.6 Ray number and type for initial guess VTCR and converged VTCR described in Section 3.4.2 at 3700 Hz. converged VTCR is the result of the three-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Figure 4.12 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.7 Quantities of interest of the example in Section 4.5.1. . . . . . . . . . . . 106 4.8 Ray number and type for initial guess VTCR and converged VTCR at [1000,2000,4000] Hz. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 4. Simulation of impact response using a frequency approach

Introduction

Aerospace and automotive industries are increasingly interested in numerical prediction techniques. Virtual testing reduces real prototype use and simplifies and speeds-up the design process. The automotive business is very competitive and making fast improvements is a key factor. For this reason, rapid and reliable numerical simulations are crucial. Differently, aerospace industries are interested in virtual testing since experimental campaigns are extremely costly or even impossible due to the uniqueness of every rocket or spacecraft.

Aircraft, spacecraft, rockets, and ground vehicles are all composed of complex shell structures. Shells are so widely used because of their high resistance and light weight. Lightness of a structure is a key characteristic of aircraft and spacecraft. Any unneeded weight increases fuel consumption and costs. Recently, this need has fostered a focused interest towards composites. These materials present extremely high stiffness-to-weight ratios because fibers are oriented along principal stress directions increasing resistance where is needed. For this reason, composites usually present orthotropic behavior. Equilibrium and constitutive equations for orthotropic shell structures are quite complex and in almost every real case no analytic solution is known. Hence, an approximated numerical prediction technique is needed.

Shock vibration problems still pose major difficulties to commercial codes since the relevant frequency range is wide, explosion loads are unknown, and structures are often composed of sandwich composite shells. The aim of the present thesis is to develop a reliable computer program capable of simulating these complex industrial problems to reduce the number of real tests required.

A vibration or acoustic response can be divided in three zones as shown in Figure 1: low-, mid-, and high-frequency.

The low-frequency range is almost unaffected by uncertainties and is characterized by the local response. The Finite Element Methods (FEM) are the most suited in this case. These approaches use a large but finite number of small elements to discretize the problem. Equilibrium equations on nodes are substituted with their respective finite difference set of equations. Once the nodal problem is solved, dynamic field variables are interpolated over the interior of elements using polynomial shape functions controlled by nodal values. A refined mesh is required to limit discretization error and pollution effect produced by the two approximations aforementioned.

The typical frequency response function in the high-frequency range is characterized by many small overlapping resonance peaks. The energy flows almost unrestrained

MF LF HF

Amplitude Frequency

Figure 1: A typical frequency response function divided in low-mid-and high-frequency zones.

among sub-domains and there are no dangerous resonance peaks. Moreover, the system is extremely sensible to uncertainties. For this reason, the local response is meaningless. The Statistical Energy Analysis (SEA) and its developments are the best suited methods in this frequency range. These techniques analyze the energy flow among big sub-domains neglecting the local response. In fact, these approaches determine only averages and variances of dynamic field variables over big sub-regions. These methods are based on some diffused field approximations that limit their applicability range to high-frequency. Moreover, their global nature and the lack of a local response investigation make them pointless at low-and mid-frequency.

Between the low-and the high-frequency lies the mid-frequency range. It shows hybrid characteristics between low-and high-frequency because it presents many high and partially overlapping resonance peaks. Hence, the local response is still required. Moreover, the system is more sensible to uncertainties than at low-frequency. For these and other reasons, high-and low-frequency methods are not fully suitable. On one hand, high-frequency methods does not provide the local response. On the other hand, computational costs become prohibitive for low-frequency techniques due to the high mesh refinement required to reduce pollution effect.

In the last decades, various approaches were proposed to study the mid-frequency range. On one hand, there are methods that extend FEM applicability range mitigating its drawbacks in various ways. The Adaptive Finite Element [START_REF] Stewart | h-Adaptive finite element computation of time-harmonic exterior acoustics problems in two dimensions[END_REF] refines the mesh where is needed. The Stabilized Finite Element Methods [START_REF] Babuška | A generalized finite element method for solving the helmholtz equation in two dimensions with minimal pollution[END_REF], Franca and Carmo 1989, Harari and Hughes 1992] add a residual error minimiza-Simulation of impact response using a frequency approach tion to reduce the pollution effect. The Domain Decomposition Approaches [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], Franca et al. 1997, Hughes 1995, Sandberg et al. 2001] reduce prohibitive computational costs parallelizing calculations. The Boundary Element Method [Banerjee 1993] solves the problem along boundaries extrapolating, afterwards, the solution in the whole domain. The Enrichment Approaches [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] increase FEM precision introducing shape functions that are more suited for the problem while decreasing the number of Degrees of Freedom (DoF). The meshless techniques [START_REF] Bouillard | Element-Free Galerkin solutions for Helmholtz problems: fomulation and numerical assessment of the pollution effect[END_REF], Lacroix et al. 2003, Suleau et al. 2000] discretize the domain in several points that have overlapping influence zones allowing the p field approximation.

On the other hand, there are techniques to extend the SEA to mid-frequency. The Energy Flow Analysis [START_REF] Ichchou | Energy models of one-dimensional, multi-propagative systems[END_REF], Lase et al. 1996] focuses on the energy density e. The Wave Intensity Analysis [Langley 1992, Langley et al. 1997] relaxes the diffuse field hypothesis. Nevertheless, the local response is not provided. The Statistical Modal Energy Distribution Analysis [Guyader et al. 1988, Totaro andGuyader 2012] adds a kinetic energy calculation in every sub-system. The Ray-Tracing Method [START_REF] Chae | Prediction of Vibrational Energy Distribution in the Thin Plate At High-Frequency Bands By Using the Ray Tracing Method[END_REF], Krokstad et al. 1968, Wilson and Hopkins 2015] approximates the vibrational field as a sum of rays (which are propagative plane waves) whose paths are followed until full damped. The Quasi-SEA method [Mace 2003] investigates the SEA problem using a modal energy approach.

There are hybrid methods that combine FEM and SEA. The Asymptotical Scaled Modal Analysis [De Rosa andFranco 2008, 2010] scales the system by means of SEA quantities and, after that, computes the vibrational response of the scaled system using a standard FEM calculation. The Hybrid Finite Element/Statistical Energy Analysis Method [START_REF] Cicirello | The vibro-acoustic analysis of built-up systems using a hybrid method with parametric and non-parametric uncertainties[END_REF], Genechten et al. 2011, Langley and Cordioli 2009, Langley and Cotoni 2007] decomposes the system in two different assemblies: one is considered at high-frequency and the other at low-frequency. The first one is treated by the SEA, the second one and the coupling conditions are addressed by the FEM.

Recently, some Trefftz methods were proposed to investigate mid-frequency. Shape functions are exact solutions of the equilibrium equations where some parameters are left unknown. The particular boundary conditions of the problem set these variables. Since these methods use Freedholm integrals to address boundary conditions, matrices are illconditioned. However, this problem does not affect the final solution as explained in [START_REF] Yeih | Numerical instability of the direct Trefftz method for Laplace problems for a 2D finite domain[END_REF]]. Typically, these approaches present a very fast convergence speed in terms of DoF compared to the FEM. The Ultra Weak Variational Formulation [START_REF] Cessenat | Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem[END_REF]] divides the domain in elements and interfaces associated with variables determined by continuity conditions among elements in weak form. T-elements discretize the domain in Trefftz elements and seek the solution imposing boundary conditions using either Lagrange multipliers [Freitas and Teixeira de Freitas 1999], or least square functional [START_REF] Jirousek | T-elements: state of the art and future trends[END_REF], Monk and Wang 1999, Stojek 1998]. The Discontinuous Enrichment Method [START_REF] Farhat | The discontinuous enrichment method[END_REF], Massimi et al. 2008a, 2010, Tezaur et al. 2014, 2008] introduces in a FEM shape function base set exact solutions of the equilibrium equations meeting boundary conditions using Lagrange multipliers. The Wave Boundary Element Method [START_REF] Perrey-Debain | Plane wave interpolation in direct collocation boundary element method for radiation and wave scattering: numerical aspects and applications[END_REF], 2004] is a BEM enhance-Simulation of impact response using a frequency approach ment where the polynomial boundary functions are multiplied by plane waves that are exact solution of the equilibrium equations. The Wave Based Method [START_REF] Atak | Coupling of Boundary Element and Wave Based Methods for the efficient solution of complex multiple scattering problems[END_REF], Deckers et al. 2009, Klanner and Ellermann 2015, Vanmaele et al. 2007, 2009] approximates the solution with plane waves and addresses boundary conditions by means of a variational residue weight.

The Variational Theory of Complex Rays (VTCR) [Ladevèze 1996] belongs to Trefftz methods. It was applied to: plate theory [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF], general shell theory [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF][START_REF] Riou | The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems[END_REF], transient dynamics [START_REF] Chevreuil | Transient analysis including the low-and the medium-frequency ranges of engineering structures[END_REF]], 3D acoustic [START_REF] Kovalevsky | The variational theory of complex rays for three-dimensional helmholtz problems[END_REF]], 2D unbounded domains [START_REF] Kovalevsky | On the use of the Variational Theory of Complex Rays for the analysis of 2-D exterior Helmholtz problem in an unbounded domain[END_REF], orthotropic plates [START_REF] Kovalevsky | A Trefftz approach for medium-frequency vibrations of orthotropic structures[END_REF], and on a wide frequency band [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF][START_REF] Ladevèze | [END_REF]. During the last decade three different versions were developed depending on the shape functions chosen: Fourier [Kovalevsky et al. 2012a], rays [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]], or sectors [START_REF] Ladevèze | Analysis of medium-frequency vibrations in a frequency range[END_REF].

For the explained advantages, this method was chosen to study shock propagation in the mid-frequency range. However, in order to study real complex industrial problems, many more improvements were implemented. The developed computer program is called Transient Analysis for PYROtechnic Shocks in Shells (TAPYROSS).

The experimental code is validated on a full-scale ground test (called HSS3+ test) conducted by the Centre National d' Études Spatiales (CNES) and Airbus Defence & Space (Airbus DS) to study the pyrotechnic detachment of the fairing of the European heavy-lift launch vehicle Ariane 5 (and 6 in the future). It is an expendable launch system used to deliver satellites into Geostationary Transfer Orbit (GTO) or Low Earth Orbit (LEO). The CNES and the European Space Agency (ESA) produce the Ariane rocket family. Airbus DS is the prime contractor for these launch vehicles. They are built in Europe and launched from the Guiana Space Center near Kourou in French Guiana. This location is particularly suitable since it fulfills the two crucial spaceport geographical requirements:

• it is close to the equator, so that the Earth spin provides some extra speed to the rocket if launched eastward;

• it has open sea to the east, so that boosters, lower stages, and rocket debris due to launch failures cannot fall on human habitations.

The rocket is composed (from the bottom to the top) by the Cryogenic Main Stage (CMS), the Storable Propellant Stage (SPS), two Solid Rocket Boosters (SRB), the payload, the SYstème de Lancement Double Ariane (SYLDA) structure 1 , and a fairing as illustrated in Figure 2. During launch, the CMS and the SRB bring the rocket in the upper atmosphere. The SPS is responsible for final satellite deployment. The payload is composed of two satellites, one inside and one on top of the SYLDA structure. The payload and all upper stages are covered at launch by a fairing which protects satellites and provide an aerodynamic shape to the rocket. It is discarded once sufficient altitude is reached (typically above 100 km). This operation is carried out by pyrotechnic charges that severe the protective nose cone from the rocket as depicted in Figure 3. Real tests illustrate that particularly intense vibrations occur. Since nose detachment happens in the upper atmosphere, acoustic vibrations are negligible due to the low air density. Conversely, structural explosion vibrations are very dangerous for the payload and must be assessed properly. Moreover, explosion loads are still unknown due to measurement difficulties.

Storable

Prop. The HSS3+ test is not an ad-hoc academic problem created to investigate performances of the chosen method. Conversely, it is a real industrial test where many small details complicate the problem. For this reason, it is the perfect testing ground for the developed computer program.

Stage

This dissertation follows the chronological order. Chapter 1 is devoted to the description of the HSS3+ test and the relevant literature analysis. Chapter 2 enunciates the general shell theory on which the various VTCR versions are developed. Chapter 3 illustrates the VTCR in the particular case of shallow shells introducing some theory and performance improvements. This approximation simplifies the weak form of the VTCR generalizing its formulation. At the end of the Chapter, some major drawbacks are unveiled. For this reason, Chapter 4 expands the VTCR and all its previously introduced enhancements to the general shell theory. In the same Chapter, many more theoretical and performance improvements are added to assess some complex features of the HSS3+ test.

The theoretical improvements introduced in Chapters 3 and 4 address:

• in-plane inertia (Section 3.1),

• particular solutions for general surface loads (Section 3.3),

• orthotropic and sandwich materials (Section 4.2),

• conic structures and shells of variable thickness (Section 4.3 ),

Simulation of impact response using a frequency approach

• moving loads (Section 4.4),

• concentrated masses and springs (Section 4.5),

Performance enhancements illustrated in Chapters 3 and 4 are:

• a quasi-symmetric ray distribution algorithm (Section 3.6),

• iterative solvers (Section 3.7),

• a fast integral computation technique (Section 4.6).

Chapter 5 illustrates the general shell-VTCR theory applied to the HSS3+ test comparing the real sensor data signals to the simulated acceleration functions. The last Chapter presents final remarks and conclusions.

Introduction

Simulation of impact response using a frequency approach Chapter 1

The industrial problem and its solution methods at mid-frequency

This Chapter presents the common background required to illustrate the Ph. D. work.A brief literature overview about mid-frequency methods is provided. The bibliography available is vast due to the relatively recent interest of the scientific community towards the mid-frequency problem. Four different strategies have been developed over the years:

• low-frequency methods improved to reduce computational costs required to reach the mid-frequency range

• high-frequency methods modified to address the local response which is required in mid-frequency

• hybrid methods that combine low-and mid-frequency method advantages

• dedicated mid-frequency methods

In the end, the VTCR is described in detail since it is the chosen method.

The literature survey

The bibliography about mid-frequency approaches is huge due to the interest of the scientific community over the years. For this reason, this Section proposes just a literature overview of the major mid-frequency methods. The study is not limited to mid-frequency methods in shells because, usually, mid-frequency approaches can be easily adapted to different PDE problems due to their general nature. The various approaches can be classified in:

• low-frequency methods improved to reduce computational costs required to reach the mid-frequency range

• high-frequency methods modified to address the local response which is required in mid-frequency

• hybrid methods that combine low-and mid-frequency method advantages

• dedicated mid-frequency methods

Low-frequency approaches

Low-frequency methods are well established and very effective at low-frequency. However, computational costs become prohibitive as the frequency rises. Moreover, the pollution effect [Deraemaeker 1999] further relegates these approaches to the low-frequency range. For these reasons, several methods were proposed to overcome such problems and effectively extend the applicability range of low-frequency methods to mid-frequency.

The following Sections illustrate such methods.

The Finite Element Method (FEM)

The present Section focuses on the FEM and its improvements. This is the standard approach to tackle low-frequency problems due to its intrinsic flexibility, its scalability, and its simple formulation. The FEM can address almost all PDE problems with very small adjustments and the error can be majored a priori. Moreover, it is easily implementable in a computer program.

The FEM [START_REF] Zienkiewicz | The Finite Element Method[END_REF] is based on a weak formulation of the equilibrium equations

b(u, v) = l(v) ∀v ∈ V = {v ∈ H 1 (Ω)|v = 0 along ∂Ω} (1.1)
where u = u(x) is the PDE vector field solution, x is the spatial variable, and v = v(x) is a generic smooth vector field equal to zero along boundaries (and on corners in case of shells). This formulation varies depending on the particular PDE problem treated but the Simulation of impact response using a frequency approach general philosophy remains the same. In shell theory it is formerly equal to the virtual work theorem discussed in Section 2.4. The domain Ω is divided into smaller sub-domains Ω e . The solution for every subdomain is approximated as a linear combination of polynomial shape functions

u(x) ≈ u h (x) = n el e=1 u e ψ e (x), x ∈ Ω e (1.2)
where ψ e are the shape functions and u e is the unknown coefficient vector of the particular sub-domain e. Once this discretized solution is introduced in the weak form, the problem is transformed in a linear set of equations. In matrix form:

Bu = l (1.3)
where B is the matrix relative to the bilinear form, l is the vector related to the linear form, and u = [u 1 , u 2 , ..., u n el ] is the vector of the unknown shape function coefficients.

The error [Deraemaeker 1999] can be majored by

ε EF ≤ α ε kh q q + β ε kL kh q 2q (1.4)
where k is the problem characteristic wavenumber, h is the maximum element size, q is the polynomial interpolation degree, L is a characteristic dimension, and α ε and β ε are some problem dependent constants. The first term of Equation (1.4) is related to the polynomial approximation while the second one concerns the pollution error. At lowfrequency the first term is dominant. For this reason, if kh is kept constant, the error remains the same as the frequency increases. Conversely, the other term keeps growing as the frequency increases. Hence, when the mid-frequency range is studied, the second term leads the error. Hence, mid-frequency is an inborn limit for the FEM. Several approaches were developed over the years to mitigate or bypass the pollution effect and to reduce computational costs. Such improvements are described in the following Paragraphs.

The adaptive FEM. Pollution error and computational costs rise when a very refined mesh is required to capture the small wave length at mid-frequency. In order to mitigate such problems, a non-uniform mesh refinement was proposed. In standard FEM the mesh size is related to the global frequency of the system. Yet, often the local frequency consistently varies among different sub-domains. For this reason, the adaptive FEM increases accuracy only where is needed. Different strategies were proposed over the years. They can be classified in three groups:

• h-refinement [START_REF] Stewart | h-Adaptive finite element computation of time-harmonic exterior acoustics problems in two dimensions[END_REF], Johnson 1990, Tie et al. 2003]. This strategy refines the mesh where high scattering zones are present.

• p-refinement [START_REF] Zienkiewicz | The finite element method for solid and structural mechanics[END_REF], Chaljub et al. 2007[START_REF] Vilotte | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF]]. This approach (sometimes called Spectral Element Method) introduces high order polynomial shape functions in the critical sub-domains.

Simulation of impact response using a frequency approach

• ph-refinement [START_REF] Demkowicz | Toward a universal adaptive finite element strategy, part 1. Constrained approximation and data structure[END_REF], Oden et al. 1989, Rachowicz et al. 1989].

It is a combination of the previous two.

All these methods are based on a priori and a posteriori local error estimators that indicate where additional accuracy is needed. Several error estimators were proposed ( [START_REF] Stewart | Explicit residual-based a posteriori error estimation for finite element discretizations of the Helmholtz equation: Computation of the constant and new measures of error estimator quality[END_REF], Stevenson 2007, Cascon et al. 2008]). Usually, they are optimal or quasi-optimal in an energetic sense.

Adaptive FEM is an excellent choice when the lower part of the mid-frequency range is studied. Conversely, it is not indicated if the mid or high part of the mid-frequency are investigated since pollution effect and computational costs are still present.

The stabilized FEM. The FEM pollution effect is related to the conditioning number of the quadratic form associated to the bilinear form of the weak formulation [Deraemaeker 1999]. During mesh refinement, when it becomes ill-conditioned, the pollution effect arises. For this reason, some approaches were proposed to stabilize the quadratic form. The Galerkin Least-Squares -FEM (GLS-FEM) stabilizes the quadratic form adding a term related to the minimization of the equilibrium residue [START_REF] Franca | Stabilized finite element methods[END_REF]]. [START_REF] Harari | Galerkin least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains[END_REF] illustrate that the technique completely erases the pollution error in 1D problems. However, successive works such as [START_REF] Thompson | A Galerkin leastsquares finite element method for the two-dimensional Helmholtz equation[END_REF] affirm that as the problem dimension increases, the pollution error is suppressed only along preferential directions. The Galerkin Gradient Least-Squares -FEM (G∇LS-FEM) is an evolution branch of the GLS-FEM [START_REF] Franca | The Galerkin gradient least-squares method[END_REF]. In this case the stabilizing term is related to the gradient of the equilibrium residue. [START_REF] Harari | Improved finite element methods for elastic waves[END_REF]] develop G∇LS-FEM for 3D, inhomogeneous, elastic wave propagation problems. The presented examples testify pollution error suppression in every direction for academic 2D problems. Yet, a very refined mesh is still required for mid-frequency. The Quasi Stabilized FEM [START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF] modifies matrices to stabilize the solution and to erase the pollution effect. 1D pollution error is suppressed and in 2D is minimized given a regular enough mesh. Yet, the implementation is quite complex.

Domain Decomposition Methods (DDM).

The pollution effect is one of the major FEM limitations. Stabilized FEM can mitigate or erase such problem very effectively. Nevertheless, very refined meshes are still required at mid-frequency. Such constraint leads to high computational costs. Recently, the scientific community focused on DDM to address such huge computational problems. The domain is split into smaller subdomains. Now, each smaller problem is affordable by a single computer. Therefore, subdomain problems are solved in parallel by several computers and, then, reconstructed. Split, sub-domain communication, and reconstruction are delicate procedures that needs to be highly optimized to ensure good performances. In fact, if communication overhead is not specifically addressed throughout the whole process, the algorithm could be outperformed by its non-parallelized counterpart.

Simulation of impact response using a frequency approach DDMs are very effective techniques to address mid-frequency problems. Yet, in the present work a different approach is preferred. DDM is a very interesting and complicate topic and a detailed analysis is out of the aims of this work. For the interested reader, Gosselet and Rey provide an excellent DDM survey in [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. The following Paragraphs illustrate the major DDM.

Component Mode Synthesis (CMS) developed in [MacNeal 1971] preliminary calculates eigenmodes for every sub-domain. After, global domain modes are investigated using local modes as a basis. A parallelization is performed to calculate sub-domain modes.

Guyan's decomposition introduced in [START_REF] Sandberg | Domain decomposition in acoustic and structure-acoustic analysis[END_REF] hierarchically classifies Degrees of Freedom (DoFs). Some DoFs are slaves and contribute to sub-domain mode computations. After, data are transmitted to the global system by a limited number of master nodes. These DoFs condense the information of their slave DoFs. Several calculators compute sub-domain problems in parallel.

Finite Element Tearing and Interconnecting (FETI) is one of the most popular DDM. It was first introduced in [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. Displacement problems of each subdomain are rearranged into a functional minimizations. Displacement continuity among sub-domains is ensured by Lagrange multipliers. Nodal variables (sub-domain by subdomain) condensation into Lagrange multipliers leads to a matrix system that determines the statically admissible internal loads. The final matrix problem is solved by an iterative algorithm to provide robustness. It was applied to acoustics in [START_REF] Magoulès | Application of a Domain Decomposition Method With Lagrange Multipliers To Acoustic Problems Arising From the Automotive Industry[END_REF], Farhat et al. 2000, Tezaur et al. 2001]. In [Mandel 2002] Mandel studies a coupled vibroacoustic problem. In [Farhat andMandel 1998, Farhat et al. 1998] the method is applied to plates and shells. In [START_REF] Klawonn | Dual-Primal FETI Methods for Three-Dimensional Elliptic Problems with Heterogeneous Coefficients[END_REF] the FETI is applied to 3D elliptic problems in general with heterogeneous coefficients.

The Multiscale FEM enriches the solution though a non-FE refinement. Several scales are added to the original FEM problem. Results are similar to the stabilized FEM since theories are developed from the same basis.

The variational multiscale FEM was first introduced by Hughes in [Hughes 1995]. The method approximates the solution at coarse scale as

u ≈ u p + u e (1.5)
where u is the complete solution, u p is the standard FEM solution, and u e is an enrichment that addresses the field behavior at a very refined scale. Enrichment functions are boundary zero-trace polynomials.

Simulation of impact response using a frequency approach

Residual Free Bubbles (RFB) lies on the same principles of the variational multiscale FEM [START_REF] Franca | Residual-free bubbles for the Helmholtz equation[END_REF]]. In this case, enrichment functions are residual-free bubbles that present zero-trace along boundaries.

The enrichment methods add non-polynomial shape functions that have a strong mechanical meaning to enrich the solution. In fact, classical polynomial shape functions used in standard FEM are not perfectly adapted for vibrational problems. In particular, such functions cannot address the very nature of the problem without approximations of some sort. These new functions are strictly related to the problem nature and need a priori qualitative knowledge of the solution. They can be added at convenience where needed. Such enrichment reduces DoFs required to reach convergence. In the same time, the matrix conditioning number grows at mid-and high-frequencies. In contrast to multiscale FEM where the enrichment is additive, in this case it is performed multiplying the standard solution with the enrichment functions. The Partition of Unity (PUM) [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF], Strouboulis and Hidajat 2006, Babuska and Melenk 1995] and the Generalized Finite Element Method (GFEM) [Strouboulis et al. 2000b,a, Fries and Belytschko 2010, Strouboulis et al. 2001] (which is a development of the PUM) belong to the enrichment methods.

The PUM is based on the following definitions and theorem here enunciated as in [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF].

Definition 1. Let Ω ⊂ R be an open set, {Ω i } be and open cover of Ω satisfying a pointwise overlap condition

∃M ∈ N, ∀x ∈ Ω, card{i|x ∈ Ω i } ≤ M (1.6)
Let {c i } be a Lipschitz partition of unity subordinate to cover {Ω i } satisfying

sup c i ⊂ closure (Ω) ∀i (1.7) i c i = 1 on Ω (1.8) c i L ∞ (R n e ) ≤ C ∞ (1.9) ∇c i L ∞ (R n e ) ≤ C G diam (Ω i ) (1.10) where C ∞ , C G are two constants. Then {c i } is called a (M,C ∞ ,C G ) partition
V := ∑ i c i V i = ∑ c i v i |v i ∈ V i ⊂ H 1 (Ω)
(1.11)

Simulation of impact response using a frequency approach is called the PUFEM space. The PUFEM space is said to be of degree m ∈ N if V ⊂ C m (Ω). The spaces V i are referred to as the local approximation spaces.

Theorem 1. Let Ω ⊂ R n e be given. Let {Ω i }, {c i }, and {V i } be as Definitions 1 and 2. Let u ∈ H 1 (Ω) be the function to be approximated. Assume that the local approximation spaces V i have the following properties: on each patch Ω i ∩ Ω, u can be approximated by

a function v i ∈ V i such that u -v i L 2 (Ω i ∩Ω) ≤ ε 1 (i) (1.12) ∇u -v i L 2 (Ω i ∩Ω) ≤ ε 2 (i) (1.13)
Then, the function

u ap = ∑ i c i v i ∈ V ⊂ H 1 (Ω) satisfies u -v i L 2 (Ω) ≤ C ∞ M ∑ i ε 2 1 (i) (1.14) ∇u -v i L 2 (Ω) ≤ 2M ∑ i C G diam(Ω i ) 2 ε 2 1 (i) +C 2 ∞ ε 2 2 (i) (1.15)
Since this property, it is possible to enrich the local approximation introducing v i functions having mechanical or topological contents.

The IsoGeometric Analysis (IGA) The IGA was first proposed in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] to combine FEM and the Computer Aided Design (CAD) programs. This approach is more geometrically based than the FEM since the polynomial shape functions are replaced with splines. The most commonly used type are the Non-Uniform Rational B-Splines (NURBS) defined recursively as:

N, 0 (ξ) = 1 if ξ i ≤ ξ < ξ i+1 0 otherwise (1.16) N, p (ξ) = N, p-1 (ξ) ξ -ξ i ξ i+p -ξ i + N i+1,p-1 (ξ) ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 p = 1, 2, 3, ..., (1.17)
where p is the NURBS degree and ξ i ∈ R are the elements of the generic knot vector ξ = ξ 1 , ξ 2 , ..., ξ n+p+1 . NURBS shares with meshless and enrichment methods some useful properties such as: the partition of unity, the non-iterpolatory nature of the basis, and the compact support. Refinements can are p-, h-, and ph-. Every refinement type displayed better performances than the standard FEM. It is illustrated that such improvements are largely due to the precise geometric description provided by spline curves. The approach was applied on structural vibrations in [Kolman 2012, Cottrell et al. 2006] and in particular on Kirchhoff-Love shells in [START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF]].

The Boundary Element Method (BEM)

The BEM is the standard for low-frequency problems in unbounded domains [Banerjee 1993, Bonnet 1999, Cruse 1969, Rizzo 1967]. First, the problem is solved along boundaries ∂Ω. Approximations are introduced at this point. In fact, boundaries are meshed and boundary constraints are imposed on boundary nodes. After, the solution is propagated inside the domain Ω thanks to the appropriate Green function g. The problem along ∂Ω is

c(x) • u(x) + ∂Ω u(y) [∇g(x, y)] n -g(x, y) [∇u(y)] ndy = 0 ∀x ∈ Ω (1.18)
where x is a spatial coordinate at which the solution is investigated, y is the spatial coordinate at which punctual forces are applied, u is the displacement solution, and c is a known function that addresses singularities. In contrast to FEM, matrices are fully populated. Hence, iterative solvers are employed to speed-up computations. Moreover, integral computations require more efforts. This is particularly true if x and y are close each other. For these reasons, even if there are much less DoFs then FEM, computational costs can be expensive or prohibitive for complex cases. In order to overcome this limitation, some improvements where proposed. They are described in the following Paragraphs.

The fast multipole BEM is the most famous BEM improvement [Bonnet 1999, Bonnet et al. 2008, Liu 2009, Nishimura 2002, Liu and Nishimura 2006, Shen and Liu 2007]. It imposes a decoupling hypothesis on x and y in Equation (1.18). This approximation holds if x and y are enough far away from each other. In this way, integrals can be computed separately and reused at convenience. Such approach greatly reduces memory and time consumption without significantly deteriorating the solution given that the approximation is applicable.

The clustering methods approximates matrix-vector multiplications extensively used by the iterative solver process . Such modification drastically reduces time and memory consumption. In [Hackbusch andNowak 1989, Sauter 1998], the proposed approach is the so called panel method. The O(q 2 ) arithmetical operations required in a standard matrixvector multiplication (B qq a q where the matrix is q × q) are reduced to O(q log d+2 (q)) where d is the dimension of the problem. The drawback is the reduced accuracy of the matrix multiplication. Usually, this problem becomes negligible in few steps since the solver is iterative. proposes an accurate survey of waves propagation in media.

The Isogeometric Analysis Boundary Element Method (IGABEM) is a migration of IGA into the general BEM theory [START_REF] Peake | Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems[END_REF], Simpson et al. 2014]. The BEM meshes just boundaries. IGA splines precisely describe boundary geometries. IGABEM merges these methods keeping the advantages of both. The synergy between BEM and IGA characteristics drastically increases BEM performances.

Simulation of impact response using a frequency approach

The meshless methods

This approach [START_REF] Idelsohn | The meshless finite element method[END_REF], Zhang et al. 2001] discretizes the domain in several points. Each one has an assigned influence region Ω e . Overlaps between influence regions are allowed. This is the main difference with the standard FEM. Such influence regions provide the support for the shape functions. The solution at any given point E is the average of the points that influence the point E. For this reason, this method is not interpolation based. It was applied to plates [START_REF] Ferreira | Natural frequencies of functionally graded plates by a meshless method[END_REF]]. The principal development is the Element-Free Galerkin Method.

The Element-Free Galerkin Method was developed in [START_REF] Bouillard | Element-Free Galerkin solutions for Helmholtz problems: fomulation and numerical assessment of the pollution effect[END_REF], Belytschko et al. 1994, Zhu and Atluri 1998] as an improvement of the Diffuse Element Method [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF]]. The unknown field u is approximated as a linear combination of ψ e (x) shape functions such that

u(x) ≈ u h (x) = n e ∑ E=1 u e ψ e (x) = u (x)ψ(x) ∀x ∈ Ω (1.19)
where ψ e (x) are shape functions constructed by a Moving Least Square procedure. Lagrange multipliers address Dirichlet boundary conditions. Works illustrate that the Element-Free Galerkin Method mitigate the pollution effect. This is particularly true when plane waves enrich the solution. It was applied to shells in [START_REF] Krysl | Analysis of thin shells by the Element-Free Galerkin method[END_REF].

High-frequency methods

The Statistical Energy Analysis (SEA), first introduced in [START_REF] Lyon | Power flow between linearly coupled oscillators[END_REF], dominates the high-frequency study. Its relatively simple theory, its inexpensive computational costs, and its effectiveness (when the basic assumptions are met) make this approach the strongly advised choice at high-frequency. SEA provides only spatial average values of the solution due to its statistical nature. The local behavior remains unknown. This is not a problem at high-frequency where holds the diffuse field hypothesis. Conversely, at mid-frequency such assumption is not fully respected. SEA spatial average energy previsions are still valid but are not sufficient since the response is still localized.

The following Sections present the standard SEA and its improvements to address midfrequency.

The SEA is commonly used as a research [START_REF] Lyon | Power flow between linearly coupled oscillators[END_REF] and an industrial tool [Troclet 1995]. This approach divides the domain in big sub-systems and, then, investigates the energy flow to extrapolate the average vibrational response. The SEA solution is the spatial average and variance [Langley andBrown 2004, Cotoni et al. 2005] of the energy for each sub-domain on a frequency band. Other interesting variables are

Simulation of impact response using a frequency approach determined a posteriori from the energy. The power balance of the sub-domain Ω i is

P in i = P diss i + ∑ j P coup i j
(1.20)

where P in i is the input power, P diss i is the dissipated power, and

P coup i j
is the exchanged power between Ω i and Ω j . Input and dissipated powers can be directly determined from problem data. The coupling term P, coup j is assumed to be proportional to the modal energy difference:

P coup i j = ωη i j E i n i - E j n j (1.21)
where n i and n j are the modal densities of Ω i and Ω j respectively, ω is the angular frequency, η i j is the coupling loss factor, and E i and E j are the energies of Ω i and Ω j respectively. The term ω considers the average over the frequency band. The SEA lies on some strong assumptions that are generally true only at very high frequency:

• the energy is transmitted only to adjacent sub-domains,

• the energy field is diffused in every sub-system, [Mace 2003] provides an excellent SEA review.

The Energy Flow Analysis is a SEA extension developed to address the local response. The problem unknowns are the energy densities e i . The energy flow I i is related to the energy density by:

I i = - c 2 gi ∇e i η i ω (1.22)
where c gi is the group velocity. Equation (1.20) becomes

P in i = ωη i e i - c 2 gi ∇e i η i ω (1.23)
This method has demonstrated its effectiveness for 1D problems [START_REF] Lase | Energy flow analysis of bars and beams: theoretical formulations[END_REF], Ichchou et al. 1997], but its application to more complex problems is difficult due to continuity conditions [START_REF] Langley | On the vibrational conductivity approach to high frequency dynamics for two-dimensional structural components[END_REF]]. Moreover the power balance provided in Equation (1.23) could be problematic in 2D. In fact, the 2D radiated field decays as 1/ √ r while the analytic theory affirms a 1/r decay factor [START_REF] Carcaterra | Thermal Analogy in Wave Energy Transfer: Theoretical and Experimental Analysis[END_REF].

Simulation of impact response using a frequency approach

The Wave Intensity Analysis is a natural SEA development where the diffuse field hypothesis is removed. In contrast to standard SEA, the field is considered directional [Langley 1992]. Propagative energy waves are associated with amplitudes. Assuming wave non-correlation, the energy E ji of the sub-system i of the propagative wave j can be expressed as a Fourier series

E ji (θ, ω) = ∑ p E jip (ω)N jip (θ)
(1.24)

The power balance between sub-systems provide the amplitudes E jip (ω). This approach presented considerably better results than SEA on plate assemblies [START_REF] Langley | Statistical energy analysis of periodically stiffened damped plate structures[END_REF]]. However, the local response is not addressed and the coupling coefficients could be hard to find.

The Statistical modal Energy distribution Analysis (SmEdA) is a very interesting SEA extension [START_REF] Guyader | Methode des coefficients d'influence energetiques[END_REF]. It relies on the SEA energetic approach. In addition, this technique adds a kinetic energy calculation for each sub-system to extrapolate the local behavior and address the mid-frequency. Energy transfer coefficients are evaluated from the coupled eigenmodes which are computed a priori using a FEM snapshot technique. Additionally, the method describes how the energy flow among sub-domains is influenced by geometry and damping. Many successful applications were presented over the years [Maxit and Guyader 2001b,a, Totaro et al. 2009, Maxit and Guyader 2003, Totaro and Guyader 2012]. In particular, plate assemblies treated with the SmEdA present very good performances. The principal drawback is the FEM snapshot computation requirement. In contrast to SEA simplicity, this part could be very expensive in terms of computational costs.

The Ray Tracing Method (RTM) also called Dynamical Energy Analysis (DEA) is a fusion of SEA with linear optic theory. It was first introduced in [START_REF] Krokstad | Calculating the acoustical room response by the use of a ray tracing technique[END_REF]] to predict acoustic performances in rooms. The vibrational response is calculated following a set of propagative waves until fully damped. Transmissions and reflections are computed using the classical Snell formula. If frequency and damping are enough elevate, the RTM is cheap and accurate. Otherwise, computational costs explode. Moreover, complex geometries are difficult to study due to their high scattering behavior. This technique was applied to acoustic [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF], Yang et al. 1998, Tanner 2009, Chappell et al. 2011], to plate assemblies [Chae andIh 2001, Chappell et al. 2014], and to 2D and 3D problems in general [START_REF] Chappell | Boundary element dynamical energy analysis: A versatile method for solving two or three dimensional wave problems in the high frequency limit[END_REF]].

The Quasi-SEA method introduces a modal approach to relax some SEA assumption and extends its applicability range to mid-frequency. In particular, indirect coupling among sub-domains are considered. Energy influence coefficients are expressed in terms of structure modes. Rain-on-the-roof excitation over the frequency band is assumed. The main drawback is the additional assumption on mode density. Broadly, the mode shapes Simulation of impact response using a frequency approach in the frequency band should be typical enough in terms of energy distribution throughout the system. Mace introduced and developed this method in [Mace 2003[Mace , 2005]].

Mid-frequency methods

Mid-frequency methods can be divided in hybrid and Trefftz approaches. [START_REF] Desmet | Methodologies for mid-frequency analysis in vibration and acoustics[END_REF]] presents a general overview of these theories.

Hybrid methods combine low-and high-frequency methods to address the mid-frequency range. The synergy reduce or bypass the limitations of both methods. In particular, SEA can be coupled with FEM, BEM, or a Trefftz method (described in Section 1.1.3). Some examples of hybrid approaches are provided in [Van Hal et al. 2005, Deckers et al. 2011, Cicirello and Langley 2013, Atak et al. 2014].

Trefftz methods approximate the solution as a linear combination of shape functions that satisfy a priori the equilibrium equations. The boundary conditions are just approximated. The Trefftz solution converges towards the exact solution under the complete metric space hypothesis [Herrera 1984]. Such assumption is called also T-completion. Usually, Trefftz approaches present fast convergence rates in terms of DoFs and illconditioned matrices. [Zienkiewicz 1997, Kita and Kamiya 1995, Pluymers et al. 2007, Liu et al. 2006] provide good reviews of Trefftz methods applied in different contexts. The following Sections describe the major Trefftz methods. The Variational Theory of Complex Rays (VTCR) is indeed a Trefftz method. Yet, its general version is introduced in a reserved separate Section Section 1.2 since it is the chosen method used in the rest of the work.

The hybrid Finite Element/Statistical Energy Analysis

This method was introduced in [Shorter and Langley 2005] and reviewed in [START_REF] Cicirello | The hybrid finite element/statistical energy analysis method[END_REF]. The method decomposes a system in two different assemblies, called the master and the slave systems, with common characteristics. The master system presents a deterministic response. Such assembly can be easily treated by a standard FEM. Conversely, the slave system shows a randomized response. This model is investigated by a SEA approach. The system coupling is performed throughout the master system. In this way the method inherits SEA advantages. In fact, the uncertainties are directly addressed by the SEA part which directly provides average and variance without any information on stochastic parameters. Moreover, this approach does not require any Montecarlo simulation and its very flexible due to FEM study. Parametric uncertainties are taken into account in [START_REF] Cicirello | The vibro-acoustic analysis of built-up systems using a hybrid method with parametric and non-parametric uncertainties[END_REF]. Other applications of the method are illustrated in [START_REF] Langley | Response variance prediction for uncertain vibro-acoustic systems using a hybrid deterministic-statistical method[END_REF], Low and Langley 2008, Langley and Cordioli 2009, Langley and Bremner 1999, Cotoni et al. 2007].

Asymptotical Scaled Modal Analysis

The Asymptotical Scaled Modal Analysis was first introduced for plates in [De Rosa andFranco 2008, 2010]. The proposed approach scales the system into a smaller model that Simulation of impact response using a frequency approach can be addressed by standard FEM. An equivalence based on the quadratic velocity links the original and the scaled system. The quadratic velocity is a global system characteristic that can be computed by an inexpensive SEA analysis. In particular, damping and geometric dimensions are scaled so that the quadratic velocity remains unchanged between the original and the scaled models. For this reason, the scaled system retains all the energetic properties of the original model. This technique is very inexpensive in terms of computational costs and can be directly coupled with a standard FEM. It was applied to metallic and composite plate assemblies.

The Ultra Weak Variational Formulation (UWVF)

The UWVF discretizes the domain in elements adding a variable for each interface. A weak formulation of the boundary conditions sets such unknowns. For each element the vibrational field is approximated as a linear combination of plane waves where a Galerkin approach derived from the weak form provides the boundary variables. Once interface variables are computed, the solution field inside each element is locally extrapolated. Even if the UWVF presents inexpensive computational costs, it is an h-method and matrices are ill-conditioned. For this reason, when large problems are investigated, a preconditioner is introduced. Cessenat and Despres in [START_REF] Cessenat | Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem[END_REF]] demonstrate that a uniform distribution of the propagation directions maximizes the UWVF matrix determinant. The UWVF is compared to the PUM1 in [START_REF] Huttunen | Comparison of two wave element methods for the Helmholtz problem[END_REF]] on 2D Helmholtz problems with irregular meshes. The authors conclude that both methods perform well with coarse meshes. However, the UWVF seems more effective at mid-frequency while the PUM seems to perform better at low-frequency . The UWVF presents worse conditioning numbers than PUM even after the pre-conditioning process. In [START_REF] Gittelson | Plane wave discontinuous Galerkin methods: analysis of the h-version[END_REF]] the authors explain that the UVWF is a special case of the Discontinuous Galerkin approaches2 using plane waves.

T-elements

The methods that use the so called T-elements divide the system in sub-domains and investigate the local problem using a Trefftz approach. Therefore, these sub-domains are called T-elements. They differ from each other on how are addressed the coupling conditions. In particular, continuity conditions are taken into account with either Lagrange multiplier [Freitas and Teixeira de Freitas 1999], or least square functional [START_REF] Jirousek | T-elements: state of the art and future trends[END_REF], Monk and Wang 1999, Stojek 1998]. In the last case T-elements are called Least-Squares T-elements (LST). Since only boundary conditions should be verified, these methods produce small matrices compared to FEM reducing computational costs. Moreover, complex geometries and solution singularities can be addressed using sufficiently small T-elements and injecting ad-hoc functions in the approximation space. [Zielinski 1997] improves the conditioning number of the matrix produced by T-elements and [START_REF] Freitas | Three-dimensional hybrid-Trefftz stress elements[END_REF] present 3D structure T-elements.

1.1.3.5 The Discontinuous Enrichment Method (DEM)

The DEM was first introduced in [START_REF] Farhat | The discontinuous enrichment method[END_REF]. It is similar to multi-scale FEM but, in contrast to multi-scale FEM which uses zero-trace functions over boundaries, the DEM enriches the shape function basis with Trefftz functions. Such functions identically satisfy equilibrium equations. Boundary and coupling conditions are met using Lagrange multipliers. In particular, these functions are equally-distributed plane waves. The infsup condition enunciated in [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] relates the Lagrange multipliers to the plane waves of the element. Hence, special elements are introduced. 2D elements were proposed in [Farhat et al. 2003a[Farhat et al. ,b, 2004a,b] ,b] while 3D elements were presented in [START_REF] Tezaur | Three-dimensional discontinuous Galerkin elements with plane waves and Lagrange multipliers for the solution of midfrequency Helmholtz problems[END_REF]]. [Farhat et al. 2004a] proved that the multi-scale FE approximation (polynomial function enrichment) does not improve accuracy when many Henholtz problems are considered. The DEM provides accurate results with coarser meshes than FEM and matrices are better conditioned than PUM [START_REF] Grosu | Studies of the discontinuous enrichment method for two-dimensional acoustics[END_REF]. Recently, the DEM was applied to acoustics [Gabard 2007], to fluid/fluid and fluid/solid interaction problems [START_REF] Massimi | A discontinuous enrichment method for three-dimensional multiscale harmonic wave propagation problems in multi-fluid and fluid-solid media[END_REF], to plate assemblies [START_REF] Massimi | A discontinuous enrichment method for the efficient solution of plate vibration problems in the medium frequency regime[END_REF], Zhang et al. 2006], and to high Péclet advection-diffusion problems [START_REF] Kalashnikova | A discontinuous enrichment method for the finite element solution of high Péclet advection-diffusion problems[END_REF]]. Moreover, [START_REF] Wang | Overview of the discontinuous enrichment method, the ultra-weak variational formulation, and the partition of unity method for acoustic scattering in the medium frequency regime and performance comparisons[END_REF] compares DEM with UWVF and PUM methods. In conclusion, this method provides accurate results with less DoFs than FEM and the pollution effect is drastically mitigated by the plane waves. Moreover, complex geometries can be addressed due to the FE mesh. Yet, the Lagrange multiplier involvement hinders flexibility.

The Wave Boundary Element Method (WBEM)

The WBEM is a direct BEM extension. It was first presented in [START_REF] Perrey-Debain | Plane wave interpolation in direct collocation boundary element method for radiation and wave scattering: numerical aspects and applications[END_REF], 2004]. In contrast to low-frequency BEM extensions, this approach introduces plane waves to enrich the solution. In particular, polynomial shape functions are multiplied by these Trefftz functions. Such functions permit coarser meshes along boundaries. There are no a priori restrictions on Trefftz shape function number and distribution algorithm. For this reason, the method is very flexible. Moreover, the WBEM reaches good accuracy levels with less DoFs than the standard BEM. Yet, the drawback is the degradation of the matrix conditioning number due to plane waves.

The Wave Based Method (WBM)

The WBM was first introduced in [START_REF] Desmet | A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems[END_REF]. The method approximates the solution as a linear combination of plane waves. In acoustics, the solution in the sub-Simulation of impact response using a frequency approach domain Ω i is expressed as

p i (x, y) = ∞ ∑ m=0 a mi cos mπx L xi e ±ı k 2 -mπx L xi y + ∞ ∑ n=0 a ni cos mπy L yi e ±ı k 2 -mπy L yi x (1.25)
where L xi and L yi are the dimensions of the smallest encompassing rectangle of the subcavity Ω i . The series are truncated to numerically implement the approach. The number of shape functions is chosen accordingly with the following criteria

n xi L xi ≈ n yi L yi ≈ T k π (1.26)
where T is the chosen truncation parameter. Desmet in [Desmet 1998] suggests to impose at least T = 2. The weighted residual variational technique addresses boundary and coupling conditions. In particular, each condition along boundaries is multiplied by its proper dual quantity and integrated along the boundary itself. In acoustic the pressure dual variable is the normal fluid speed through the boundary( and vice versa). This Galerkin weak form leads to a matrix system. Its solution is a vector of complex wave amplitudes. The solution uniqueness property is ensured only if damping is present. Numerical tests affirm that if boundary and propagation media impedances are close each other, the formulation is more stable [START_REF] Pluymers | Trefftz-based methods for time-harmonic acoustics[END_REF]]. Typically, p-refinements perform better than h-refinement. It is very similar to the Variational Theory of Complex Rays illustrated in detail in Section 1.2. The differences lies in the used shape functions and in the way boundary conditions are addressed. The present approach was applied to several different problem types. In particular, it was applied to 2D and 3D acoustics [START_REF] Desmet | A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems[END_REF], to plate assemblies [START_REF] Vanmaele | An efficient wave based prediction technique for plate bending vibrations[END_REF][START_REF] Vanmaele | An efficient wave based prediction technique for dynamic plate bending problems with corner stress singularities[END_REF], to thick shells [START_REF] Klanner | Wave Based Method for the steady-state vibrations of thick plates[END_REF], to poroelastic materials [START_REF] Deckers | A novel wave based method for the efficient 2D dynamic modelling of the poro-elastic Biot equations[END_REF][START_REF] Deckers | Efficient treatment of stress singularities in poroelastic wave based models using special purpose enrichment functions[END_REF][START_REF] Deckers | A Wave Based Method for the efficient solution of the 2D poroelastic Biot equations[END_REF], to structural-acoustic radiation analysis [START_REF] Pluymers | On the use of a wave based prediction technique for steady-state structural-acoustic radiation analysis[END_REF]], 3D Helmholtz problems in general [ Van Genechten et al. 2012], and to unbounded domains [START_REF] Genechten | A Trefftz-based numerical modelling framework for Helmholtz problems with complex multiple-scatterer configurations[END_REF]].

The Variational Theory of Complex Rays (VTCR)

This Section completes the literature survey introducing the state-of-the-art of the VTCR before the present thesis. The reader interested to the latest version of the VTCR theory applied to thin shells can see Chapter 4.

The VTCR was first introduced in [Ladevèze 1996]. Since this approach is very versatile, over the years it was applied to:

• plate theory [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF]],

• general shell theory [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF][START_REF] Riou | The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems[END_REF],

Simulation of impact response using a frequency approach • transient dynamics [START_REF] Chevreuil | Transient analysis including the low-and the medium-frequency ranges of engineering structures[END_REF]],

• 3D acoustic [START_REF] Kovalevsky | The variational theory of complex rays for three-dimensional helmholtz problems[END_REF]],

• 2D unbounded domains [START_REF] Kovalevsky | On the use of the Variational Theory of Complex Rays for the analysis of 2-D exterior Helmholtz problem in an unbounded domain[END_REF]],

• orthotropic plates [START_REF] Kovalevsky | A Trefftz approach for medium-frequency vibrations of orthotropic structures[END_REF]],

• pyrotechnic shocks [Bézier 2014],

• a wide frequency band [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF], Ladevèze and Riou 2005, Ladevèze et al. 2012].

The VTCR is a Trefftz method specifically developed to address mid-frequency problems. The system is divided in big sub-domains and the introduced shape functions identically verify the equilibrium equations. A weak formulation addresses boundary, corner, and coupling conditions. The following Sections describe the shape function types and the weak form. In particular, these Sections report the theory exposed in [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF], Riou et al. 2004, Kovalevsky et al. 2012a, 2013] and properly reworked to be consistent with the rest of the thesis.

The shape functions

Be u the solution of a linear partial differential equation problem. By definition, u verifies equilibrium equations and boundary conditions. u can be expressed without approximations in 2D and 3D as

u(θ) 2D = 2π 0 a(θ) 2D e ık(θ) r rel dθ + u p (1.27) u(θ, ξ) 3D = π 0 2π 0 a(θ, ξ) 3D e ık(θ,ξ) r rel dθ sin(ξ)dξ + u p (1.28)
where r rel is a generic spatial vector variable relative to the sub-domain, a(θ) (or a(θ, ξ)) is the unknown amplitude function, u p is a particular solution, and k(θ

) (or k(θ, ξ)) is the wave direction vector. It is further decomposed in k(θ) 2D = k k(θ) 2D k(θ, xi) 3D = k k(θ, ξ) 3D (1.29)
where k is the wavenumber determined by the dispersion equation and k(θ) 2D (or k(θ, ξ) 3D in 3D) is a unit direction vector. In acoustics, only propagative waves are required and the dispersion equation is k = ω/c where c is the speed of sound [Kovalevsky et al. 2012a].

Conversely, in thin shell theory this is not sufficient [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF]]. In fact, the shape function set should be enriched with semi-evanescent waves and the dispersion equation is extrapolated from the equilibrium equations neglecting in-plane inertia. For

Simulation of impact response using a frequency approach this reason, in-plane waves are disregarded. Propagative and semi-evanescent shape functions differ on the formula for the unit direction vector k. Section 3.1 describes in detail the various shape function types required in shell theory. The amplitude function a(θ) (or a(θ, ξ) in 3D) is discretized. This is the only VTCR approximation. Over the years, three different VTCR versions were proposed depending on the particular discretization method used. They are qualitatively illustrated in Figure 1.1. Each algorithm approximates the amplitude function with a different converging series. For the sake of rigorousness, each illustrated series approaches infinity. However, they are truncated due to computational implementation reasons. The Ray-VTCR was introduced in solid mechanics in [Ladevèze 1996]. Plane waves were proposed as shape functions. In this case, the generic solution u becomes

a(θ) 2D ≈ ∞ ∑ i=1 a i δ(θ -θ i ) (1.30) a(θ) 3D ≈ ∞ ∑ j=1 ∞ ∑ i=1 a i j δ(θ -θ i , ξ -ξ j ) (1.31)
where θ i (and ξ j in the 3D case) is the sampled version of the continuous function θ and a i (or a i j ) is an element of an unknown amplitude vector determined by a weak form (see Section 1.2.2 for an introduction and Section 3.1 for a detailed description). This is the oldest and most versatile VTCR version. These functions can address 2D and 3D scenarios, bounded and unbounded systems, and shells and acoustics problems with simple adjustments. In particular, Dirac deltas erase the integrals in Equations (1.27) and (1.28) providing great flexibility. For this reason, the integrals of the weak form can be computed analytically if some conditions are verified. Typically, the computational costs required to create the matrices are comparable or higher than computational costs of the final matrix inversion due to the small number of DoFs and the high number of required Simulation of impact response using a frequency approach numerical integrals. Hence, this characteristic can dramatically reduce computational costs.

The sector-VTCR was presented in [START_REF] Ladevèze | Analysis of medium-frequency vibrations in a frequency range[END_REF]] in plates and in [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]] for thin shells. This algorithm approximates the amplitude function with constant angular sectors

a(θ) 2D ≈ ∞ ∑ i=1 a i h i (θ) (1.32)
where h i (θ) are rectangular functions

h i (θ) =    0 if θ < θ i 1 if θ i ≤ θ < θ i+1 0 if θ i+1 ≤ θ (1.33)
In contrast to ray-VTCR, these functions do not simplify integrals and their extension to 3D scenarios is difficult. In fact, there is no literature about 3D problems tackled with the sector-VTCR. Moreover, if the amplitude is assumed continue, the mean value theorem applied to Equation (1.27) affirms that the integral over a sector is equivalent to a particular plane wave whose direction lies in the angular sector. Hence, in a real test where many DoFs are required regardless to the adopted discretization method, the ray-VTCR outperforms the sector-VTCR due to the numerically evaluated integrals required by the latter. For this reason, this approximation algorithm has been abandoned in the last few years.

The Fourier-VTCR was first introduced in 2D [ Kovalevsky et al. 2012a] and 3D acoustics. The amplitude is a scalar function and it is approximated with a Fourier series

a(θ) 2D ≈ ∞ ∑ i=-∞ a i e ıiθ (1.34) a(θ, ξ) 3D ≈ ∞ ∑ l=0 l ∑ m=-l a lm y m l (θ, ξ) (1.35) (1.36)
where y m l (θ, ξ) denotes the non-negative spherical harmonics of order l and index m, given by

y m l (θ, ξ) = 2l + 1 4π (l -m)! (l + m)! p m l (cos(ξ)) e ımθ (1.37) with p m l ( ) = (-1) m (1 -2 ) m 2 ∂ m p l ( ) ∂ m with -l ≤ m ≤ l (1.38)
Simulation of impact response using a frequency approach where p l ( ) are the Legendre polynomial of the generic function and p m l ( ) are the associated Legendre functions. Even if this approach does not bypass the integrals of Equation (1.27) (or Equation (1.28)), it offers some key properties. The double (or triple in the 3D scenario) integrals in the weak form3 can be simplified to simple integrals. Moreover, an energetic preconditioner can be constructed to reduce ill-conditioning problems. However, integrals of the weak form must be computed numerically and the extension to shell theory is difficult due to the complex formulation of the semi-evanescent waves.

The weak formulation

A weak form addresses boundary, corner, and coupling conditions. The focus is on the generic boundary condition

u = u along ∂Ω (1.39)
where u is the solution, u is the boundary condition, and ∂Ω is the boundary. Since the solution is defined in a Banach space, an integral norm is defined. In formula

, w f = ( )ds (1.40)
where and are two generic functions in the Banach space. Calculus of variations affirms that the boundary condition in Equation (1.39) is equivalent to u e , v e w f = ue, v e w f ∀v e

(1.41)

where v e is a generic function defined on a sub-domain Ω e that has the boundary ∂Ω e .

Broadly, instead of directly investigating the boundary condition, the weak form tests it studying the norms produced by the solution and the boundary condition with test functions. If this infinite set of equations is met, then the boundary condition is satisfied in strong form too. It is in general for all the boundary, corner, and coupling conditions

b(u, v) = l(u, v) ∀v (1.42)
where b(u, v) is the bilinear form and l(u, v) is the linear form composed by the left and the right side of Equation (1.41) respectively. These forms are related to an energetic quantity that has a strong physical meaning. In particular, in shells they are obtained from the boundary version of the virtual work theorem. In Section 3.1 they are further specialized for the shell case reported in Equation (2.58). In order to obtain a single and unique solution set, the number of test functions is equal to the number of variables. Moreover, shape functions themselves are used as test function being the VTCR is a Galerkin method. Injecting the shape functions in Equation (1.42) leads to a set of linear equations. In matrix form

Ba = l (1.43)
where B and l are a matrix and a vector related to the bilinear and the linear forms respectively and a is the amplitude vector. B presents Fredhold integrals of the first kind. [Sourcis 2009, Liu et al. 2006] demonstrate that these integrals generate ill-conditioning problems. This characteristic destabilizes the solution of the matrix inversion which is the amplitude vector. However, the final solution remains stable because of the final multiplication of the amplitude vector by the shape functions. Such behavior is typical of the Trefftz methods and can be detrimental if the method is coupled with a Reduced Order Model (ROM) to extend the study to a frequency band [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF]. For this reason, preconditioners and iterative solver are suggested.

Simulation of impact response using a frequency approach Chapter 2

The general shell theory

The Koiter orthotropic shell theory is presented to mathematically formulate the vibrational problem. It is akin to the one provided in [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF]] and [ Van der Heijden et al. 2008]. At the same time, the notations are introduced. Some energetic quantities are defined and an interesting virtual work theorem version is defined. They will be useful in the following Chapters. Finally, the shallow shell approximations required in Chapter 3 are introduced.

Equilibrium equations

The general reference example is presented in Figure 2.1. The focus is on a generic subdomain Ω e of the frame structure in Figure 2.1. The number of sub-domains is n e . Various boundary, corner, coupling, and surface conditions are split in Figure 2.2 to better identify the various constraints. The symbol refers to a generic quantity. From now on, every quantity is referred to Ω e and the subscript e is omitted for simplicity. Nevertheless, it occasionally reappears when is required. The ∂ Ω symbol refers to a generic boundary where a condition is applied. In the particular case of a boundary shared among subdomains, Γ is used instead. In the same way, for conditions applied on corners, a symbol

∂∂ Ω is used. The generic corner shared among sub-domains is indicated with C . The over-line symbol indicates that a quantity is known (i.e. a value of a boundary constraint). n indicates an outward normal unit vector of a boundary. Ω is subject to loads, Simulation of impact response using a frequency approach displacements constraints, and continuity conditions along boundaries (Figure 2.2a) and on corners (Figure 2.2b) as well as a distributed load per unit surface g (Figure 2.2d). Without loss of generality, a displacement constraint u = [v , w] along ∂ u Ω can be divided in in-plane v and out-of-plane w components 1 . In the same way, a load per unit length p = b , q along ∂ p Ω can be divided in in-plane b and out-of-plane q components. A rotation condition w, n is imposed along ∂ w, n Ω while a bending moment per unit length m is applied along ∂ m Ω. Corners are subject to out-of-plane displacements constraints w C on ∂∂ w C Ω and point forces q C on ∂∂ q C Ω. Coupling conditions are applied on C and along Γ, in order to ensure stresses and displacements continuity among sub-domains (Figure 2.2c). Since the VTCR is a frequency approach, the present discussion focuses on equilibrium equations in the frequency-domain. In order to pass from the time-to the frequencydomain and vice versa, the Fourier transform and anti-transform are used

ω = +∞ -∞ t e -ıωt dt (2.1) t = 1 2π +∞ -∞ ω e ıωt dω (2.2)
where ω = 2π f is the angular frequency, f is the frequency, the subscripts ω and t refer to the frequency-and the time-domain respectively, and ı = √ -1 is the imaginary unit. The equilibrium equations in the time-domain are formally equal to their corresponding frequency-domain versions except for the time derivatives. In fact,

t , t = ıω ω (2.3)
Since the discussion is on the frequency domain, the subscript ω is hereafter omitted.

1 is the transpose operator.

Simulation of impact response using a frequency approach Curvilinear coordinates α, β, ẑ define the shell geometry. α, β are unit vectors orthogonal each other and tangent to shell middle surface in every point. ẑ is the out-ofplane coordinate and it is determined in each surface point by ẑ = α × β where × is the cross product. The generic surface is described by the vector field

r = r(α, β) (2.4)
where r : R 2 → R 3 is a generic 3D position vector. Lamé surface parameters L α , L β and curvature radii R α , R β are

L α = r , α r, α (2.5) L β = r , β r, β (2.6) R α = L 2 α ẑ r, αα (2.7) R β = L 2 β ẑ r, ββ (2.8)
where commas indicate directional derivatives. These definitions will be useful in the following Sections. Since the general shell theory is too complex to be directly treated by the VTCR, hereafter we restrict the study to surfaces with constant Lamé parameters and curvature radii. Displacement field can be restricted to (Koiter's kinematics assumptions)

Simulation of impact response using a frequency approach

u z = u -zφ (2.9) φ = ∇w -R • v (2.10) R = 1 R α 0 0 1 R β (2.11)
where u z is the displacement thorough shell thickness, u, v and w are respectively total, inplane and out-of-plane displacements of the middle surface and R is the curvature matrix.

The solution is researched in D = {u, N, M}. u is a finite-energy displacement and N and M are finite-energy generalized stress tensors. D i satisfies the in-plane equilibrium equation

∇ • N -R (∇ • M) + g αβ + ρhω 2 v = 0 (2.12)
and the out-of-plane equilibrium equation

∇ • (∇ • M) + R : N + g z + ρhω 2 w = 0 (2.13)
where distributed loads g are divided in in-plane and out-of-plane components in this way

g = g α , g β , g z (2.14) g αβ = g α , g β (2.15)
Matrices N and M are stresses and stress moment resultants along thickness respectively. Constitutive relations, also called stress-strain relations, are, for orthotropic materials, Simulation of impact response using a frequency approach

M = -D : (∇∇w + RE) (2.16) N = B : (E -Rw)
(2.17)

E = [∇v] sym = 1 2 ∇v + ∇v (2.18) B =   B α ν αβ B β 0 ν βα B α B β 0 0 0 B G   (2.19) B α = E α h 1 -ν αβ ν βα
(2.20)

B β = E β h 1 -ν αβ ν βα (2.21) B G = Gh (2.22) D =   D α ν αβ D β 0 ν βα D α D β 0 0 0 D G   (2.23) D α = E α h 3 12(1 -ν αβ ν βα )
(2.24)

D β = E β h 3 12(1 -ν αβ ν βα )
(2.25)

D G = Gh 6 (2.26) E α = E α0 (1 + ıη α ) (2.27) E β = E β0 1 + ıη β (2.28) G = G 0 (1 + ıη G ) (2.29)
where B and D are Hooke's plane stress operators concerning in-plane and out-of-plane stresses respectively, ρ is the density, h is the shell thickness, E α0 and E β0 are the Young moduli along directions α and β respectively, η α and η β are the relative damping coefficients of the Young moduli, ν αβ and ν βα are the Poisson's ratios (α, β) and (β, α) respectively, G αβ is the in-plane shear modulus, η G is its specific damping coefficient, : is the inner matrix product operator, [ ] sym = 1 2 ( + ) is the symmetric part operator, and N and M are the stress and stress moment resultants tensors respectively. The sub-space of D associated with homogenized conditions (g = 0) is denoted as D 0 = {u 0 , N 0 , M 0 }. This definition will be useful in Chapter 4.

Simulation of impact response using a frequency approach

The boundary conditions

For the sake of clarity, subscripts e and e that refer to Ω e and another sub-domain Ω e respectively are not omitted in this Section. In order to present a well-posed problem, three conditions must be imposed along each boundary and one on each corner. Boundary and corner conditions presented in Figure 2.1 can be classified in this way:

1. an in-plane condition, either a displacement constraint or a load per unit length (v e or b e ),

2. an out-of-plane condition, either a displacement constraint or a load per unit length (w e or q e ), 3. either a rotation or a bending moment per unit length (w e , ne or m e ),

4. an out-of-plane condition on corners, either a displacement constraint or a point load (w C e or q C e ).

Boundary conditions, in the most general case are

v e = v e
b e = [b αe , b βe ] = (N e -R e M e ) ne = b e - n Γe ∑ e Γe =1
b e Γe (2.33)

q e = (∇ • M e ) ne + ∇ t e M e ne • te = (∇ • M e ) ne + t e M e ne , te = q e - n Γe ∑ e Γe =1
q e Γe (2.34) w Simulation of impact response using a frequency approach

m e =
C e = w C c ∃w C c w C e w C c (2.36) q C e = t e1 M e ne1 + t e2 M e ne2 = q C c - n C e ∑ e C e =1

Energies

In shell theory and in case of forced vibrations some useful energy quantities can be defined. These are the strain energy E S , the kinetic energy E K , and the dissipation energy

E D E S = 1 2 n e ∑ e=1 Ω
E D =ℑ {E S } . (2.40)
where H is the hermitian operator. The operator ℑ { } extracts the imaginary part of . These will be extensively used.

Virtual work theorem

This Section enunciates the virtual work theorem for shells subject to forced vibrations. After some rearrangements, a useful theorem version is illustrated.

Theorem 2 (The virtual work theorem). For deformable solids the following equation holds

L in = L ex (2.41)
where L in and L ex are internal and external virtual works.

For the sake of brevity, its proof is skipped. In case of forced vibrations it can be developed by means of the previously defined energies

2 (E S -E K ) = L ex (2.42)
The quantity E S -E K can be rearranged through residues along boundaries and on corners. This equivalence lies on a divergence vector calculus identity. It can be specialized in two slightly different forms

∇ • (Ab) = (∇ • A) b + A : ∇b (2.43) ∇ • (cb) = c (∇ • b) + b∇c (2.44)
where A, b, and c are a generic matrix, vector, and scalar respectively. Such identities can be coupled with the divergence theorem to obtain the following set of equations Simulation of impact response using a frequency approach

Ω i A : ∇bds = ∂Ω i (Ab) nds - Ω i (∇ • A) bds (2.45
)

Ω i b∇cds = ∂Ω i (cb) nds - Ω i c (∇ • b) ds.
(2.46)

The strain energy E S can be developed using these two equations. In particular, the first term of the part related to in-plane stresses and displacements becomes

Ω N H : Eds = ∂Ω (N n) H vds - Ω (∇ • N) H vds (2.47)
where the following equality is implicitly used

N H : E = N H : ∇v (2.48)
which holds since N is symmetric by definition. The last term can be substituted using the in-plane equilibrium Equation (2.12)

Ω N H : Eds = ∂Ω (N n) H vds + Ω ρhω 2 v H vds+ + Ω g H αβ vds - Ω R(∇ • M)vds (2.49)
In the same way, the first part related to the out-of-plane stresses and displacements of the strain energy E S using Equation (2.45) becomes

- Ω M H : ∇∇wds = - ∂Ω (M n) H ∇wds + Ω (∇ • M) H ∇wds (2.50)
The last term can be further developed introducing the identity expressed in Equation (2.46)

Ω (∇ • M) H ∇wds = ∂Ω [(∇ • M) n] H wds - Ω [∇ • (∇ • M)] H wds (2.51) Using Equation (2.13) - Ω M H : ∇∇wds = ∂Ω [(∇ • M) n] H w -(M n) H ∇wds+ + Ω (R : N) H w + ρhω 2 w H w + g H z wds
(2.52)

The first term of the right side relative to the integral along the boundary can be further developed. In fact, by definition

∂Ω [(∇ • M) n] H w -(M n) H ∇wds = - ∂Ω m H w, n ds+ + ∂Ω [(∇ • M) n] H w -t M n H w, t ds (2.53)
Simulation of impact response using a frequency approach

The term -∂Ω t M n H w, t ds can be integrated by parts

- ∂Ω t M n H w, t ds = -t M n H w , t s=0 s=0 + ∂Ω t M n , H t wds (2.54)
The term t M n H w , t is a continuous function along boundaries that jumps on corners. Thus

-t M n H w , t s=0 s=0 = -q H C w C ∂∂Ω (2.55)
The second part of the term relative to the out-of-plane stresses and displacements using Equation (2.45) becomes where L exg is the virtual work of distributed loads g.

- Ω M H : REds = - ∂Ω R M H n vds + Ω R (∇ • M) H vds (2.

Shallow shell approximations

Shallow shells are commonly used in many engineering applications. Since any shell can be sub-divided in small shallow shell elements, these approximations provide a general and flexible albeit simple theory.

Definition 3 (Geometry approximation). The intrinsic shallow shell geometry is identical to its underlying plane geometry. Broadly, shallow shells are considered as plates where some curvature corrections are introduced. For the sake of clarity, we operate the following change of variables Simulation of impact response using a frequency approach α ≈ x (2.59)

β ≈ y (2.60)
to underline that the coordinate system is Cartesian since the underlying area is flat by definition. In fact, x and y are Cartesian coordinates of the underlying area and x and ŷ their respective unit direction vectors. The shell mid-surface is defined as

r = r(x, y) (2.61)
There is no ambiguity since there is a bijection correspondence between mid-surface points and their rectangular projections on the underlying plane. According to [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF]

, a shell is considered shallow if r, 2 x ≤ 0.05 (2.62) r, 2 y ≤ 0.05 (2.63)
Therefore an angle of 0.224 rad (≈ 13 o ) is the limit angle between the tangent plane of the mid-surface and its referring plane. In this case, Lamé surface parameters L x , L y and curvature radii R x , R y are approximated as

L x ≈ 1 (2.64) L y ≈ 1 (2.65) R x ≈ r, -1 xx (2.66) R y ≈ r, -1 yy (2.67)
Definition 4 (Kinematic approximation). In-plane displacements (u and v) are neglected in the kinematic expressions for curvature and twist changes. Mathematically, Equation (2.16) becomes

M = D : (∇∇w + RE ) (2.68)
Definition 5 (Stress approximation). Transverse shear forces are neglected in the in-plane equilibrium Equation (2.12) and in the definition of b.

Mathematically,

∇ • N- $ $ $ $ $ R (∇ • M) + g αβ + ρhω 2 v = 0 (2.69) b e = (N e -$ $ $ R e M e ) ne (2.70)
Simulation of impact response using a frequency approach

The previous approximations greatly simplify the theory since the equilibrium equations Equations (2.13) and (2.69) are partially uncoupled and the coordinate system is always Cartesian. The strain energy E S becomes

E S = 1 2 n e ∑ e=1 Ω
N H : Eds = ∂Ω (N n) H vds + Ω ρhω 2 v H vds+ + Ω g H αβ vds - $ $ $ $ $ $ $ $ Ω R(∇ • M)vds (2.72)
Coupling this result with Equations (2.68) and (2.70), leads to an equation of the virtual work theorem which is formally equal to Equation (2.58) albeit b is modified.

Simulation of impact response using a frequency approach Chapter 3

The ray-VTCR applied to shallow shell theory

This Chapter presents the VTCR applied to shallow shells. Since any shell can be considered piecewise shallow, these assumptions provide great flexibility. First, the theory is enunciated demonstrating existence and uniqueness properties of the VTCR solution. After, some improvements regarding theory and performances are introduced. For the theory, we add:

• in-plane waves,

• particular solutions, For performances we propose:

• a quasi-symmetric ray distribution method,

• a study on p-h-and ph-refinements,

• the iterative solvers.

In the end, pros & cons of the VTCR applied to the shallow shell theory are discussed. The computer program TAPYROSS is benchmarked on some simple academic examples.

Simulation of impact response using a frequency approach

Shape functions

The reference shallow shell theory and notation is presented in Section 2.5. As introduced in Section 1.2, the VTCR is a Trefftz method. Equilibrium Equations (2.13) and (2.69) are identically satisfied by the particular chosen shape functions while boundary and corner conditions are only met in weak sense. Previous works (see Section 1.2) suggest three shape function sets1 that distinguish the three VTCR versions:

• ray-VTCR,

• sector-VTCR,

• Fourier-VTCR.

The chosen VTCR version was the ray-VTCR since:

• the theory is very versatile and can handle propagative and evanescent waves at once;

• the literature about the ray-VTCR applied to thin shells is huge;

• this particular shape function set allows to analytically compute weak form integrals if some geometric conditions are satisfied (see Section 3.1).

For the sake of clarity, let us recall the generic VTCR solution formula first exposed in Equation (1.27) for 2D problems specialized for shallow shell structures

u(θ) = 2π 0 a(θ)e ık(θ) r rel dθ + u p (3.1)
where u is the displacement field, k is the wave vector, u p is a particular solution discussed in detail in Section 3.3, r rel = {x, y} is a spatial variable relative to the sub-domain, and a is an amplitude vector. It is discretized as

a(θ) = ∞ ∑ i=1 a i δ(θ -θ i ) (3.2)
where δ(θθ i ) are Dirac delta functions and a i is the discretized version of a(θ). Actually, the series is truncated to permit numerical computation. Substituting Equation (3.2) in Equation (3.1) leads to

u h = ∞ ∑ i=1 a i e ık i r rel (3.3)
where the subscript h refers to the homogeneous part of u, k i is the discretized version of k(θ). Dirac deltas effectively erase the integral in Equation (3.1) greatly simplifying weak form computations. k i and a i are divided in

k i = k i ki (3.4) a i = a i âi (3.5)
where ki and âi are unit direction vectors of the wave vector k i and the amplitude a i respectively, k i is an associated wavenumber determined by the equilibrium equations, and a i is an unknown scalar variable set by the weak form. Sections 3.1.1 and 3.1.2 describe the determination of k i , ki , and âi . Conversely, the discussion about a i is postponed to Section 3.2.

3.1.1 The unit direction vector of the wave vector [Ladevèze 1996] proves that both propagative and evanescent plane waves are required.

Their qualitative behaviors are illustrated in Figure 3.1. The difference between them lies in the definition of the unit direction vector ki . It is a 2D vector imposed a priori by the discretization process. Propagative waves are defined as

ki = T i • 1 0 (3.6)
where T i is a rotation matrix

T i = cos(θ i ) -sin(θ i ) sin(θ i ) cos(θ i ) (3.7)
and θ i ∈ [0; 2π).

Simulation of impact response using a frequency approach Evanescent waves use the following formulation as suggested in [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]]

ki j = T i • 1 + cos 2 (φ j ) ı cos(φ j ) (3.8)
where φ j ∈ [-π; π] is a different discretized angle in the complex domain. This angle controls the ratio r eo between oscillatory and evanescent parts of the wave. In this case it is

r eo = cos(φ j ) 1 + cos 2 (φ j ) (3.9)
By definition r eo ∈ [-1 2 ; 1 2 ], therefore, if faster-oscillating evanescent waves exists, they are omitted. In fact, there is no a priori consideration that allows such restriction. For this reason, we modified the formula

ki j = T i • cosh(φ j ) ı sinh(φ j ) (3.10)
In this case the ratio r eo is

r eo = sinh(φ j ) cosh(φ j ) (3.11)
This new version has no ratio limits while keeping ki j (θ i , φ j ) a unit vector. Such characteristic erases numerical difficulties that arise during the wavenumber determination process. This is further discussed in Section 3.1.2.

An important simplification in the evanescent wave discretization can be introduced when straight boundaries or polynomial sub-domains are involved. It is important because the vast majority of sub-domain boundaries lies in this category.

Definition 6 (Straight boundary). A boundary is straight if is true at least one of the following sentences:

• it is along one of the principal directions of the local Lagrange coordinate system defined on the sub-domain,

• it is straight and the local coordinate system of the sub-domain is Cartesian.

Definition 7 (Polygonal sub-domain). A sub-domain is polygonal if its perimeter is composed of straight boundaries solely.

Evanescent waves in polygonal sub-domains are crucial only along specific boundaries due to their intrinsic evanescent behavior. In particular their effects can be always neglected but when their wave vector orthogonally intersects a boundary. For this reason, all evanescent waves that do not point towards a boundary can be discarded. Therefore, in polygonal sub-domains the index i of the evanescent waves is restricted to i ∈ [1, 2, . . . , n b ] Simulation of impact response using a frequency approach where n b is the boundary number. This drastically reduces the number of DoFs. Moreover, straight boundaries permits very fast weak form integral calculations further reducing computational costs as discussed in Section 4.6. Section 3.6 exposes in detail the used distribution algorithm of ki .

The wavenumber and the unit direction vector of the amplitude

The wavenumber k i and the unit direction vector âi are set so that the equilibrium Equations (2.12) and (2.13) are identically satisfied. Let us focus on the wavenumber. This complex scalar is determined injecting its related shape function (where ki was set a priori) into Equations (2.12) and (2.13) with homogeneous conditions since u p addresses the particular solution. This operation leads to a linear set of equations that in matrix form is

Z i âi = 0 (3.12)
where Z i = Z i (k i ) depends on the wavenumber. Since âi is a unit vector, then

det[Z i ] = 0 (3.13)
This is called the dispersion equation. It sets the wavenumber k i relative to the unit direction vector ki . The other unit direction vector âi is determined injecting k i in Equation (3.12) and imposing that:

âi ∈ ker(Z i ) (3.14) âi = 1 (3.15)
where is the euclidean norm and ker( ) is the kernel operator. A solution exists because the injected k i satisfies the dispersion equation. Without approximations the dispersion Equation (3.12) is an eighth-degree equation that has eight solutions. Half of them can be discarded since concern regressive waves already addressed by the opposite unit direction vector ki = -ki . Three of the remaining four wavenumbers has ℜ{k i } >> ℑ{k i } while one of them has ℜ{k i } << ℑ{k i }. Wave types are identified combining the characteristics of the wavenumber k i , its related unit direction vector âi = [a xy , a z ] , and the introduced unit vector type ki as illustrated in Table 3.1.

In previous works such as [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]] equilibrium equations were approximated to simplify the dispersion equation. In particular, the in-plane inertia was neglected. As side effect, even in-plane waves were neglected. Broadly, the group velocity in homogeneous and isotropic plates is by definition Simulation of impact response using a frequency approach Wave type k i and âi characteristics Propagative in-plane

   propagative ki ℜ{k i } >> ℑ{k i } a xy >> a z Propagative out-of-plane    propagative ki ℜ{k i } >> ℑ{k i } a xy << a z Evanescent out-of-plane    evanescent ki ℜ{k i } << ℑ{k i } a xy << a z
c g ip = E ρ(1 -ν 2 ) (3.16) c g oop = ω 4Eh 2 3ρ(1 -ν 2 ) (3.17)
where c g ip and c g oop are respectively the in-plane and the out-of-plane group velocities. At low-and mid-frequency c g oop << c g ip . Their decay ratios in plates as well as shells qualitatively follow the same rule. Since in [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]] the focus was on the stationary behavior, neglecting in-plane waves was completely legit. This simplification leaded to a fourth-degree dispersion equation that produced an evanescent and a propagative wave. Both of them were out-of-plane. Hence, the wave identification process was unnecessary and the DoF number was drastically reduced since there were no in-plane waves. In the HSS3+ test the focus is on the transient dynamics which is dominated by in-plane waves. For this reason, in-plane waves should be addressed and the simplification is inappropriate to our case. In Sections 3.4.3 and 3.4.4 corrections for in-plane inertia are benchmarked on some relevant academic tests. Numerical difficulties can arise during wavenumber computation if the direction vector introduced in the dispersion equation is not unitary. In that case

k i = k i ki (3.18) ki = c ki (3.19) ki = k i c , (3.20)
since the direction vector k i remains unchanged. ki is the modified wavenumber, ki is the direction vector introduced in the dispersion equation, and c = ki is its norm. Eventually, the product c ki should be performed to compute k i . When c >> 1 (or c << 1) this

Simulation of impact response using a frequency approach operation could be inaccurate due to truncation error. For this reason, the direction vector introduced in the dispersion equation should be unitary.

The weak form

This Section describes the weak variational formulation and its key properties

D d = {u d , N d , M d } = {u s2 -u s1 , N s2 -N s1 , M s2 -M s1 } is
+ w H C ed q C ed ∂∂ q C ed Ω e -q H C ed w C ed ∂∂ w C ed Ω e + n Γe n Γe + 1 Γe b H ed v ed + q H ed w ed -m H ed w ed , ne ds - n C e n C e + 1 q H C ed w C ed C e - 1 n Γe + 1 Γe v H ed b ed + w H ed q ed -w ed , H ne m ed ds + 1 n C e + 1 w H C ed q C ed C e = 0 ∀ {u ed , N ed , M ed } ∈ D ed (3.
D s1 = D s2 (3.27)
The existence property directly descends from the uniqueness property. If the solution of the forced vibration problem respects the regularity conditions of finite energy vector fields enunciated in Section 2.1, it is also solution of the weak variational formulation.

If the solution D s is a discrete space approximation, the uniqueness property implies the existence property.

Computational properties of the weak form

Injecting the shape functions in Equation (3.21) produces a linear set of equations which is in matrix form

Ba = l (3.28)
where B, l, and a are related to the bilinear form, the linear form, and the amplitude vector respectively. a is composed by the unknown variables a i . The matrix is always square because the number of shape functions is equal to the number of test functions being the VTCR a Trefftz method. The uniqueness property coupled with η > 0 (by definition) implies that [ℑ {B M }] sym is positive definite being B M the main diagonal block matrix of B. The matrix related to the bilinear form is neither symmetric nor skew-symmetric and presents ill-conditioning problems. Section 3.7 further discusses these problems and proposes a way to mitigate their effects. The bilinear form is a block matrix as illustrated in Figure 3.2. A generic block B e e considers the effects of the shape functions of the sub-domain Ω e along boundaries and on corners of Ω e by means of the test functions of the sub-domain Ω e .

Shape functions (which are test functions too) defined on Ω e produce effects only along ∂Ω e and on ∂∂Ω e by definition. Hence, when ∂Ω e ∩ ∂Ω e = 0 and ∂Ω e ∩ ∂Ω e = 0, the sub-matrices B ee and B e e are null. In general, a block matrix B e e is a rectangular matrix n t × n s where n s is the shape function number of Ω e and n t is the test function number of Ω e . Main diagonal blocks are always square matrices since test and shape Simulation of impact response using a frequency approach Every quantity can be extrapolated from the generic shape or test function along a boundary using Equations (2.32) to (2.35). Since shape and test functions are plane waves, the generic quantity related to a shape function u j is

= L u j (3.31)
where L [ ] is an operator relative to . In particular it depends on a and k of the introduced function and the outward unit normal vector n, and the curvilinear coordinate along the boundary s. The same reasoning can be performed for a test function δu i . Let us denote the dual quantity of . Integrals are always in the following form

∂ Ω e ∩∂ Ω e L [δu i ] H L δu j ds (3.32)
Simulation of impact response using a frequency approach

Corners require just one condition either w C or q C since w C is the dual quantity of q C and vice versa. Be C either w C or q C and C its dual quantity and be L C [ ] and L C [ ] their respective operators. The corner quantity produced by u j and δu i on corners is Matrix construction could be very expensive in terms of computational costs since the weak form requires a large amount of integrals that, in general, should be evaluated numerically. Typically, it is the computation bottleneck. As mentioned before, along straight boundaries integrals can be computed analytically. Section 4.6 focuses on analytic integrals illustrating in detail a new effective algorithm to compute and stock them. It drastically reduces computational costs.

L C [δu i ] H L C u j (3.

Particular solutions for general surface loads

The present Section exposes a very general method to address surface loads of any sort. Previously, only particular surface load types could be addressed.

The particular solution u p of Equation (3.1) part arises if distributed loads g are present. Even if the HSS3+ test does not require a particular solution since g = 0 everywhere, the particular solution component is important to benchmark various VTCR improvements. In fact, in literature is present an analytic solution which is treated in detail in Section 3.4.4.

The present method addresses every possible distributed external load g given that g satisfies some minor assumptions. Almost every real scenario can be tackled. The focus is on the inhomogeneous Equations (2.13) and (2.69) where g = g(x, y). Let us suppose that g : R 2 → R 3 can be approximated with a 2D Fourier series:

g = ∑ i, j∈Z
g Fi j e ık Fi j r rel (3.34)

g Fi j = 1 d x d y d x /2 -d x /2 d y /2 -d y /2 ge -ık Fi j r rel dxdy (3.35) k Fi j = 2π i d x , j d y (3.36)
where the subscript F is related to the Fourier series, d x and d y are characteristic dimensions of a 2D sub-domain Ω Fe , and r rel is a 2D position vector relative to the geometric Simulation of impact response using a frequency approach center of Ω Fe . The sub-domain Ω Fe is the 2D space portion over which the Fourier series approximates g. For this reason, it is

Ω e ⊆ Ω Fe (3.37)
In case of flat or cylindrical sub-domains it is a rectangle or a cylinder sector respectively. For computational reasons, the series in Equation (3.34) is symmetrically truncated around zero. Since the distributed load is a Fourier series and the equilibrium Equations (2.12) and (2.13) are a linear set of differential equations, the solution of the inhomogeneous problem becomes

u p = ∑ i, j∈Z
u pFi j e ık Fi j r rel (3.38)

Z Fi j u pFi j = g Fi j (3.39)
where Equation (3.39) is the matrix form of the equilibrium equations for each injected g Fi j and the series in Equation (3.38) is actually truncated. In particular, Z Fi j is a 3 × 3 complex matrix. For computational purposes, the indexes i and j are summarized in i s . In this way, u p is approximated as

u p ≈ n i s ∑ i s =1
u pFi s e ık Fi s r rel (3.40)

where n i s is a truncation parameter that takes into account the truncation parameters of i and j. The created solution u p satisfies the equilibrium equations but not boundary, corner, and coupling conditions. The weak form used to determine the homogeneous solution u h addresses residues of u p along boundaries and on corners. In particular, a generalized linear form l p between the Fourier functions and the test functions of Ω e is computed. The final linear set of equations in Equation (3.28) becomes:

Ba = l -l p (3.41)
For the sake of generality, the term l p can be further developed as

l p = B p c p (3.42)
where B p is a generalized bilinear form between the Fourier functions e ık Fi s and the test functions of Ω e and c p is a complex vector. Its lengths is equal to the number of Fourier functions. The generic element c pi s of c p is series of products and multiplication between the elements of u pFi s and k Fi s depending on the particular boundary and coupling conditions applied along ∂Ω e and on the corner and coupling conditions applied on ∂∂Ω e .

We have to remark that generalized bilinear form B p is generated computing integrals between the Fourier functions of u p and the shape functions of the homogeneous solution.

For this reason, B p is not strictly a Galerkin operation. All routines for the computation Simulation of impact response using a frequency approach of B can be recycled to compute B p since both Fourier functions of u p and test functions of u h are imaginary exponentials. This last step allows to compute B p once. When the particular solution changes, only c p and a should be rebuild. In particular, B, B p , and l remains unchanged. This drastically reduces computational costs of successive computations since matrix construction is a relevant part of total computational costs. In this way the VTCR can address every possible particular solution given that its Fourier series approximation exists. In real scenarios such assumption is always true.

Numerical tests

In order to present consistent test results, error indicators, software, and convergence criteria are alike among different tests. For this reason, they are introduced in the following paragraphs.

Error indicators

They are based on the displacement field since it unequivocally defines the solution. In order to summarize results in one relevant visual comparison, displacement magnitude portraits are investigated. In the following Sections two comparison types are performed:

• a VTCR solution with a FEM reference,

• two different VTCR solutions.

The first one is a cross-method confrontation while the second one is not. At midfrequency a small difference in theories can lead to slightly different frequency responses. Since at this frequency range there are many wave lengths per sub-domain, a small difference in wavenumbers can lead to different displacement magnitude portraits that have almost the same energy. For this reason, two different error indicators are defined for each comparison type. Both are based on the kinetic energy E K defined in Equation (2.39). Direct comparisons of VTCR and FEM displacement magnitude portraits could be non-optimal due to small theory differences that can lead to different solutions at mid-frequency. Therefore, in this case, an error indicator based on a comparison between total energies is used

err FEM = |E K (u FEM ) -E K (u V TCR )| E K (u FEM ) , (3.43)
where u FEM and u V TCR are displacement fields of FEM and VTCR respectively and | | denotes the absolute value.

In order to either highlight theoretical improvements or study VTCR convergence, a VTCR -VTCR comparison is required. In this case a more strict error indicator can be Simulation of impact response using a frequency approach defined since the two solutions share the same theory. It is based on the energy of the displacement magnitude portrait difference

err V TCR = E K (u V TCR m -u V TCR n ) E K (u V TCR m ) , (3.44)
where u V TCR n and u V TCR m are two different VTCR solutions, n and m denote ray numbers, and m > n.

Software and convergence criteria

VTCR solutions are compared with ABAQUS R tests which are considered FEM references. In low-frequency range the rule-of-thumb to set mesh size is

h 0 = λ 10 , (3.45)
where h 0 is the maximum edge length of a mesh element and λ is the typical phenomenon wavelength. Equation 3.45 can be rearranged in terms of wavenumber k

h 0 k = 2π 10 . (3.46)
Ihlenburg in acoustic [Ihlenburg 1998] and Deraemaeker, Babuška, and Bouillard for general Helmholtz problems [Deraemaeker 1999] proved that this relation is not valid at mid-frequency due to high-scattering behavior. In particular, they affirmed that pollution error becomes predominant as the wavenumber increases. They suggested a corrected version of the rule-of-thumb

h 2 0 k 3 = 4π 2 100L (3.47)
where L is a characteristic dimension of the problem considered. Since in shell theory many wave types are present, the following rule-of-thumb is used

h 2 0 k 3 MAX = 4π 2 100L , (3.48)
where k MAX is the greatest wavenumber. By definition, in shells in-plane stiffness is much greater than out-of-plane stiffness; thus k MAX is always the wavenumber of the propagative out-of-plane wave along the direction relative to the highest Lamé parameter. VTCR is implemented in MATLAB R . In the following tests VTCR convergence is studied using Equation 3.44 where u V TCR m is an initial guess. The rays are ordered by importance:

• propagative rays of the particular solution,

• propagative out-of-plane rays, Simulation of impact response using a frequency approach

• evanescent out-of-plane rays,

• propagative in-plane rays.

In fact, the propagative rays of the particular solution define the VTCR input. For this reason, they are crucial. After, since out-of-plane waves are the fastest oscillating ones, these rays are the most scattering and should be considered first. Anyhow, every other sequence is possible. For the sake of simplicity, the number of propagative normal and shear stress carrying rays are kept alike. Hereafter, the displayed number refers to their sum. Therefore, it is always an even number. Moreover, all the displayed VTCR DoF numbers are implicitly intended per sub-domain. A posteriori convergence-seeking is a five-step process:

1. A VTCR solution is computed with an initial guess of the ray numbers. Ray numbers (per element and per type) are intentionally overestimated. However, in case the discretization is still too coarse, new DoFs can be added without recomputing matrices. This matrix-recycling property is related to a new quasi-symmetric ray distribution algorithm. It is discussed in Section 3.6.

2. Propagative rays of the particular solution are increased up to the chosen convergence criterion (based on err V TCR ) keeping the other ray types equal to the reference.

3. Propagative out-of-plane rays are increased until convergence keeping the number of propagative rays of the particular solution equal to their converged value while evanescent and in-plane rays are equal to the reference.

4. Evanescent out-of-plane rays are increased up to convergence. As before, propagative rays are kept equal to their converged values while the remaining ray types are set equal to the reference.

5. Finally, propagative in-plane rays are investigated keeping the other ray types equal to their converged values.

Converged rays are used after the first two steps since interactions between rays can slightly modify the convergence error. As side effect, this yields plateaus in the convergence graphics at the level of the error threshold chosen. This process is cheap in terms of computational time because matrices can be reused thanks to the quasi-symmetric ray distribution algorithm. In fact, once matrices are computed for the initial guess they are reused to analyze convergence.

The same workstation is used to perform all the tests (VTCR and FEM) to produce consistent performance data. Its characteristics are reported in Figure 3.3 presents geometry and boundary conditions. Material properties, frequency, and shell thickness are reported in Table 3.3. The vibrational problem is depicted in .

Since no analytic solution is known in this case, the VTCR solution is compared with a FEM reference. The VTCR convergence curve is reported in Figure 3.4. Since no outof-plane waves are required, the convergence process described in Section 3.4.2 can be reduced to just one step. Table 3.4 summarizes VTCR convergence results. We consider that the solution is at convergence when err V TCR ≤ 0.01. Figure 3.5 illustrates VTCR and FEM displacement results. Table 3.5 compares VTCR and FEM performances as well as DoFs number at convergence. In this particular academic test, err FEM ≤ 10 -6 at convergence. Finally, Figure 3 Simulation of impact response using a frequency approach Remark. The present academic example proves that in-plane waves are taken into account. At the same time, some boundary conditions are verified to work as intended.

VTCR and FEM performances are comparable even if the computed frequency is relatively low. In Section 4.6 a fast computation algorithm is presented to drastically improve VTCR performances. Finally, Figure 3.6 demonstrate that the VTCR solution is unaffected by the particular domain division chosen. Conversely, performances degrade as the number of sub-domains increase. Such behavior is investigated in Section 3.7.1 in detail.

Numerical test: a punctual force applied on a shallow shell

This academic example presents a shallow shell subject to a punctual force on the surface.

It is relevant because in this particular case a quasi-analytic solution is known [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF]. It is quasi-analytic since it is based on a Fourier series (which should be truncated for computation purposes) and in-plane inertia is neglected. According to [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF], it is legit in this case. Since in this case the reference is not based on FEM, some interesting studies can be performed. The present Section proposes a triple comparison between the VTCR, the FEM, and the quasi-analytic reference to further validate TAPYROSS. Furthermore, the present example studies the differences introduced by the shallow shell assumptions since the FEM reference uses the complete shell theory. Geometric dimensions and boundary conditions are depicted in Figure 3.7. Some

Simulation of impact response using a frequency approach important quantities are reported in Table 3.6.

g u = 0 v = 0 N x = 0 N y = 0 m = 0 w = 0 w = 0 m = 0 u = 0 N y = 0 w = 0 m = 0 v = 0 N x = 0 m = 0 w = 0 (
[ [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF] reports the reference solution. The introduced approximations are:

• shallow shell assumptions,

• the in-plane inertia is neglected,

• there are no in-plane point loads.

It uses the Airy's stress function ϕ. It acts like a potential field of the in-plane stress resultant tensor N. It is defined as

h thickness 3 mm R x radius of curvature along x 2.0304 m R y radius of curvature along ŷ ∞ f
∇ 2 ∇ 2 ϕ + Eh∇ 2 k w = 0 (3.52) -∇ 2 k ϕ + D∇ 2 ∇ 2 w = g z -ω 2 ρhw (3.53)
where ∇ 2 = , xx + , yy is the lagrangian operator and

∇ 2 k = , xx R y + , yy R x .
The point load is approximated with a Fourier series

g z = ∞ ∑ i=1 ∞ ∑ j=1 C i j sin(λ i x) sin(µ j x) C i j = 4 d x d y d x d y δ(x -p) sin(λ i x) sin(µ j x)dydx λ i = iπ d x µ j = jπ d y (3.54)
where δ(x-p) is a Dirac delta function. The reference solution in the form of {ϕ(x, y), w(x, y)} is

a λ i µ j = (λ 2 i + µ 2 j ) 2 (3.55) b λ i µ j = ( λ 2 i R y + µ 2 j R x ) 2 (3.56) B i j = Ehb λ i µ j a λ i µ j + Da λ i µ j -ρhω 2 -1 C i j (3.57) A i j = Ehb λ i µ j a λ i µ j B i j (3.58) w = ∞ ∑ i=1 ∞ ∑ j=1 B i j sin(λ i x) sin(µ j x) (3.59) ϕ = ∞ ∑ i=1 ∞ ∑ j=1 A i j sin(λ i x) sin(µ j x) (3.60)
Simulation of impact response using a frequency approach Table 3.7: Ray number and type for initial guess VTCR and converged VTCR. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.4 where the error threshold is err V TCR ≤ 0.01.

The performance comparison between VTCR and FEM is reported in Table 3.8. The propagative rays of the particular solution are not taken into account in the VTCR DoF number since they are determined a priori and are an input of the homogeneous problem. The FEM DoFs number is set by the convergence criterion exposed in Section 3.4.2 since the FEM is affected by pollution [Deraemaeker 1999] in the mid-frequency range.

The triple comparison is depicted in Figure 3 difficult, these displacement fields are omitted.

Remark. Boundary conditions are always met regardless of the particular method. Moreover, the VTCR solution perfectly matches the reference. Conversely, the FEM solution is slightly different due to the shallow shell approximations which are not introduced in the FEM study. However, displacement and stress resultant levels are of the same order of magnitude. It is interesting to point out that at higher frequencies these differences increase since there are more wavelengths per sub-domain. This suggests that the shallow shell approximations may produce unacceptable errors in the mid-frequency range. This, along with other shallow shell drawbacks are further discussed in Section 3.7.1.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach

Numerical test: complex frame structure

The present academic test studies a complex frame structure. A triple-joint structure is investigated and VTCR solution compared with a FEM reference. Geometric dimensions and boundary conditions are reported in Figure 3.10 while the other important quantities are stated in Table 3.9. The test is performed at low-frequency since FEM computational costs are prohibitive at mid-frequency.

The multi-step convergence analysis described in Section 3.4.2 is depicted in Figure 3.11. The error threshold is set to err V TCR ≤ 0.01. Figures 3.11b and 3.11c prove that in this case in-plane and evanescent waves are required. Its results are reported in Table 3.10 as well as the ray numbers and types of its initial guess.

FEM and VTCR solutions are compared in Figure 3.12. Their computational costs as well as DoF numbers are reported in Table 3.11.

Remark. Figure 3.12 confirms that the VTCR correctly addresses the three-shell-joint. In fact, the VTCR-FEM error defined in Section 3.4.2 is in this case err FEM ≈10%. The small remaining difference is due to the shallow shell approximations introduced in the VTCR that are absent in the FEM reference. The performance comparison in Table 3.11 assert that the VTCR and FEM perform alike at low-frequency. Since the pollution effect is negligible at this frequency, FEM time performance greatly improves reaching its VTCR counterpart. Yet, VTCR still outperforms FEM of one order of magnitude in terms of memory consumption even at low-frequency. This is due to the great DoF number difference between these two different approaches. Table 3.11: DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 3.4.5.

h
Simulation of impact response using a frequency approach

Performances

In Section 3.4 some performance improvements are implicitly used. Even if some of them are essential to study a numerical problem, their description is postponed in the present Section for the sake of classification and clarity. Section 3.5.1 analyses VTCR performances over a frequency band in comparison with a FEM reference. Section 3.6 describes the quasi-symmetric ray distribution algorithm already used in Sections 3.4.3 and 3.4.4 to study VTCR convergence with negligible computational efforts. Section 3.7 proposes iterative solvers to mitigate the ill-conditioning behavior of the VTCR matrices. Section 3.7 studies three iterative solvers on some numerical tests to find the best suited for the VTCR. The VTCR refinement process described in Section 3.4.2 is based on the assumption that increasing the shape function number is better than the sub-domain number. [START_REF] Riou | The multiscale VTCR approach applied to acoustics problems[END_REF] demonstrates such property in the acoustic domain. Section 3.7.1 studies the two approaches on a numerical test to determine which one performs better in the shallow shell theory.

Numerical test: performances over a frequency band

The present Section investigates VTCR performances over a frequency band [2000; 8000] Hz comparing VTCR and FEM results. The considered academic problem is the shallow shell subject to a point force exposed in Section 3.4.4. In fact, in this case the known quasi-analytic solution provides a neutral comparison reference. Since it is computed using neither the VTCR nor the FEM, it is affected by neither pollution effects nor illconditioning problems. In order to measure performances degrees of freedom (DoFs) cannot be directly compared since methods are intrinsically different. Nevertheless, for the sake of completeness, they are reported with time and memory consumption in Remark. At mid-frequency FEM suffers of pollution effect as proved in [Deraemaeker 1999, Ihlenburg 1998]. DoFs, time, and memory required increase faster as frequency grows. Conversely, VTCR is unaffected by pollution effect as illustrated by Figure 3.13. Since VTCR remains stable, computational cost differences grow with frequency. At 2000 Hz consumption differences are already of one order of magnitude and increase with frequency. Due to low memory requirements VTCR could have been run on a laptop while FEM needs a workstation. Moreover, the difference in computational time is so massive that VTCR on a laptop would still have outperformed FEM on a workstation. The fact that these results are obtained comparing an unoptimized MATLAB R program against a FEM implemented in the commercial code ABAQUS R confirms even further VTCR effectiveness.

Simulation of impact response using a frequency approach

Quasi-symmetric ray distribution and matrix recycling

Precedent numerical tests used the multi-step convergence process proposed in Section 3.4.2. Implicitly, a quasi-symmetric ray distribution method was employed to speed up the process. This algorithm permits to recycle matrices. Since VTCR uses an exiguous DoF number, matrix construction is a relevant (predominant in some cases) operation in terms of computational costs while matrix inversion is a relatively inexpensive step. For this reason, the present algorithm drastically reduces computational costs when matrix can be reused.

The angle θ i ∈ [0, 2π) controls propagative wave direction. The parameter φ j governs ratio r eo . Since it is discretized, φ j is distributed over [-φ , φ ] where φ is set a priori and can be always changed. For the sake of simplicity, the range is kept symmetric. In previous works [START_REF] Ladevèze | Calculation of medium-frequency vibrations over a wide frequency range[END_REF], Riou et al. 2004, Kovalevsky et al. 2014, Riou et al. 2008] θ i and φ j were distributed symmetrically over the unit circle and the half complex unit circle respectively. This choice imposes a complete matrix recomputation as the number of rays change. For example, let us suppose we computed a VTCR solution of a vibrational problem with three propagative rays and we found it unsatisfactory. We want to add one or two rays. In order to keep a symmetric ray distribution, previous rays must move as illustrated in Figure 3.14. For this reason previous matrices cannot be reused and must be recomputed anew as ray number change. Since for VTCR matrix creation can be relevant in terms of computational time, this effect can be computationally Simulation of impact response using a frequency approach expensive. Conversely, in the present approach a quasi-symmetric ray distribution is used. In this algorithm previous rays are fixed as new ones are added. The first ray can be placed in any direction. After that, new rays are inserted in gaps among previous rays in the most possible symmetric way. Figure 3.15 illustrates the algorithm for the first eight rays. This method permits complete recycling of old matrices when more rays are added. The drawback is that, for a given ray number, its distribution could be asymmetric. For this reason, this algorithm could require more rays to reach convergence than a symmetric distribution method. Let us compare a quasi-symmetric and a symmetric distribution of n r rays. The shift angle ∆θ i of ray i in the quasi-symmetric distribution with respect of the symmetric one is

∆θ i = γ i 2π n r (3.61)
where γ i ∈ [0, 1] is a discrete parameter. Since it varies with n -1 r , the difference between Simulation of impact response using a frequency approach algorithms decreases as the ray number n r increases. In practice, when convergence is reached the difference is already negligible. Moreover, in this case

F n r ⊆ F m r n r ≤ m r (3.62)
where F n r and F m r are shape function spaces with n r and m r shape functions respectively. Therefore, except for pollution and resonance effects, an error indicator is a neverincreasing function with the number of shape functions.

For the sake of clarity the present discussion was implicitly developed considering a propagative wave distribution. Yet, it can be extended to evanescent wave distribution. As explained in Section 3.1.1, there are two parameters ki j = k(θ i , φ j ) that control the unit direction vector of the wave vector. θ i indicates the direction along which the wave exponentially increases. In general it is distributed as for the propagative waves. In case of polygonal sub-domains only the θ i that point towards a boundary are retained. The reasons behind this simplification are illustrated in Section 3.1.1. φ j leads the ratio between the oscillatory and evanescent behaviors. This parameter is a real scalar that can theoretically assume any value. For practical reasons, it is limited into a symmetric range φ ∈ [-φ, φ]. The quasi-symmetric distribution is identical to the algorithm proposed in Figure 3.15 for θ i except for the distribution boundaries.

Iterative solvers

The present Section describes in detail the ill-conditioning problem focusing on its causes, its effects on the solution, and the chosen way to mitigate it. Typically, VTCR suddenly converges when ill-conditioning problems appear. For example, let us consider the academic test exposed in Section 3.4.4. Figure 3.16 illustrates this phenomenon superposing the conditioning number on the convergence graphic depicted in Figure 3.8b. In order to further investigate this correlation, the definition of amplitude portrait is introduced.

Definition 9 (Amplitude portrait). In case of propagative waves, the amplitude portrait is a polar graphic where is plotted the discrete function ap(i) = ki a i . In case of evanescent waves, the amplitude portrait plots the discrete function ap(i) = a j cos(φ j ), sin(φ j ) since the focus is usually on the ratio r eo .

In both cases a set of arrows or segments is used since the amplitude function is discretized using Dirac deltas. Figure 3.17 compares two amplitude portraits of two different discretizations (beyond convergence, 100 and 200 propagative out-of-plane rays respectively) and their respective solutions.

Even if both cases meet convergence criteria and their final solutions |w 100 | and |w 200 | are indistinguishable with the naked eye, their amplitude portraits are completely different. Even the order of magnitude of the rays changes between the two discretizations. [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF] suggests to couple VTCR with a Reduced Order Model (ROM) method to tackle wide frequency band problems reducing computational costs. Typically, these methods are very sensible to ill-conditioning problems and amplitude portrait instabilities. This astonishing behavior is related to the Fredholm integrals used in the weak form2 . [START_REF] Yeih | Numerical instability of the direct Trefftz method for Laplace problems for a 2D finite domain[END_REF]] mathematically describes this effect. Broadly, from the point of view of a computer, where cond[] and det[] are two operators that compute the conditioning number and the determinant respectively. Equation (3.63) implies that the matrix is rank deficient and the problem solution is no more unique. The present work uses iterative solvers to mitigate ill-conditioning issues. Typically, these techniques impose to the problem solution an additional minimization condition into the Krylov subspace of B. This additional requirement restores the uniqueness property of the solution.

In order to decide which iterative solver is the best suit for the VTCR, we study the amplitude portraits of the VTCR solutions of the problems presented in Sections 3.4.3 and 3.4.4 as well as another academic problem specifically proposed to highlight solver differences. The study focuses on solutions that are beyond convergence thresholds. Evanescent rays are omitted for the sake of simplicity. Since displacement magnitudes are unaffected, they are not reported. The solvers considered are:

• pinv. This algorithm dermines the Moore-Penrose pseudoinverse of matrix. It is suggested for ill-conditioning problems since it normalizes to one the smallest singular values. The result is a relatively well-conditioned pseudoinverse [Courrieu 2008].

• lsqr. It is based on the Lanczos tridiagonalization [START_REF] Paige | LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares[END_REF]. Simulation of impact response using a frequency approach

• gmres. It uses the Arnoldi's method to compute an orthonormal basis of the Krylov sub-space. The method restarts if stagnation occurs [START_REF] Saad | GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]].

• backslash. It is the standard direct MATLAB R solver. It is considered for reference.

The plate subject to an in-plane load in Section 3. Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach

q = 1 N/m Ω 2 w, n = 0 Ω 1 Ω 3 Ω 4 w, n = 0 w, n = 0 q = 0 w = 0
The shell subject to an point force in Section 3.4.4 focuses on the out-of-plane behavior. The propagative waves of the particular solution do not affect the bilinear form. Thus, the analysis focuses only on the out-of-plane propagative waves. The compared discretizations are at 100 and 200 DoFs which are both beyond convergence (65 DoFs). Figure 3.23 reports the results. Remark. The various academic example presented illustrate the effects of the solvers on the amplitude portraits. As expected, the order of magnitude of the amplitudes computed with the standard direct solver (backslash) explodes in every studied case. Conversely, the iterative solvers remains in the same range. In the plate subject to an in-plane load and in the shell subject to a point force all the iterative solvers qualitatively keep the same amplitude portrait as the rays increase. Conversely, in the plate subject to an out-of-plane load pinv and gmres behave alike while lsqr proposes solutions where horizontal waves are predominant. For this reason, lsqr seems to be the best suited method. However, further tests can revert this conclusion since it is neither definitive nor based on many tests.

Numerical test: p-and ph-refinement comparison

Until now, we implicitly supposed that increasing the number of shape functions was the best way to refine the VTCR solution. This assumption is plausible because a study in Simulation of impact response using a frequency approach [START_REF] Riou | The multiscale VTCR approach applied to acoustics problems[END_REF] supports it in the acoustic domain. However, there is no such analysis in literature for the VTCR applied to shells. For the sake of rigorosity, the present Section compares p-and ph-refinement performances. The p-refinement method keep the sub-domain number as low as possible and increases the shape function number. It was implicitly used in the multi-step convergence analysis in Section 3.4.2. Conversely, the ph-iteration is a two-steps procedure:

1. The domain is divided in n e sub-domains.

2.

A p-refinement is performed for every set-up.

The numerical test investigated is the shallow shell subject to a punctual force reported in Section 3.4.4 because the quasi-analytic reference solution increases the study precision. The shape functions of the particular solution are kept fixed for the sake of simplicity. The p-refinement was already presented in Figure 3.8. It can be considered the first step of the ph-refinement. For this reason, it is reproposed in Figures 3.25a and 3.27a for comparison purposes.

In this numerical example four configurations are chosen n e = [1, 2, 4, 8] for the phrefinement. The structure presents two different curvature radii that lead to two domain sub-division types (and their hybrid version). They are illustrated in Figure 3.24.

Ω 1 Ω 2 Ω 1 Ω 2 Ω 3 Ω 4 Ω 1 Ω 2 Ω 3 Ω 7 Ω 4 Ω 6 Ω 5 Ω 8 (a) Sub-divisions along x. Ω 1 Ω 2 Ω 1 Ω 2 Ω 3 Ω 4 Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6 Ω 7 Ω 8 (b) Sub-divisions along ŷ.
Figure 3.24: Different sub-division methods of the numerical example described in Section 3.4.4 and resumed in Section 3.7.1 to compare p-and ph-refinements.

Simulation of impact response using a frequency approach Remark. The required rays to reach convergence per sub-domain decreases as the the domain is divided. Broadly, decreasing the characteristic dimension of a sub-domain is equivalent to decreasing the frequency since there are less wave lengths per sub-domain. Therefore, if the sub-domain becomes smaller, the frequency needed to present a mediumfrequency behavior increases. Conversely, the total required number of rays increases as the domain is divided. Therefore, there is a serious advantage in convergence rate if the number of sub-domains is kept as low as possible. Such result is consistent with the study in [START_REF] Riou | The multiscale VTCR approach applied to acoustics problems[END_REF] developed in acoustic domain.

Conversely, the sub-division along x creates sub-domains that present projection planes that are slightly rotated with respect of the original projection plane as illustrated in detail in Figure 3.26. For this reason, a small angle ϑ e appears between the generic sub-domain Ω e and the successive Ω e+1 . If sub-divisions are uniform, it decreases as n e increases in this way

ϑ e = π - ϑ n e (3.64)
where ϑ is the angle between surface unit normal vectors at the extreme boundaries (x = [-0.35, 0.35] m). Shallow shell assumptions impose that ϑ ≤ 0.448 rad. For the same reason, sub-domain projection planes approach the effective shell surface as n e increases. The ph-refinement is depicted in Figure 3.27. Since in this case the solution changes as the sub-domains increase, displacements are illustrated in Figure 3.28 along with the reference solution to highlight differences.

Remark. The sub-division along x presents a serious VTCR solution deterioration as illustrated in Figure 3.28. Convergence graphics depicted in Figure 3.27 use a VTCR solution with the same sub-division as reference. They prove that the VTCR converges to a solution. Yet, it is different from the reference one. Probably, introducing an abrupt projection plane rotation generates virtual auto-equilibrated loads along the boundary that degrade the solution. These loads are related to θ e . Hence, the effect of each load decreases as n e increases. Conversely, their number increases. The example illustrates that their global effect degenerates the solution as n e increases.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach For the sake of simplicity and symmetry, the complete cylinder is divided in 16 subdomains. Sub-divisions, geometry, and boundary conditions are illustrated in Figure 3.30 while material properties and the studied frequency are reported in Table 3.12.

Figure 3.31 compares the VTCR solution with a FEM reference. The VTCR displacement magnitude is a reference for the multi-step convergence process. Further refinements (omitted for the sake of simplicity) illustrate that the VTCR is already at convergence. Since there is no correlation between VTCR and FEM solutions, convergence and performance studies are unnecessary. In this case the VTCR completely failed to determine the correct result. Table 3.13 reports ray number and type.

Remark. Since the material is homogeneous and isotropic and loads and geometry are axisimmetric, the solution should be axisimmetric too. The FEM reference reported in Figure 3.31b confirms it. Similarly to the refinement along the x axis in Section 3.7.1, the θ ψ ψ

Surface Underling plane

Figure 3.29: Section of a generic cylinder sector and its underling plane. Angles are trigonometrically related since θ = 2ψ.

Simulation of impact response using a frequency approach VTCR solution is heavily degraded. Probably, the auto-equilibrated loads along boundaries are predominant and the real solution, if correctly computed, is hided. However, this is only a possible explanation. This example is presented only to highlight the problem.

N x = 0 m = 0 N y = 1 N/m q = 1 N/m u = 0 m = 0 3.5 m (a) Boundary conditions. 2.75 m Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6 Ω 7 Ω 8 Ω 9 Ω 10 Ω 11 Ω 12 Ω 13 Ω 14 Ω 15 Ω 16 (b) Sub-division.
In particular, we do not pretend to mathematically explain such effect since further studies are required.

Conclusions

On one hand, the shallow shell approximations provide many advantages to the VTCR. The equations are simpler. This can be a key feature for the future VTCR developments required to analyze the HSS3+ test. Previously [START_REF] Riou | Extension of the variational theory of complex rays to shells for medium-frequency vibrations[END_REF]], the weak form required a tuning phase if the coordinate system was not Cartesian (i.e. cylindrical). This is no more necessary with the present theory since the surface is projected to the underling plane and the geometry is always described in a Cartesian coordinate system. Any shell can be divided in shallow shell sub-domains. This property provides great flexibility to the VTCR.

On the other hand, some serious drawbacks arise. Section 3.7.1 demonstrated that the VTCR performs best when the sub-domain number n e is kept as low as possible. In general, the shallow shell theory requires a finer domain sub-division than the shell theory. As explained in Section 3.7.2, 16 sub-domains are necessary to analyze a complete cylinder. The same problem analyzed with the shell theory requires just 4 sub-domains3 . Moreover, Sections 3.7.1 and 3.7.2 illustrated that the shallow shell theory does not produce good results if the sub-divisions are along a curved boundary. Probably, it is due to virtual auto-equilibrated loads among sub-domains produced by the angles between adjacent projected planes. Despite the fact that each load decreases as the domain is sub-divided, the solution degrades since the boundary number (and the number of virtual auto-equilibrated loads) increases.

For this reason, in Chapter 4 we return to the VTCR applied to the shell theory. All the shallow shell improvements are incorporated in the theory along with many more enhancements.

Chapter 4

The ray-VTCR applied to shell theory

The present Chapter develops the VTCR version for the general shell theory which is included in TAPYROSS. Improvements already introduced for the shallow shell theory are also taken into account. Further enhancements are presented to tackle the industrial test HSS3+. They concern:

• orthotropic and sandwich materials,

• conic structures and shells of variable thickness,

• moving loads,

• concentrated masses and springs.

Finally, a new algorithm that drastically increase performances is introduced to compute the weak form integrals. Academic examples are analyzed throughout the Chapter to benchmark the introduced improvements.

General VTCR applied to the shell theory

The general shell theory is illustrated in Chapter 2. The general VTCR applied to shell theory is formally similar to its shallow shell counterpart presented in Sections 3.1 and 3.2.

Simulation of impact response using a frequency approach However, in this case the coordinate system is curvilinear and the in-plane equilibrium Equation (2.12) for shells and the definitions of M and b present additional coupling terms. The discussions in Sections 3.1 and 3.2 are still formally valid in the shell theory framework given that the substitutions

x → α (4.1) y → β (4.2)
are performed and the quantities M and b are defined as in Equations (2.16) and (2.33). In previous works the weak form required to be adjusted for the chosen coordinate system. This was due to the different gradient and divergence formulations among different coordinate systems. However, when the Lamé parameters are constant we have

∇( ) = ([ , α , , β ]L) (4.3) ∇ • = [1, 1]∇( ) (4.4) L = 1 L α 0 0 1 L β (4.5)
where is a generic quantity and L is the Lamé matrix. In this case, the weak form remains the same and formally equal to Equation (3.21) for any curvilinear coordinate system. Conversely, shape functions are modified introducing a correction matrix into the ki formulation

ki = L -1 ki (4.6)
This upgrade keeps ∇2 ∇ 2 w1 in Equation (2.13) equal to k 4 . This characteristic stabilizes wavenumber computation increasing its accuracy similarly to the discussion in Section 3.1.2 about the norm of the direction vector introduced in the dispersion equation.

Various theoretical and performance improvements presented in Chapter 3 are included in the new TAPYROSS version with the introduced modifications for curvilinear coordinate systems. The following academic examples benchmark the code. In order to further improve performances, a new way to compute weak form integrals along straight boundaries is implemented and included in TAPYROSS. For the sake of classification, its described in Section 4.6. However, it is used in the following examples to test it on various academic cases.

Numerical test: a punctual force applied on a shallow shell

The present example studies a punctual force applied on a shallow shell. It is alike the test analyzed in Section 3.4.4. The reference is a quasi-analytic solution which is known in this particular case. It is quasi-analytic because it is embedded in the shallow shell approximations, the in-plane inertia is neglected, and the double infinite series in Equations (3.59) and (3.60) are stopped for numerical computation purposes.

This Section performs a triple comparison between VTCR, FEM, and the quasianalytic reference. Differently, performances are compared among the FEM, the TAPY-ROSS version with the shallow shell approximations, and its new version. At the end, the domain is sub-divided in four sub-domains to determine if the deterioration effect illustrated in Sections 3.7.1 and 3.7.2 is solved.

The multi-step VTCR convergence process is alike that one in Section 3.4.4 since theories are similar. In-plane waves are unnecessary in this case since the out-of-plane behavior is predominant. For the sake of clarity and simplicity, the error threshold is err V TCR ≤ 0.01 and In order to investigate the deterioration behavior highlighted in Sections 3.7.1 and 3.7.2, the domain is divided as in Figure 4.2 and the VTCR solution recomputed. The number of functions of the particular solution is kept equal to the previous case and applied only on the sub-domain where the point force is applied. As explained in Section 3.4.2, this choice produces a plateau in the convergence graphics of the other convergence steps after the error threshold is reached (err V TCR ≤ 0.01). The p-refinement of the propagative Simulation of impact response using a frequency approach Remark. The triple comparison in Figure 4.1 shows that the new VTCR is indistinguishable from the FEM solution. In fact, err FEM ≈ 10 -3 . This further sustains that the difference between the FEM and the quasi-analytic reference is related to the approximations of the latter instead of the FEM problems at mid-frequency. Performances in Table 4.2 point out that the improvement introduced in Section 4.6 produces a dramatic spare of computational costs. In particular, time consumption decreased of one order of magnitude compared with the old TAPYROSS version. Memory consumption was halved. This is a key factor when more complex problems such as the the HSS3+ test are involved.

Figure 4.4 demonstrates that the deterioration effect in Sections 3.7.1 and 3.7.2 is solved relaxing the shallow shell approximations. This result sustains our theory that it was due to virtual auto-equilibrated loads created by the angles among underlying planes.

Simulation of impact response using a frequency approach

Numerical test: a supported cylinder section

This problem is similar to the example provided in Section 3.4.5. However, shells are not shallow in this case. The study aims to determine the importance of in-plane waves in a problem where the out-of-plane behavior seems to be predominant. Moreover, the structure presents a triple joint to benchmark the coupling conditions. After, the domain is divided in five sub-domains to investigate if the deterioration effect highlighted in Sections 3.7.1 and 3.7.2 is solved even in this general case. At the end, performances are investigated on a frequency band f ∈ [2000; 4000] Hz. Geometry and boundary conditions are illustrated in Figure 4.5. Three sub-domains are connected by the same edge. The first two are cylinder parts while the last one is a plate. All boundaries are clamped but left edge where an out-of-plane oscillatory distributed load p = [1, 0, 0] N/m is applied.

This academic example studies the in-plane inertia effects when a cylinder is subject to out-of-plane loads. For this reason, a triple comparison among VTCR with and without in-plane inertia and a FEM reference is performed at fixed frequency. Material properties, frequency, and shell thicknesses are reported in Table 4.3.

According with Section 3.4.2, a initial guess VTCR with many rays is computed for convergence analysis. Error threshold is err V TCR ≤ 0.01 for every step. Convergence results are summarized in Table 4.4.

Since without in-plane inertia there are no in-plane propagative waves, convergence analysis for VTCR without in-plane inertia is equal to VTCR with in-plane inertia convergence arrested to second step. For the reasons reported in Section 3.4.2, the second Simulation of impact response using a frequency approach u = 0 u, n = 0 Table 4.4: Ray number and type for initial guess VTCR and converged VTCR with inplane inertia described in Section 3.1.1. converged VTCR is the result of the three-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Figure 4.6 where the error threshold is err V TCR ≤ 0.01. NS and SS stands for Normal and Shear Stresses respectively.

p = 1 N/m Ω 2 Ω 1 Ω 3 1 m 1 m
Simulation of impact response using a frequency approach and the third convergence analysis illustrate a plateau when the error threshold is reached.

The triple comparison among displacement magnitudes of VTCR with and without inplane inertia and a FEM reference is reported in Figure 4.7. The error between VTCR with in-plane waves is err FEM ≈ 8%. Conversely, VTCR without in-plane waves completely fails to find the correct solution.

Figure 4.9 compares the VTCR solutions with different domain sub-divisions. In particular, Figure 4.9a is divided in 3 sub-domains as in VTCR convergence processes for each frequency step are omitted for the sake of simplicity. However, the error threshold is err V TCR ≤ 0.01. DoFs cannot be directly compared because methods are intrinsically different. Nevertheless, for the sake of completeness, DoFs required over the frequency band are illustrated in Figure 4.10a.

Remark. The triple comparison in Figure 4.7 affirms that even in a problem where the out-of-plane behavior should be predominant, in-plane waves are crucial to converge to the correct solution. For this reason, in general they cannot be discarded a priori. At this Simulation of impact response using a frequency approach Ω 1 Simulation of impact response using a frequency approach [2000; 4000] Hz of the numerical example described in Section 4.1.2. VTCR and FEM are compared considering time and memory consumptions. DoFs required are reported for the sake of completeness, a direct comparison is unmeaningful due to theory differences.

Ω 2 Ω 3 Ω 4 Ω 5
Simulation of impact response using a frequency approach frequency the response is very sensitive to wave number variations. Since the FEM uses a limited DoF number, the discretized structure is usually stiffer than the continuous one.

Conversely, the VTCR shape functions solve precisely the equilibrium Equations (2.12) and (2.13). For this reason the FEM and VTCR solutions do not match perfectly. However, the difference is relatively small (err FEM ≈ 8%). Figure 4.9 further confirms that the degradation effects described in Sections 3.7.1 and 3.7.2 are solved relaxing the shallow shell approximations.

At mid-frequency FEM suffers of pollution effect. DoFs number and time and memory consumptions become prohibitive. Conversely, VTCR is unaffected as supported by Figure 4.10. Since VTCR remains stable, computational cost differences grow as frequency increases. Figure 4.10b and Figure 4.10c illustrate that time and memory consumption differences are already of some orders of magnitude and increase with frequency. Due to differences in memory and computational costs, VTCR could have been run on much less powerful machine (i.e. a laptop) than FEM (which needs a workstation). The performance enhancement in Section 4.6 exacerbates even further such effect.

Orthotropic and sandwich materials

Orthotropic materials are addressed correcting the k i formulation. Differently, sandwich shells are tackled using an homogenization approach.

Orthotropic materials

Since the material presents different wave speeds2 along α and β, the ki should change with the direction as illustrated in Figure 4.11.

[ [START_REF] Kovalevsky | A Trefftz approach for medium-frequency vibrations of orthotropic structures[END_REF]] suggested a multiplicative correction matrix O to tackle orthotropic materials in plates

ki = O ki (4.7)
where

O = 4 ρhω 2 D -1/4 x 0 0 D -1/4 y (4.8)
In this formulation O = O(ρhω 2 ). This dependency can lead to numerical difficulties during calculation of k i and âi because O ki is no more a unit vector. Such numerical trouble is further described in Section 3.1.2. The present approach proposes a slight ki (a) Wave distribution in an isotropic material.

O = 8 D α D β D -1/4 α 0 0 D -1/4 β (4.10)
In this case coefficients are dimensionless. Such modification generalizes the formulation increasing precision in computation of k i and âi . More details about determination of these two parameters are illustrated in Section 3.1.2. This theoretical improvement is benchmarked in Section 4.2.2.

Numerical test: a supported orthotropic cylinder section

This Section tests VTCR corrections for orthotropic materials explained in Section 4.2. A VTCR with in-plane inertia solution is compared with a FEM reference at fixed frequency. Geometry and boundary conditions are alike the example in Section 4.1.2. However, in this case the material is a typical aerospace composite. Table 4.5 reports its characteristics as well as frequency studied.

A FEM reference solution is calculated using the rule-of-thumb reported in Equation (3.48). VTCR convergence is studied by a three-step process described in Section 3.4.2. It is reported in Figure 4.12. Table 4.6 illustrates initial guess VTCR and converged ray numbers summarizing convergence analysis results. The error threshold is err V TCR ≤ 0.01.

Displacement magnitude comparison is depicted in Figure 4.13. Even if VTCR and FEM solutions are not equal the error based on kinetic energy is ≈ 7%.

Performances are omitted since results are almost identical to Figure 4.10 at 3700 Hz.

Simulation of impact response using a frequency approach f frequency 3700 Hz

E θ1 = E θ2 = E z3 Young moduli 125 GPa E y1 = E y2 = E y3 Young moduli 60 GPa G θy1 = G θy2 = G zy3 Shear moduli 18 GPa ν θy1 = ν θy2 = ν zy3 Poisson's ratios 0.3 ρ 1 = ρ 2 = ρ 3 densities 2000 kg/m 3 η θ1 = η θ2 = η z3 damping factors 0.001 η y1 = η y2 = η y3
damping factors 0.001 η G1 = η G2 = η G3 damping factors 0.001 Table 4.5: Orthotropic material properties and frequency examined of the numerical example described in Section 4.2. Simulation of impact response using a frequency approach

Remarks The small error confirms that orthotropic materials are correctly addressed and included in TAPYROSS. This improvement allows the study of a great variety of composite shell structures effectively extending VTCR applicability to common composite aerospace and automotive structures.

Sandwich materials

Sandwich materials are homogenized. This is completely legit due to the thin shell approximation already introduced in Section 2.1 for the Koiter shell theory. The present homogenization is akin to that one illustrated in [START_REF] Ventsel | Thin plates and shells: theory: analysis, and applications[END_REF]. The study is limited to symmetric skin-core-skin orthotropic sandwich shells since the HSS3+ test presents only this sandwich kind. Figure 4.14 depicts a generic symmetric skin-coreskin sandwich shell and its various characteristics. η α , η β , and η G for both skins and core are already incorporated in the various E and G reported in Equations (2.27) to (2.29). The quantities required in the Koiter thin shell theory are Simulation of impact response using a frequency approach h = 2h s + h c (4.11)

ρ = 2ρ s h s + ρ c h c h (4.12) B α = 2E αs h s 1 -ν αβs ν βαs + E αc h c 1 -ν αβc ν βαc (4.13) B β = 2E βs h s 1 -ν αβs ν βαs + E βc h c 1 -ν αβc ν βαc (4.14) B G = 2Gh s + Gh c (4.15) D α = 2E αs 3(1 -ν αβs ν βαs ) h c 2 + h s 3 - h c 2 3 + E αc h 3 c 12(1 -ν αβc ν βαc ) (4.16) D β = 2E βs 3(1 -ν αβs ν βαs ) h c 2 + h s 3 - h c 2 3 + E βc h 3 c 12(1 -ν αβc ν βαc ) (4.17) D G = 4G s 3 h c 2 + h s 3 - h c 2 3 + G c h 3 c 6 (4.18)
This technique is the standard method to address sandwich shells. Even the FEM elements that follow the Koiter shell theory use it. For this reason, numerical tests on this approximation are superfluous.

Conic structures and shells of variable thickness

This Section describes the additional approximations required to tackle conic structures and shells with variable thickness.

Conic structures

Conic structures pose major theory problems. Figure 4.15 illustrates a generic conic structure with its curvilinear coordinates α and β and its constant apex angle ξ.

Lamé parameters and curvature radii are

L α = β sin(ξ) (4.19) L β = 1 (4.20) R α = β tan(ξ) (4.21) R β = ∞ (4.22) Since L α = L α (β) and R α = R α (β)
, the theory exposed in Section 2.1 is no more valid. In this case the equilibrium equations are much more complex than Equations (2.12) Simulation of impact response using a frequency approach

ξ β r α L α R α (a) Generic conic structure. L αb L αt L αm R αb R αm R αt (b) Conic sector.
Figure 4.15: A generic conic structure and a conic sector. The apex angle ξ, the curvilinear coordinates α and β, the Lamé parameter L α , and the curvature radius R α are highlighted. L β and R β are omitted because trivial. The conic sector illustrates the bottom and top values of L α and R α . Section 4.3 further discusses conic structures. and (2.13) since various L β and R β derivatives appear. Consequently, even the developed VTCR theory is not applicable. Despite our efforts to develop a comprehensive and relatively simple VTCR theory for conic structures, the formulation remains too complex and still requires major approximations. Moreover, conic structures in the HSS3+ test are similar to cylinders since L β and R β do not considerably vary over cone sectors. For these reasons, we approximate them as constant

L α ≈ L αm = L αt + L αb 2 (4.23) R α ≈ R αm = R αt + R αb 2 (4.24)
where L αm (R αm ) are the average values of L α (R α ) being L αt (R αt ) and L αb (R αb ) its top and bottom values respectively as illustrated in Section 4.3.1.

Numerical test: a conic structure

The present Section tests the approximation introduced to address the conic structures. Geometry and the boundary conditions are illustrated in Figure 4.16. A cylinder and a cone sector are connected. The left boundary is subject to an edge transverse load of 1 N/m. The right boundary is clamped. The remaining edges are free. The structure is studied at different frequencies and with various cone apex angles θ a . They are reported in Table 4.7 with the material properties. As explained in Sections 3.4.2 and 3.5.1, VTCR and FEM solutions diverge as the frequency increases due to small theory differences. For this reason, the angle θ a = 0 is studied as reference of this effect.

The VTCR is compared with a FEM reference in the top view since it is the best suited to displays the surface. The FEM reference follows the rule-of-thumb discussed in Section 3.4.2. FEM DoFs are reported in Figure 4.22a. VTCR convergence processes are depicted in Figure 4.17 and their results are reported in Table 4.8. They should be performed for each frequency and for each θ a . However, there are no appreciable differences when θ a varies. Therefore, for the sake of simplicity, only one convergence study per frequency step is illustrated. The error indicator is set to err V TCR ≤ 0.01. Solution comparisons are reported in Figures 4.18 Remark. The present study highlights the effects of the apex angle θ a , the frequency, and their interaction. Let us focus on Figure 4.22. The VTCR solution degrades as the apex Simulation of impact response using a frequency approach 4.8.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach angle grows since it is directly related to the difference between the exact geometry and its approximation. When the structure does not present cones, the VTCR and FEM solutions diverge as the frequency increases. This further confirms the remarks in Sections 3.4.2 and 3.5.1. Conversely, when θ a = 0, the accuracy slightly increases with frequency. This effect is more accentuated when as θ a increases. This means that the cone approximation is more accurate as the frequency increases. At high frequency the inertia terms in Equations (2.12) and (2.13) become predominant on the coupling factors. For this reason, locally, the structural behavior tends toward the plate. The solutions in Figures 4.18 to 4.20 illustrate that the surface discontinuity between the cylinder and the cone sectors stops or drastically modify vibrations at low frequency. At higher frequencies this effect is less evident or negligible. Even the HSS3+ test exhibits this behavior where the most important picks are at mid frequency. Further information on the HSS3+ test are provided in Chapter 5. The reason is that at higher frequencies many modes are present and the energy flows throughout the whole structure unhindered. It is the same effect on which the energy based methods lie (see Section 1.1.2).

Performances in Figure 4.22 confirms once more that VTCR largely outperforms FEM at mid-high frequency. The performance improvement in Section 4.6 further increased this effect.

Shells of variable thickness

Shell sectors of variable thickness pose theory issues similar to conic structures. In fact, if h = h(α, β) then B and D derivatives appear in the equilibrium equations. For this reason, Equations (2.12) and (2.13) and the VTCR for thin shell structures are no more valid.

Simulation of impact response using a frequency approach Similarly to the approximation for conic structures, variable thickness is approximated as constant and equal to its integral average

h m h ≈ h m = L α L β h(α, β)dαdβ L α L β dαdβ (4.25)
4.3.4 Numerical test: a point force on a cylinder sector of variable thickness

This test investigates the approximation about shells of variable thickness in Section 4.3.

The LVA is the only HSS3+ component that present a variable thickness is the LVA. Thus, the study is restricted to shells of linear variable thickness for the sake of simplicity. The proposed vibrational problem is similar to the numerical test in Section 4.1.1 since FEM was accurate in this case. However, curvature and point force are modified to study a problem closer to the LVA geometry. In order to reduce the number of tests, we will consider various cases that present the same average thickness. This implies that the VTCR solution for all these cases is the same. For this reason, two different thickness function types are investigated along α and ŷ respectively

h = h(c h , α) = a h (c h )α + b h (c h ) (4.26) h = h(c h , y) = a h (c h )y + b h (c h ) (4.27)
where

a h = a h (c h ) = 2hc h l (2 + c h ) (4.28) b h = b h (c h ) = 2hc h 2 + c h (4.29)
h is the average thickness, l is the appropriate characteristic dimension (l α or l y ), and c h ∈ [0, 1] is a thickness parameter so that

c h = h(l ) -h(0) h(0) (4.30)
Varying only c h different thickness functions can be assessed. All these functions present the same average h and, therefore, the same VTCR solution. Comparing all these different vibrational problems with the VTCR solution the present test characterizes the error introduced approximating the thickness as constant as the thickness function changes. Figure 4.23 exposes geometry and boundary conditions. Table 4.9 reports frequency, material properties, and thickness parameters. Position and direction of the point force are provided in the global cartesian coordinate system highlighted in Figure 4.23b. VTCR convergence process is omitted because it is qualitatively similar to that Simulation of impact response using a frequency approach Remark. Figures 4.24 and 4.25 illustrate that the variable thickness drastically changes the vibrational result. However, Figure 4.26 affirms that the kinetic energies remain similar. When h = h(c h , α) err FEM ≤ 35%. Conversely, when h = h(c h , α) err FEM ≤ 30%.

g u = 0 v = 0 N x = 0 N y = 0 m = 0 w = 0 w = 0 m = 0 u = 0 N y = 0 w = 0 m = 0 v = 0 N x = 0 m = 0 w = 0 α (
Simulation of impact response using a frequency approach Figure 4.24: FEM solutions and the VTCR approximation of the vibrational problems described in Section 4.3.4 for c h = [0, 1 4 , 1 2 , 3 4 , 1] when the thickness is described by Equation (4.26).

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach R α radius of curvature along α 0.61 m R y radius of curvature along ŷ ∞ h average thickness 3 mm l α surface dimension along α 1.2215 rad Simulation of impact response using a frequency approach

l y surface dimension along ŷ 1 m c h thickness parameter [0, 1 4 , 1 2 , 3 4 ,

Moving loads

The present Section describes how moving loads are addressed in the VTCR. Section 4.4.1 presents some key properties of the frequency approaches. In particular, a moving load in the time domain is transformed in a non-uniform load (along the moving load path) in the frequency domain. Section 4.4.2 illustrate how these boundary loads are addressed in the VTCR.

Key properties of the frequency approaches

The focus is on the generic system illustrated in Figure 4.27. A box represents the system. A state-variable vector x c3 determines its characteristics. It is composed by the smallest variable subset required to represent the system at any given time. It receives inputs u and produces outputs y. Inputs, outputs, and state-variables are all signals. Signals are functions of time or frequency depending on the used approach. All these signals are related by first-order differential equations. For the sake of brevity and clarity, the discussion is limited to linear, time-invariant, causal systems since the HSS3+ specimen presents these characteristics. Linearity directly descents from equilibrium and constitutive relations. Causality derives from the fact that the system is real. The time-invariant hypothesis is not so trivial. For this reason, Section 5.2 is reserved to demonstrate this system property. In literature, there are two main ways to analyze a structure transient response: timeand frequency-domain approaches. The following paragraphs compare these methods providing some key notes about control theory.

Time-domain approaches study the system evolution in time. The standard state-space representation of a continuous time-invariant model is

x c , t (t) = A c x c (t) + B c u c (t) y c (t) = C c x c (t) + D c u c (t) (4.31)
where A c is the state matrix, B c is the input matrix, C c is the output matrix, and D c is the feedthrough (or feedforward) matrix. By definition, A c is always square. Conversely, other matrices are in general rectangular. In our specific case, the state-space representation of Equations (2.12) and (2.13) is

x c = [ u , u ] (4.32) u = u, t (4.33) u c = g (4.34) A c x c =   M -1 c ∇ • N -R (∇ • M) ∇ • (∇ • M) + R : N u   (4.35) B c = M -1 c 0 3 (4.36) C c = I 6 (4.37) D c = 0 6 (4.38) M c = ρhI 3 (4.39)
where I 3 is a 3 × 3 identity matrix, I 6 is a 6 × 6 identity matrix, 0 3 is a 3 × 3 null matrix, 0 6 is a 6 × 6 null matrix, and the product A c x c is presented instead of just the matrix A c for the sake of simplicity. Since Equations (4.37) and (4.38)

y c = x c = [ u , u ] (4.40)
After, the problem should be completed to consider boundary, corner, and coupling conditions. Typically, numerical implementations approximate a continuous differential problem into its discrete state-space representation

x c , t (t i+1 ) = A c x c (t i ) + B c u c (t i ) (4.41) y c (t i ) = C c x c (t i ) + D c u c (t i ) (4.42)
Such formulation is easily implementable. It can tackle non-linearities and timevariant systems. For these reasons it is very general. However, when fast-oscillating dynamics occur, the time-step should be very small to catch all the relevant frequencycomponents. Moreover, since the time is approximated with time-steps, the method presents a small error. Eventually, such discrepancy between continuous and discrete systems propagates throughout the integral Equation (4.41) increasing over time. This effect can deteriorate results after some time-steps. Typically, when initial conditions are uncertain, Monte Carlo simulations are required. This technique runs many simulations where random inputs (inside boundary values) are introduced and analyze results. Each simulation should be fully computed since results depend on initial conditions. When the system is linear, causal, and time-invariant, a useful strategy permits to drastically reduce Simulation of impact response using a frequency approach computational costs. Let us consider the particular case where a system is subject to a Dirac delta δ c (t) as input. Its output is h c (t). By definition, the convolution operator is

⊗ (t) = +∞ -∞ (τ) (t -τ)dτ = +∞ -∞ (t -τ) (τ)dτ (4.43)
The convolution between δ(t) and a generic input function u c (t) is

δ(t) ⊗ u c (t) = u c (t) (4.44)
since the Dirac delta is the convolution neutral element. For this reason, since the system is linear,

y c (t) = h c (t) ⊗ u c (t) (4.45)
where y c (t) is the output of u c (t). Hence, once h c (t) is computed, the solution is obtained just computing a convolution. However, the actual calculation of the convolution operator is complex. A c , B c , C c , and D c are unchanged since the system is time-invariant. The system in Equation (4.47) can be rearranged to erase state-variables and obtain just one equation

y cs (s) = G c u cs (s) (4.48)
where G c is the transfer function. It is

G c = C c (sI -A c ) -1 B c + D c (4.49)
where I is an identity matrix of the same dimension of A c . Even in this case the problem should be completed with boundary, corner, and coupling conditions. It becomes alike the reference problem exposed in Section 2.1 where s is implicitly restricted to s = ıω. Yet, the computation of the transfer function by mean of Equation (4.49) can be difficult. A fast Simulation of impact response using a frequency approach and simple way to compute a transfer function is to study the system response to a Dirac delta input δ(t) c . In fact, its Laplace transform is a step function cs (s), Equation (4.48) becomes

y cs (s) = G c cs (s) = G c (4.50)
since the step function is the multiplication neutral element. This is the equivalent to the convolution in the previous paragraph. The advantage is that in the frequency-domain the convolution operator is transformed in a multiplication. After the computation of Equation (4.48) the output y cs can be retransformed into the time-domain using the Laplace anti-transform

c = 1 2πı lim ω→∞ -∞ -∞ cs e -(σ+ıω)t dω (4.51)
where σ is a real number so that the integration path is in the convergence region of cs .

In this approach, differential equations are transformed in algebraic equations which are much simpler. Moreover, the transfer function does not change as the input function varies. Hence, input uncertainties can be taken into account with almost no additional computational efforts. However, these methods works only when the system is linear and causal. Time-variant systems can still be tackled but the problem resolution is much more complex and depends on the specific addressed problem. This thesis uses a frequency approach (based on the VTCR theory) since the HSS3+ test is linear, causal, and time-invariant and because the input signal (which are the explosion loads) are mostly uncertain due to the great amount of power released and its rapidity. These difficulties are treated in detail in Section 5.2.

Non-uniform conditions

The HSS3+ specimen is subject to shock loads along explosion paths. As explained in detail in Section 4.4, since the system is linear, time-invariant (this property is proved in Section 5.2), and causal, the frequency response of these detonation inputs is easily determined once the vibrational behavior of the structure subject to Dirac delta inputs is analyzed. Dirac delta inputs are transformed into non-uniform boundary conditions in the frequency domain. Hence, a method to tackle these conditions is required.

The present approach addresses every possible non-uniform boundary condition if some minor assumptions are met. Moreover, in case boundary conditions change among various simulations or are uncertain, all matrices can be reused skipping the matrix buildup process. Usually, it is the VTCR bottleneck. Therefore, this method drastically reduces computational efforts in these cases which are common in engineering. Be = (s) a generic non-uniform boundary condition where s is a curvilinear coordinate along the boundary and is one of the interesting quantities listed in Equation (3.29). In a similar way to Section 3.3, let us assume (s) is approximable with a 1D Fourier series where (e ık F s ) is a vector of harmonic functions of the Fourier series and F is its vector of complex scalar coefficients. As explained in Section 3.2, its relative linear form subvector l e of the weak variational formulation is in matrix form

l = ∂ Ω (e ık F s L [δu] H ) F ds = ∂ Ω e ık F s L [δu] H ds F (4.57)
where the dyadic product e ık F s L [δu] H produces a rectangular matrix n t × n F . n t is the number of test functions and n F is the number of harmonic functions of the Fourier series. In this way, when a different boundary condition is introduced, all matrices can be reused. In fact, only F changes. Typically, since the VTCR DoF number required to reach convergence is small, the resolution of Equation (3.28) is very cheap in terms of computational costs. Conversely, the matrix build-up operation can demand relevant computational efforts due to its numerous integrals required. Therefore, the present approach can dramatically reduce computational costs when boundary conditions change among different simulations. Moreover, this method can address every possible non-uniform boundary condition given that its Fourier series approximation exists. Boundary conditions in real test scenarios almost always meet this assumption. As minor advantage, the integrals required by this approach are equal to the integrals of the bilinear form since plane waves and harmonic functions of the Fourier approximation are formally alike. For this reason, only one computer routine is necessary to compute all the required integrals drastically saving code-typing and debug time. Remark. Figure 4.29 demonstrates that the theory introduced in Section 4.4.2 is correctly implemented. In Figure 4.30 the curves does not match perfectly since the VTCR is a Trefftz method which approximates boundary conditions to privilege equilibrium equations. u always presents a better approximation since in-plane waves vibrate less than Simulation of impact response using a frequency approach out-of-plane ones due to the difference between B and D. The solution is along the x axis because the imaginary part is an odd function. This tilts the waves introducing an odd delay along the boundary. At low-frequency (800Hz) this effect is less appreciable due to the long wave length. In particular, the u seems to be symmetric at 800 and 1800 Hz. Yet, at 3400 Hz is becomes noticeable.

Concentrated masses and springs

In the HSS3+ test numerous concentrated masses are present along boundaries ρ b = ρ b (s) and on corners ρ C due to shell sector joints. On one hand, these masses are approximated as constant along boundaries. On the other hand, joint stiffnesses are considered infinite. For this reason, only constant concentrated masses are present in the HSS3+ test. Yet, since concentrated masses and springs behave alike, concentrated springs are included in this discussion.

Be a stress constraint among b α , b β , q, m along a boundary. As described in Equations ( 2

µ b = µ b (ω) = ρ b ω 2 + κ b (4.61)
The displacement e along a boundary is unknown a priori since it is a part of the VTCR solution. In particular, it is the VTCR solution computed along the boundary. In the weak form, Equation (4.60) is multiplied by every test function L [δu i ] and integrated along the boundary. Since Equation (4.60) is a stress constraint, L [ ] is an operator that computes the displacement quantity along the boundary. The inertia-spring term µ b e and e belong to B ee which is a sub-matrix of the main block diagonal. The terms Simulation of impact response using a frequency approach e Γe compose off-diagonal blocks and take into account coupling conditions. e is a part of the linear form. Weak form integrals related to µ b are

µ b ∂ µ b Ω e L [δu i ] H L [u j ]ds (4.62)
where ∂ µ b Ω e is the portion of ∂Ω e where µ b = 0 and µ b is extracted from the integral because constant. It is interesting to highlight that in this case the displacement quantity is multiplied by itself instead of its dual quantity . It is correct because µ b e is a stress quantity. Be B eeµ b the part related to µ b of B ee . Since the VTCR is a Galerkin method, B eeµ b is, in matrix form,

B eeµ b = µ b ∂ µ b Ω e L [u e ] L [u e ]
H ds (4.63

)
where there is a dyadic product of two vector of functions L [u e ] where shape function amplitudes are a j = 1. B eeµ b is a symmetric, positive-definite bilinear form by definition. This ensures uniqueness and existence properties of the weak form even in this particular case. Concentrated masses and springs on corners follow the same discussion except that, at the end, integrals are substituted with sums on corners as in Equation (3.33) with respect of Equation (3.32). Section 4.5.1 benchmarks this VTCR improvement.

Numerical test: a plate with a linear distributed mass

Geometry and boundary conditions are illustrated in Figure 4.31. The left edge of a plate is subject to an out-of-plane edge load p, top and bottom boundaries are completely fixed (u = 0, w, n), and the right side is free and presents a linear distributed mass m l along the boundary. Table 4.11 reports material properties and the frequency studied. In this case no analytic solutions are known. For this reason, a FEM solution is used as reference. Since the focus is on the mid-frequency, the maximum FEM mesh size is subject to Equation (3.48). The multi-step convergence procedure of the VTCR described in Section 3.4.2 is depicted in Figure 4.32 and its results are summarized in Table 4.12. Propagative in-plane waves are unnecessary since the problem is completely out-of-plane. The error threshold is err V TCR ≤ 0.01. Results are compared in Figure 4.33. The difference between the two solutions is err FEM ≈ 6%. VTCR and ABAQUS R performances are reported in Table 4.13.

Remark. The VTCR solution almost perfectly matches the FEM reference using much less computational resources. This validates the theoretical improvement introduced in Section 4.5. Moreover, the performance comparison in Table 4.13 confirms once more the effectiveness of the VTCR method in mid-frequency.

Simulation of impact response using a frequency approach Table 4.13: DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 4.5.1.

Simulation of impact response using a frequency approach

A fast integral computation technique

Weak form integrals along straight boundaries4 can be computed analytically reducing computational costs. The present Section analyzes in detail these integrals proposing a new effective algorithm to compute and store these integrals. Such improvement reduces computational costs of some order of magnitude as is illustrated in the previous numerical tests in Sections 4.1.1,4.1.2,4.2.2,4.3.2,4.3.4,4.4.3 and 4.5.1. For the sake of clarity, let us recall the standard weak form integral

∂ Ω e ∩∂ Ω e L [δu i ] H L δu j ds (4.64)
where is the imposed condition among u, b α , v, b β , w, q, w, n , m, is its dual quantity, δu i is a test function, u j is a shape function (where a j = 1), and L [ ] and L [ ] are two operators that extract the quantities and respectively from the function introduced. They depend on â, k, and n. In general n = n(s) where s is a curvilinear coordinate along the boundary. Let us focus on a straight boundary ∂ i b Ω e of Ω e shared with Ω e where i b is the boundary index. In this case n is constant along the boundary. Therefore L [δu i ] = c i e ık i r reli (4.65)

L δu j = c j e ık j r rel j (4.66)

where c i and c j are two constant scalars. For this reason, weak form integrals can be simplified as

c H i c j ∂ i b Ω e
e ı(k j r rel j -k H i r reli) ds (4.67)

The linear integral ∂ i b Ω e e ı(k j r rel j -k H i r reli) ds can be computed analytically since the boundary is straight. Be p 1 and p 2 the boundary extremes. Be p 1rele and p 2rele and p 1rele and p 2rele their positions relative to the sub-domains Ω e and Ω e respectively. Be where p 21 is the curvilinear distance and L is either L α or L β depending on which principal direction lies the boundary (in case of a Cartesian coordinate system L is unnecessary).

Therefore, Equation (4.67) becomes

L p 21 c H i c j e ı(k j p 1rele -k H i p 1rele ) 1 0 e ı(k j -k H i )p21ς dς = =L p 21 c H i c j e ı(k j p 1rele -k H i p 1rele )    e ı ( k j -k H i ) p 21 -1 ı(k j -k H i )p21 if k j -k H i = 0 1 if k j -k H i = 0 (4.71)
These integrals compose rectangular n t × n s matrices Q e e∂ Ω e ∩∂ Ω e relative to the boundary ∂ Ω e ∩ ∂ Ω e . In particular, they are the same for any boundary condition. Conversely, the term c i depends on the boundary, the particular condition , and the function δu i but not u j and vice versa. Hence, these quantities can be stored in rectangular n s × 8 matrices C ∂ i b Ω e as depicted in Figure 4.34. They are relative to the sub-domain Ω e and its boundary ∂ i b Ω e being i b a boundary index. where C ∂Ω e and C ∂Ω e are matrix columns relative to and respectively. In this way, during B construction, sub-matrices can be built on-demand just by multiplying matrices. This algorithm provides three major advantages:

u 1 u 2 u j u n s . . . . . . u v w w, n b α b β q m C ∂ i b Ω e
• if boundary conditions are modified, matrices can be reused providing great flexibility to the method;

• time consumption is reduced of some order of magnitude since integrals should be computed once for all boundary conditions;

Simulation of impact response using a frequency approach

• even memory storage consumption is greatly reduced since the elements of C ∂Ω e , C ∂Ω e , and Q e e∂ Ω e ∩∂ Ω e are much less than the elements of a 3D matrix n t × n s × 8 which stores the integrals for any possible boundary condition5 .

Remark. For the sake of simplicity, no specific academic tests are performed to test the present technique. Conversely, it was used in all the previous numerical examples proposed throughout this Chapter. Performances increased of at least one order of magnitude.

In particular, the comparison between Section 3.4.4 and Section 4.1.1 in Table 4.2 highlighted that the VTCR computational time dropped from 4 to 0.3 seconds in a problem where the standard FEM reference required 18 seconds. At the same time the VTCR accuracy increased. These results further confirm that the VTCR greatly outperform the FEM at mid-frequency.

Conclusions

In this chapter the VTCR and its improvements introduced in Chapter 3 are extended to the general shell theory. The new TAPYROSS version is benchmarked on two relevant numerical examples in Sections 4.1.1 and 4.1.2. In particular, results in Section 4.1.1 confirm that the new TAPYROSS version does not present the major issues illustrated in Sections 3.7.1 and 3.7.2. This supports the explanation provided at the end of Chapter 3. In fact, since the VTCR for shells uses a curvilinear coordinate system to describe the problem geometry, the angle between adjacent cylinder sectors is null. In this way no virtual auto-equilibrated loads arise and the solution is unaffected by the particular domain sub-division.

The remaining Sections present the other improvements required to study the HSS3+ test. Section 4.2 introduces a corrective matrix O in the k formulation to address orthotropic and sandwich materials. Sandwich structures are homogenized. The correction for orthotropic materials is similar to the solution illustrated in [START_REF] Kovalevsky | A Trefftz approach for medium-frequency vibrations of orthotropic structures[END_REF]. However, the present dimensionless formulation simplifies the computation of the wavenumber k increasing the robustness of the method. Section 4.3 discusses the approximations required to tackle conic structures and shells of variable thickness. Numerical tests in Sections 4.3.2 and 4.3.4 illustrate that, in some cases, the approximated results can be quite different from the real solution. However, error indicator based on the kinetic energy is always err FEM ≤ 35%. Section 4.4 describes that moving loads (i.e. a pyrotechnic explosion) are transformed in non-uniform boundary loads in the frequency domain. These loads are taken into account in the VTCR using a flexible and inexpensive method. Section 4.5 expand the VTCR formulation to take into account concentrated masses and springs. These theoretical improvements does not noticeably increase computational costs and permits to reuse matrices when the concentrated masses or springs change or when a non-uniform boundary condition is unknown. This is an useful property when boundary conditions are unknown (i.e. a pyrotechnic explosion). Section 4.6 introduces a fast computation algorithm for the polygonal sub-domains. When applicable, it drastically reduces VTCR computational costs. This is confirmed by the numerous numerical tests in which it was used.

After all these tests and enhancements, TAPYROSS is ready to study the HSS3+ test. Chapter 5 compares VTCR and real data results to determine damping coefficients and explosion input and to finally validate the computer program on a real, complex industrial test.

Chapter 5

HSS3+ test

The present Chapter is devoted to the HSS3+ test. Real data and simulations are compared to better understand the propagation of pyrotechnic shocks in launchers. The aims of the present analysis are threefold:

• to model the explosion input,

• to determine the average damping coefficient,

• to validate TAPYROSS on a complex industrial test case.

Sections 5.2 and 5.2.2 discuss some theoretic key points about explosions. Section 5.3 qualitatively analyzes the real data. Section 5.3.4 characterizes the explosion input. Section 5.3.5 studies the damping coefficient. Section 5.4 compares VTCR and real data frequency responses analyzing results.

General overview of the HSS3+ test

The present Section illustrates the HSS3+ test. It is restricted to a qualitative description of the structure and the pyrotechnic explosions due to copyright. However, the interested reader can see [Perez 2011a, Bézier 2013]. The HSS3+ test investigates the Ariane 5 vibrational response during fairing detachment. Typically, such operation happens at high altitude where air-structure interactions Simulation of impact response using a frequency approach are negligible. The test is performed on the ground to greatly increase accuracy and simplify the study since there are no engine-induced vibrations and the specimen is not flying. Air-structure interactions are assumed negligible even in the HSS3+ test for the sake of simplicity. Figures 5.1a and 5.1b depict the test. For the sake of clarity, Figure 5.1c shows the HSS3+ cross section and the sensor locations. The PLF cylinder represents the fairing. It is pre-divided in two parts to allow fairing detachment. When it is severed from the structure by pyrotechnic explosions the PLF cylinder parts fall and the specimen vibrates.

The following Paragraphs describe the HSS3+ structure stressing out some key points.

The composite cylinder was developed for the Ariane 5 fairing QM2 validation tests. It has been modified multiple times over the years. The present study focuses on the original version. It is a complete cylinder constituted of sandwich composite shell panels.

It is pre-divided in two parts to permit faring detachment. Before the explosions these sections are connected to the main structure. After that these parts are severed from the specimen they fall. The ceiling hanging prevents them to touch the ground to not produce unintended vibrations.

The ACY 5400 is another sandwich composite shell cylinder. The top boundary is connected to the SYLDA and the flange HSS3+ creating a triple-joint.

The flange HSS3+ is a small metallic flange divided in four sections (90 o each). The explosive is positioned along the top boundary between the PLF cylinder and the flange itself.

The PLF cylinder is a metallic cylinder that simulates the fairing.

The SYLDA is a sandwich composite shell structure composed of two cones an a cylinder. It sustains satellite 1 on the top of the ACU and encompasses satellite 2. Conversely from Ariane 5, there is no second satellite inside the SYLDA structure in the HSS3+ test to simplify the structure.

The ACU is composed of the LVA and the PAF. It is a conic interface that connects the SYLDA with the payload.

The LVA is a sandwich composite shell cone with variable thickness.

The PAF is a metallic structure that connects the LVA with the payload.

The payload is simulated by a complex metallic frame structure with standard weight and first resonance peak due to the variability of the payload characteristics.

Simulation of impact response using a frequency approach Figure 5.1: HSS3+ structure before and after the pyrotechnic explosions, its cross section with sensor locations, and a particular of the explosion paths.

Simulation of impact response using a frequency approach Joints are simplified as perfect connections to greatly reduce computational costs. In order to take into account the effects of the added material required to create the joint, concentrated masses and radial displacement constraints are added.

The explosion paths are illustrated in Figure 5.1d. The detonations are symmetric around the rocket to not interfere with its trajectory during launch. The explosions are much faster than the shell speed of sound in the interesting frequency range. This characteristic is further discussed in Section 5.2.

Sensors are positioned all around the rocket in asymmetric positions. They are concentrated around the triple-joint in the explosion proximity and near the payload as depicted in Figure 5.1c.

Characteristics of the pyrotechnic shocks

The present Section discusses some key properties of the pyrotechnic detonations.

Consequences of a supersonic explosion speed

Let us focus on the fairing detachment scheme illustrated in Figure 5.2. The horizontal section in Figure 5.3a illustrates that four symmetric pyrotechnic explosions severe it from the rocket following the enumerated explosion paths. Detonations are ignited in pairs in A and B, propagate along the arrows, and terminate in C and D. Figure 5.3b depicts an explosion at a generic instant t. In general, the problem geometry changes as the explosion evolves. This dramatically complicate the input formulation. However, the extreme rapidity of the explosion process permits to avoid such problem. By definition, the speed of sound is the group velocity c g Simulation of impact response using a frequency approach Shell present very similar values. c gip is constant with frequency while c goop increases. Table 5.1 reports the highest speeds of sound for in-plane and out-of-plane waves along α and y at 15000 Hz where α and y correspond to the curvilinear coordinates α and β of a generic cone introduced in Section 4.3.1. Table 5.1: The highest speeds of sound for in-plane and out-of-plane waves along α and y at 15000 Hz.

Pyrotechnic explosions progress at v pyro ≈ 7100 m/s. In the HSS3+ test they are faster than c gip and c goop in the interesting frequency range. Thus, explosions present two distinct Mach cones for in-plane and out-of-plane waves. By definition, their relative Mach numbers M ip , M oop and cone half-angles α Mip , α Moop are

M ip = v pyro c ip (5.4) M oop = v pyro c oop (5.5) α Mip = arcsin(M -1 ip ) (5.6) α Moop = arcsin(M -1 oop ) (5.7)
where v pyro is the explosion speed. Hereafter the subscript pyro refers to the pyrotechnic explosions. Let us focus on Figure 5.3b. Explosion effects are relegated into the Mach cone. On one hand, along the future explosion path (outside the Mach cone) there are neither loads nor displacements. On the other hand, along the past explosion path (inside the Mach cone) the structure is already decoupled. Thus, the problem in Figure 5.2 is equivalent to that one illustrated in Figure 5.4 for what concerns explosion vibrations.

In this case the rocket has no fairing and presents four explosion loads along the upper boundary of the flange HSS3+. Explosion loads are assumed

b pyro = b pyro (α ± v pyro R α t) (5.8) q pyro = q pyro (α ± v pyro R α t) (5.9) m pyro = m pyro (α ± v pyro R α t) (5.10)
where α is a curvilinear coordinate and the sign ± depends on the particular explosion considered. For the sake of clarity, we recall the definitions of the boundary loads {b, q, m} described in Section 2.2. Let us consider Figure 5.5a. The boundary loads are Simulation of impact response using a frequency approach Pyrotechnic Explosions (a) HSS3+ before the explosions.

(b) HSS3+ after the explosions. where the subscript e refers to the boundary b of the sub-domain Ω e , te is a tangent unit vector, e Γe is the index relative to an other sub-domain that shares with Ω e the boundary Γ e , and n Γe is their total number (Ω e excluded) as is shown in Figure 5.5b. Section 4.4 explains that the frequency response y cs of any possible input can be trivially calculated using

y cs = G c x cs (5.14)
where x cs is the Fourier transform of the input load and G c is the system transfer function.

It is numerically equal to the system response to a Dirac delta input. Thus, the vibrational problem produced by any possible input condition can be determined computing only the particular case of a Dirac delta load (in the time domain). In order to study this case the load should be transformed in the frequency domain. The Fourier transform of a Dirac delta load that moves at the speed of the explosion along the detonation path is

δ cω = exp(± ıωR α α v pyro ) (5.15)
where the subscript cω refers to the Fourier transform. The functions {b pyro , q pyro , m pyro } remain unknown. However, in literature (i.e. [Klein 2015]) they are supposed to be similar to triangle functions ∧ or Dirac deltas (which are degenerated triangle functions) due to the extreme rapidity of the explosion process. The most general triangle function is

∧ =          0 τ < 0 h ∧ τ a ∧ 0 ≤ τ < a ∧ h ∧ (a ∧ +b ∧ -τ) b ∧ a ∧ ≤ τ < a ∧ + b ∧ 0 a ∧ + b ∧ ≤ τ (5.16) ∧ cω = h ∧ a ∧ b ∧ ω 2 [(a ∧ + b ∧ ) exp(ıa ∧ ω) -a ∧ exp(ı(a ∧ + b ∧ )ω) -b ∧ ]
(5.17 where

h ∧ 2 (a ∧ + b ∧ ) is the triangle area. Moreover, lim ω→0 ∧ cω = h ∧ 2 (a ∧ + b ∧ ) (5.19)
Equation (5.19) proves that at low-frequency (compared to a ∧ and b ∧ ) the explosion input can be approximated as a constant gain. Conversely, at higher frequency the shape of the explosion input function could be crucial. Since {a ∧ , b ∧ , h ∧ } are unknown only the comparisons in Section 5.4 can determine which is the prevalent behavior in the HSS3+ test at the interesting frequency range. Simulation of impact response using a frequency approach equivalent to the time required by the explosion to reach the point B. Thus,

Delay computation

d s = a s + b s tan -1 (α Mip ) v pyro (5.21)
When the waves pass multiple materials Equation (5.21) becomes

d s = a s v pyro + v -1 pyro ∑ i m b i m tan -1 (α Mipi m ) (5.22)
where i m is the structure component index, b i m is its height, and α Mipi m is its cone halfangle.

Real data and first analysis

The present Section discuss some important characteristics of the real data. Due to structure and load symmetries, the solution has two planes of symmetry. Thus, sensors are positioned all around the rocket in non-symmetric positions and at various altitudes as illustrated in Figure 5.8. The used sensors are very fast accelerometers with a sampling frequency of 125 kHz. They are positioned in pairs to capture axial and radial accelerations1 . The studied time window is of the order of milliseconds. Remark. The frequency response in Figure 5.9 is dominated by in-plane vibrations around 12 kHz. A smaller vibration cluster is located around 4 kHz. Conversely, in Figure 5.10 the out-of-plane waves around 4 kHz lead the response and smaller vibration picks are located around 12 kHz. This suggests that the vibration cluster at 12 kHz can be related to the explosion input while the picks at 4 kHz can be associated to the structure vibrational response. This hypothesis is further investigated in the following Sections.

Major characteristics of the real data

The time-frequency plots

A time-frequency plot shows how the frequency response evolves as the time passes. A time window2 is translated on a time-band (the ordinate) passing by many time-steps. For each time-step a frequency response is computed. The envelop of the frequency responses Remark. Figure 5.11 highlights the behavior detected in Section 5.3.1. The frequency cluster at 12 kHz is heavily damped and disappear after the first instants. These frequencies are crucial near the explosion. However, they are less important near the payload. Conversely, the peaks around 4 kHz are crucial near the satellite since they are slightly damped. These results confirm the remarks discussed in Section 5.3.1. Section 5.3.5 further discusses the damping.

Further characterizations the vibration clusters

In order to distinguish the explosion input from the structure modal behavior, two different time windows are investigated. They are ten times smaller that the total signal time length. Figure 5.14: Time and frequency domain responses of an accelerometer near the explosion on a small time window at the end of the signal.

Simulation of impact response using a frequency approach However, a particularized study of these frequencies is not performed here. However, they are not essential for designers since they are heavily damped. After a transient, these shock waves disappear transferring their energy to the structure modes at 4 kHz. These modes are mainly out-of-plane. Since an explosion is a rapid phenomenon, its Fourier transform at low-and mid-frequency should be similar to a constant gain. Thus, the VTCR Dirac delta response is sufficient to predict the frequency response around 4 kHz up to a multiplicative constant. For this reason, VTCR simulations in Section 5.4 focus on this frequency range to validate the computer program on a real industrial test.

Input analysis

Section 5.3 concludes that the explosion input is characterized by a vibration cluster around 12 kHz. Conversely, peaks around 4 kHz are related to the structure modal behavior. Section 5.3.4 depicts the qualitative behavior of the absolute value of the input function in the frequency domain. Since Equation (4.50), the function is similar to a constant gain at low-and mid-frequency. The VTCR frequency response multiplied by the correct gain is sufficient to determine its frequency response at 4 kHz.

Simulation of impact response using a frequency approach Around 12 kHz the input presents a vibration cluster which leads the frequency response. Usually, at high-frequency the damping coefficient η slightly decreases with frequency [START_REF] Troclet | Prediction of the vibroacoustic response of the Ariane 4 fairing and the equipment bay[END_REF]]. However, the complete structural damping is ηk which greatly increases as the frequency rises. Thus, these frequencies are irrelevant for engineering purposes. Moreover, specifications in [START_REF] Perez | Mechanical Environment. In Ariane 5 user's Man[END_REF]] assert that vibrations after 10 kHz can be neglected, confirming our statement. For this reason, they are not considered in Section 5.4. Damping is further discussed in Section 5.3.5.

Modelling of the damping coefficient

The determination of the damping coefficient is indeed a complex task. Klein in [Klein 2015] studies the HSS3+ test comparing real data results with the simulations performed Simulation of impact response using a frequency approach using the RADIOSS method which is a FEM approach developed to study shocks. After a detailed evaluation of the time-frequency graphics (i.e. where K RADIOSS is the stiffness matrix, M RADIOSS is the matrix of inertia, and α RADIOSS and β RADIOSS are two coefficients that are determined fitting real data results with the model. α RADIOSS addresses the structural damping while β RADIOSS considers the additional damping introduced by the structure-air interaction. In the end, he affirms that α RADIOSS is preponderant. For this reason and since the VTCR can address the structural damping only, we assume that the damping due to the structure-air interaction is negligible.

Troclet in [START_REF] Troclet | Prediction of the vibroacoustic response of the Ariane 4 fairing and the equipment bay[END_REF]] affirms that at low-frequency the structural damping is constant. Conversely, at mid-and high-frequency it varies with frequency using the following semi-empiric formula For this reason, the total damping ηk increases as the frequency increases for inplane waves and remains almost stable or increases for out-of-plane waves. In particular, explosion waves are heavily damped at high-frequency. For this reason, the VTCR study in Section 5.4 focuses on the peaks at 4 kHz neglecting the vibration cluster at 12 kHz. Since there is no established literature on the cutting edge materials used in the Ariane 5, these parameters are uncertain. For this reason, this Section presents the algorithm used in Section 5.4 to extrapolate them from real data.

Low-frequency techniques provide a damping coefficient that is constant with frequency focusing on the time response. However, the interesting frequencies (around 4 kHz) lie in the mid-frequency range and its frequency dependency could be crucial [START_REF] Troclet | Prediction of the vibroacoustic response of the Ariane 4 fairing and the equipment bay[END_REF].

In order to illustrate the method, let us focus on a sensor near the satellite. Figure 5.17a is its total response in the time domain. where t is the time distance between the two time window start points. This process is repeated for many sensors increase the reliability of the method. Finally, a η and b η are extrapolated from η.

VTCR and real data comparisons

The present Section compares the VTCR solution with real data. As mentioned before, the interesting frequency range is around 4 kHz. The frequency cluster at 12 kHz can be neglected since these vibrations are heavily damped in time and space. For this reason, the chosen frequency range is [1500;5000] Hz. Acceleration responses in the frequency domain are extracted from the VTCR solution and compared with accelerometer data. The considered time window on which the frequency response of the accelerometers is computed is located where the modal behavior is well established. Matched explosion input and damping coefficients are reported in Table 5.2. The matching processes are illustrated in Sections 5.3.4 and 5.3.5. For the sake of clarity, presented VTCR results are already calculated using the correct values. Table 5.4: Ray number and type for initial guess VTCR and converged VTCR at the highest interesting frequency (5000 Hz) for the studied vibrational problem is illustrated in Section 5.4. The error threshold for the convergence process is err V TCR ≤ 0.01.

Figure 5.20 compare in-plane and out-of-plane frequency responses for sensors near the explosion, the satellite, and in the middle respectively. For the sake of completeness, three graphic types of the same results are illustrated:

• time-frequency plots,

• semilog comparisons,

• Shock Response Spectrums (SRS).

The time-frequency plots of the interesting frequency range

The present Section focus on the left-high quarter of the time-frequency plots in Simulation of impact response using a frequency approach

Comparison of the graphics of the frequency responses

These graphics compare VTCR and real data frequency responses in normal and semilog scale. Let us define a norm a k of the acceleration signal that is related to the kinetic energy

a k = f a 2 2ω 2 d f (5.32)
where the integral is computed on a frequency band f . The difference between VTCR and real data is evaluated using an index err f based on it Remark. The errors reported in Table 5.6 do not seem related to the frequency. Some produce a bump, while others increase, decrease, or remain stable on the whole frequency band. Conversely, the error behavior is associated to the particular frequency response of the measured point. This sensitivity to the local behavior is typical of the mid-high frequency responses where small model uncertainties can drastically change structure behavior. Thus, in order to further improve the VTCR, a technique to take into account uncertainties should be incorporated into the method.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach

The comparisons of the Shock Response Spectra

A Shock Response Spectrum illustrates how a Single Degree Of Freedom (SDOF) system (like a mass on a spring) responds to that input. It reports the peak acceleration of a SDOF of variable resonance frequency. Peak acceleration response is represented on the ordinate. The natural frequencies of the variable SDOF oscillator is reported on the abscissa. Broadly, it highlights how much a SDOF would be exited by the input function if it had a natural resonance frequency ω SRS . In signal theory it reports the maximum overshoot of the signal (which is the acceleration) modulated by a variable low-pass filter. Figure 5.23 compares the SRS for VTCR and real data results. Figure 5.24 reports the very same results in semilog scale to highlight magnitude differences.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Figure 5.24: SRS comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle. A semilog scale is used to highlight magnitude differences.

Simulation of impact response using a frequency approach

Remark. VTCR and real data Shock Response Spectra present similar behaviors. In particular, Figure 5.24 highlights that excitation levels are similar. Yet, the frequency response seems to be slightly shifted in frequency and peak differences are still noticeable.

Many reasons can explain these relatively small discrepancies:

• Model uncertainties can play a crucial role at mid-frequency. Small differences between the model and the real test can drastically change the frequency responses.

Since there is just one real full-scale test, we cannot statistically compare results to increase reliability.

• The model was simplified. Sensors, small equipment, joints were not taken into account.

• The VTCR approximations required to address conic structures and shells with variable thickness are detrimental for the solution accuracy as illustrated in Section 4.3.

5.4.3.1 SRS sensitivity study of the VTCR simulations to the approximations introduced to address variable-thickness shells

In order to investigate the SRS sensitivity to the approximations introduced to tackle conic structures and shells with variable thickness, a new, slightly different VTCR simulation is performed. Every parameter remains the same except the approximation in Equation (4.25) for variable-thickness shells that is modified in the following way h ≈ h m = 0.9 L α L β h(α, β)dαdβ L α L β dαdβ (5.42)

The SRS comparisons between the two VTCR simulations are illustrated in Figures 5.25 and 5.26. Remark. Figures 5.25 and 5.26 indicate that a slight change in the approximation formula can relevantly modify the VTCR SRS. Peaks can shift of some hundreds of Hertz and their amplitude can double. In particular, Figures 5.26a to 5.26d present peak (or peak clusters) that increased or doubled while Figures 5.26c to 5.26f show peaks that shifted of hundreds of Hz. Sensors near the cone of variable thickness tends to present a relevant frequency shift. Conversely, sensors far away from the approximated cone do not show this behavior. This confirms that even a slight modification of the average thickness introduced into the VTCR model can significantly change the VTCR SRS. Thus, the values introduced to approximate these complex shell sectors should be carefully evaluated since the SRS results are very sensitive to thickness changes. In particular, this suggests that the simple averaging algorithm proposed in Equation (4.25) could be non-optimal.

Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach

Conclusions

Results in Section 5.4 illustrate that VTCR results somewhat match real data. Often, real data peaks can be detected on the VTCR response even if slightly shifted in frequency and the average excitation level is similar. However, differences are still relevant. The following reasons can explain these discrepancies:

• The VTCR approximations required to address conic structures and shells with variable thickness are detrimental for the solution accuracy as illustrated in Section 4.3. In particular, Section 5.4.3.1 further investigates such eventuality. It confirms the general conclusions of Section 4.3. The frequency response (even the SRS) are highly sensitive to the approximations introduced to address conic structures and shells of variable thickness. This indicates that the solution accuracy can be greatly improved changing the averaging methods proposed in Section 4.3.

• Model uncertainties can play a crucial role at mid-frequency. For example, the measuring equipment was neglected and the particular geometry of the joints was simplified. Section 5.4.2 highlights that even small differences between the model and the real test can drastically change the frequency responses. Therefore, the development of a method to address uncertainties can drastically increase the VTCR accuracy.

Simulation of impact response using a frequency approach

Conclusions

Recently, numerical prediction techniques become a fundamental tool for aerospace and automotive industries. Virtual testing dramatically reduces costs since it reduces the number of required real prototype tests simplifying and speeding-up the design process. This is especially true for aerospace industries where experiments are extremely costly. Virtual testing of shock propagation is still an open question. Usually, the relevant resonance peaks lie in the low-and mid-frequency. The mid-frequency range presents a mixed behavior between low-and high-frequency. The local response is still required and the system is more sensible to uncertainties that at low-frequency. On one hand, highfrequency methods does not provide the local response. On the other hand, computational costs of the standard low-frequency methods are prohibitive. For these reasons, standard commercial codes are not fully suitable to analyze these problems. Moreover, shock loads are unknown and structures are often complex composite shells.

The present thesis developed a reliable computer program (called TAPYROSS) capable of simulating these complex industrial problems to reduce the number of real tests required. It is based on the VTCR which is a Trefftz method specifically optimized to tackle mid-frequency problems. It is a frequency approach where the response is computed on some frequency steps. A key feature of the method is that the structure response should be computed only once for any possible shock input dramatically reducing computational costs. Many new theoretical and performance improvements were introduced to tackle complex industrial problems. The theoretical improvements introduced in Chapters 3 and 4 address:

• in-plane inertia (Section 3.1),

• particular solutions for general surface loads (Section 3.3),

• orthotropic and sandwich materials (Section 4.2),

• conic structures and shells of variable thickness (Section 4.3 ),

• moving loads (Section 4.4),

• concentrated masses and springs (Section 4.5),

The performance enhancements illustrated in Chapters 3 and 4 are:

• a quasi-symmetric ray distribution algorithm (Section 3.6), Simulation of impact response using a frequency approach

• iterative solvers (Section 3.7),

• a fast integral computation technique (Section 4.6).

In the end, TAPYROSS was tested on the results of the HSS3+ test. It is a real fullscale ground test conducted by the Centre National d' Études Spatiales (CNES) and Airbus Defence & Space (Airbus DS) to study the pyrotechnic detachment of the fairing of the European heavy-lift launch vehicle Ariane 5 (and 6 in the future). The HSS3+ test is not an ad-hoc academic problem created to investigate performances of the chosen method. Conversely, it is a real industrial test where many small details complicate the problem. For this reason, it is the perfect testing ground.

The scrupulous examination of the frequency response of the real data in Chapter 5 highlighted two vibration cluster at 4 and 12 kHz. The peaks around 12 kHz are due to the explosion waves that propagate throughout the rocket. These shock waves are mainly inplane and are heavily damped. Conversely, the vibration cluster around 4 kHz represents the structure modal behavior. It is mainly out-of-plane and persists in time. For these reasons, the proposed input function that characterizes the explosion loads presents a peak around 12 kHz and is similar to a gain at 4 kHz. The study of the real data provided also a qualitative value for the damping coefficient using a mixed time-and frequency-domain approach. The VTCR study was compared with real data in Section 5.4. Even if differences exist, the VTCR detected the most important frequency peaks in the interesting frequency range confirming once more its effectiveness in the mid-frequency range.

Sections 4.3.2 and 4.3.4 and Chapter 5 in particular on the HSS3+ test highlighted that the VTCR simulations are highly sensitive to the approximations introduced to address conic structures and shells with variable thickness. The average values of thicknesses and radii introduced into the VTCR model should be carefully chosen to fit the real data. The present thesis used integral means to evaluate these values. Accuracy could be increased using different averaging algorithms. As a possible future VTCR development, we suggest the study of an optimal (in some sense) algorithm to determine these approximation values. Another important aspect that should be addressed in a future VTCR development are the model uncertainties since they are relevant in the mid-frequency range. From the point of view of performances, we suggest the development of:

• a more robust computation algorithm to further mitigate the ill-conditioning problems of the method,

• a more efficient computation algorithm to study large frequency bands.

Further theoretical and experimental studies of the explosion dynamics are recommended to better understand the vibration cluster at 12 kHz. As suggested before, it is determined by the local physical processes that occur during explosions. These phenomena concern and are not limited to:

• large displacements and deformations, Simulation of impact response using a frequency approach • thermal stresses,

• elastoviscoplastic materials.

The aim of the present thesis was the extension of the VTCR to the complex industrial problems. Even if many improvements were introduced, many more are still required to fully simulate and understand the rapid-decaying processes that occur during a shock. In particular, the non-linearities that can be easily described in the time domain pose major difficulties for our frequency domain approach. The challenge is to develop an adhoc algorithm that iterates between time and frequency domain and presents acceptable computational costs. In this regard, we suggest the coupling of the VTCR with a ROM method.

  response function divided in low-mid-and highfrequency zones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 The Ariane 5 composed of the Criogenic Main Stage (CMS), the upper stage called the Storable Propellant Stage (SPS), the Solid Rocket Boosters (SRB), the payload (composed of two satellites), and the fairing. . . . 5 3 Fairing detachment scheme. . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1 Shape function types in 2D: rays (Dirac delta series), amplitude sectors, and Fourier series. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.1 Generic frame structure described in Section 2.1. A Koiter's shell Ω e is connected to some shell elements along Γ e . On the structure are applied loads (p e , m e ) and displacement contraints (u e , w e, ne ) along boundaries.In the same way, corners C e are subject to out-of-plane loads (q C e ) and displacements (w C e ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Boundary, corner, coupling and surface constraints of the generic frame structure of Figure 2.1 subdivided and hightlighted. . . . . . . . . . . . . 31 2.3 n Γe + 1 sub-domains (Ω e included) sharing the same boundary Γe. . . . . 35 2.4 n C e + 1 sub-domains (Ω e included) that share the same corner C e. . . . . 35 3.1 Qualitative behavior of propagative and evanescent waves described in Section 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 A bilinear form B and its generic block B e e . A detailed description of how the bilinear form is constructed is provided in Section 3.1. . . . . . . 50 3.3 Geometric dimensions and boundary conditions of the steel rectangular plate presented in Section 3.4.3. . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Last step of the VTCR convergence process. The VTCR solution converges at 70 rays since err V TCR ≈ 0.01. More details are reported in Section 3.4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5 VTCR and FEM displacements of the example exposed in Section 3.4.3. . 58 3.6 VTCR displacement magnitudes in case of different domain sub-divisions. The example is described in Section 3.4.3. . . . . . . . . . . . . . . . . . 59 3.7 Geometric dimensions and boundary conditions of the curved shallow shell described in Section 3.4.4. . . . . . . . . . . . . . . . . . . . . . . 60 Simulation of impact response using a frequency approach 3.8 VTCR convergence process. It is a multi-step approach described in Section 3.4.2. The specific problem treated is in Section 3.4.4. . . . . . . . . 3.9 Comparison of |w| and |N| of the vibrational problem explained in Section 3.4.4. It is based on absolute values to highlight the differences. . . . 3.10 Boundary conditions and geometric dimensions of the complex frame structure described in Section 3.4.5. . . . . . . . . . . . . . . . . . . . . 3.11 VTCR convergence process. It is a multi-step approach where ray number types are increased one-by-one up to convergence. It is described in Section 3.4.2. The specific problem treated is in Section 3.4.4. . . . . . . 3.12 Displacement magnitude comparison of the vibrational problem of Section 3.4.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Performance study of the numerical example described in Section 3.4.4. . 3.14 Effect of symmetric distribution policy during discretization refinement example explained in Section 3.6. . . . . . . . . . . . . . . . . . . . . . 3.15 Quasi-symmetric ray distribution for the first eight rays as described in Section 3.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Convergence and conditioning number versus DoFs of a typical VTCR test. The figure is further discussed in Section 3.7. . . . . . . . . . . . . . 3.17 Amplitude portraits and displacement magnitude plots of the same problem with different DoF numbers. Convergence criteria are met and both final solutions are indistinguishable. Yet, the amplitude potrait is instable as DoF number changes. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Amplitude portraits of the nornal-stress-carrying propagative in-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs. . . 3.19 Amplitude portraits of the shear stress-carrying propagative in-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs. . . 3.20 Boundary conditions and sub-domain division of the rectangular plate presented in Section 3.4.3. . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Amplitude portraits of the four sub-domains of the out-of-plane waves computed with pinv and lsqr at 100 and 200 DoFs. . . . . . . . . . . . . 3.22 Amplitude portraits of the out-of-plane waves computed with gmres and backslash at 100 and 200 DoFs. . . . . . . . . . . . . . . . . . . . . . . 3.23 Amplitude portraits of out-of-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs. . . . . . . . . . . . . . . . . . . 3.24 Different sub-division methods of the numerical example described in Section 3.4.4 and resumed in Section 3.7.1 to compare p-and ph-refinements. 80 3.25 ph-refiment with the sub-division method depicted in Figure 3.25. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. With the increasing of sub-domains, convergence is reached with less rays per sub-domain but more total rays. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Sub-domain projected areas tilting effect. It is described in detail in Section 3.7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation of impact response using a frequency approach List of Figures vii 3.27 ph-refiment with the sub-division method depicted in Figure 3.27. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. Even if the VTCR solution deteriorates as the number of sub-divisions increases, each p-refinement converges to a different solution. It is illustrated in Figure 3.28. Such behavior is explained in Section 3.7.1. . . . . . . . . . . . . . . . . . . . . . . . . . 3.28 Final VTCR displacement solutions generated using the sub-division illustrated in Figure 3.27. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. Even if the VTCR solution deteriorates as the sub-divisions increase, each p-refinement converges (to a different solution). Such behavior is explained in Section 3.7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Section of a generic cylinder sector and its underling plane. Angles are trigonometrically related since θ = 2ψ. . . . . . . . . . . . . . . . . . . . 3.30 Geometry and boundary conditions of the example illustrated in Section 3.7.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.31 VTCR and FEM solutions of the problem described in Section 3.7.2. . . .4.1 Comparison of |w| and |N| of the vibrational problem explained in Section 3.4.4 and reproposed in Section 4.1.1. It is based on absolute values to highlight the differences. . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Domain sub-division of the problem exposed in Section 4.1.1. . . . . . . 4.3 Convergence analysis of the out-of-plane rays of the numerical test in Section 4.1.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Comparison between VTCR with 1 sub-domain and 4 sub-domains of the problem exposed in Section 4.1.1. . . . . . . . . . . . . . . . . . . . . . 4.5 Boundary and geometric dimensions of the frame structure described in Section 4.1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 three-step convergence process of VTCR with in-plane inertia of the numerical example described in Section 3.1.1 at 2000 Hz. This process is explained in general in Section 3.4.2. . . . . . . . . . . . . . . . . . . . . 4.7 Triple comparison of displacement magnitudes of VTCR with and without in-plane inertia and a FEM reference of the numerical example described in Section 4.1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 5-sub-domain division of the supported half cylinder described in Section 4.1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Comparison between VTCR solution with different domain sub-divisions.

  4.2. . . . . . . . . . . . . . . . . . . . . 4.13 Comparison of displacement magnitude portraits. . . . . . . . . . . . . . 4.14 Generic structure of a symmetric skin-core-skin orthotropic sandwich shell with fundamental material properties. Section 4.3 describes the homogenization process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 A generic conic structure and a conic sector. The apex angle ξ, the curvilinear coordinates α and β, the Lamé parameter L α , and the curvature radius R α are highlighted. L β and R β are omitted because trivial. The conic sector illustrates the bottom and top values of L α and R α . Section 4.3 further discusses conic structures. . . . . . . . . . . . . . . . . . . . . . 4.16 Geometry and boundary conditions of the numerical test presented in Section 4.3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Three-step VTCR convergence process at [1000,2000,4000] Hz. It is further described in general in Section 3.4.2. The error threshold is err V TCR ≤ 0.01. Convergence results and the initial guess VTCR are reported in Table 4.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 1000 Hz with various θ a . . . . . . . . . . . . . . . . . . . . 4.19 Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 2000 Hz with various θ a . . . . . . . . . . . . . . . . . . . . 4.20 Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 4000 Hz with various θ a . . . . . . . . . . . . . . . . . . . . 4.21 err FEM with frequency and θ a of the problem in Section 4.3.2. err FEM is presented in Section 3.4.2. . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Performance study of the numerical example described in Section 4.3.2. . 4.23 Geometric dimensions and boundary conditions of the shell sector described in Section 4.3.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 FEM solutions and the VTCR approximation of the vibrational problems described in Section 4.3.4 for c h = [0, 1 4 , 1 2 , 3 4 , 1] when the thickness is described by Equation (4.26). . . . . . . . . . . . . . . . . . . . . . . . . 4.25 FEM solutions and the VTCR approximation of the vibrational problems described in Section 4.3.4 for c h = [0, 1 4 , 1 2 , 3 4 , 1] when the thickness is described by Equation (4.27). . . . . . . . . . . . . . . . . . . . . . . . . 4.26 err FEM between FEM and VTCR solutions for c h = [0, 1 4 , 1 2 , 3 4 , 1] when thickness is described by either Equation (4.26) or Equation (4.27). . . . . Simulation of impact response using a frequency approach 4.27 A Generic system in control theory. The box with the state-variables represents the system. Inputs modify state-variables that control the output. Input, output, and state-variables are called signals. A more detailed description is provided in Section 4.4. . . . . . . . . . . . . . . . . . . . . 4.28 Boundary conditions of the numerical test presented in Section 4.4.3. Conditions on u and w along the left edge are alike for the sake of simplicity. The imaginary part leads the delay. . . . . . . . . . . . . . . . . . 4.29 VTCR u and w displacements at different frequencies of the example illustrated in Section 4.4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.30 VTCR u and w displacements at different frequencies. . . . . . . . . . . . 4.31 Geometry and boundary conditions of the numerical test presented in Section 4.5.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.32 VTCR convergence process. It is a multi-step approach described in Section 3.4.2. Since the problem is completely out-of-plane, in-plane propagative rays are unnecessary. The specific problem treated is in Section 4.5.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.33 VTCR and FEM displacement magnitudes of the example exposed in Section 4.5.1. The difference between the two solutions is err FEM ≈ 6%. . . 4.34 C ∂ i b Ω e matrix. Its definition is provided in Section 4.6. . . . . . . . . . . .

  5.5 A generic frame structure and the generic case of a boundary shared among multiple sub-domains. . . . . . . . . . . . . . . . . . . . . . . . . 5.6 General triangle function ∧ and its parameters {a ∧ , b ∧ , h ∧ }. . . . . . . . . 5.7 The pyrotechnic explosion as percieved by a sensor in [a s , b s ] . . . . . . . 5.8 Sensor positions of the HSS3+ test. . . . . . . . . . . . . . . . . . . . . . 5.9 Time and frequency domain responses of an accelerometer near the explosion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Time and frequency domain responses of an accelerometer near the satellite.146 5.

  . . . . . . . . . . . . . . . . . . . . . . . 5.25 SRS comparisons between VTCR simulations where the variable-thickness shells are approximated with Equation (4.25) (reference model) or Equation (5.42) (modified model) for a sensor near the explosion, the satellite, and in the middle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.26 SRS comparisons between VTCR simulations where the variable-thickness shells are approximated with Equation (4.25) (reference model) or Equation (5.42) (modified model) for a sensor near the explosion, the satellite, and in the middle. A semilog scale is used to highlight magnitude differences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation of impact response using a frequency approach List of Tables 3.1 Wave identification using the characteristics of k i , ki , and âi = [a xy , a z ] . . 3.2 Characteristics of the workstation used for the academic tests. . . . . . . 3.3 Quantities of interest of the example in Section 3.4.3. . . . . . . . . . . .

  3.2 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.9 Quantities of interest of the example in Section 4.3.4. . . . . . . . . . . . 117 4.10 Quantities of interest of the example in Section 4.4.3. . . . . . . . . . . . 123 4.11 Quantities of interest of the example in Section 4.5.1. . . . . . . . . . . . 128 4.12 Ray number and type for initial guess VTCR and converged VTCR. converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 4.5.1 where the error threshold is err V TCR ≤ 0.01. . . . . . . . . . . 129 4.13 DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 4.5.1. . . . . . . . . 129 5.1 The highest speeds of sound for in-plane and out-of-plane waves along α and y at 15000 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.2 Matched explosion input and damping coefficients for the HSS3+ test. The matching processes are illustrated in Sections 5.3.4 and 5.3.5. . . . . 153 5.3 Characteristics of the cluster used to compute the VTCR solution in Section 5.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.4 Ray number and type for initial guess VTCR and converged VTCR at the highest interesting frequency (5000 Hz) for the studied vibrational problem is illustrated in Section 5.4. The error threshold for the convergence process is err V TCR ≤ 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . 153

Figure 2 :Figure 3 :

 23 Figure 2: The Ariane 5 composed of the Criogenic Main Stage (CMS), the upper stage called the Storable Propellant Stage (SPS), the Solid Rocket Boosters (SRB), the payload (composed of two satellites), and the fairing.

Figure 1 . 1 :

 11 Figure 1.1: Shape function types in 2D: rays (Dirac delta series), amplitude sectors, and Fourier series.

Figure 2 . 1 :

 21 Figure2.1: Generic frame structure described in Section 2.1. A Koiter's shell Ω e is connected to some shell elements along Γ e . On the structure are applied loads (p e , m e ) and displacement contraints (u e , w e, ne ) along boundaries. In the same way, corners C e are subject to out-of-plane loads (q C e ) and displacements (w C e )

Figure 2 . 2 :

 22 Figure 2.2: Boundary, corner, coupling and surface constraints of the generic frame structure of Figure 2.1 subdivided and hightlighted.

Figure 2

 2 Figure 2.3: n Γe + 1 sub-domains (Ω e included) sharing the same boundary Γe.

Figure 2

 2 Figure 2.4: n C e + 1 sub-domains (Ω e included) that share the same corner C e.

  (a) Propagative wave. (b) Evanescent wave.

Figure 3

 3 Figure 3.1: Qualitative behavior of propagative and evanescent waves described in Section 3.1.

  A generic block B e e of the bilinear form B.

Figure 3

 3 Figure 3.2: A bilinear form B and its generic block B e e . A detailed description of how the bilinear form is constructed is provided in Section 3.1.

  Figure 3.3: Geometric dimensions and boundary conditions of the steel rectangular plate presented in Section 3.4.3.

Figure 3

 3 Figure 3.4: Last step of the VTCR convergence process. The VTCR solution converges at 70 rays since err V TCR ≈ 0.01. More details are reported in Section 3.4.3.

  (a) VTCR u displacement. (b) FEM u displacement. (c) VTCR v displacement. (d) FEM v displacement. (e) VTCR displacement magnitude. (f) FEM displacement magnitude.

Figure 3

 3 Figure 3.5: VTCR and FEM displacements of the example exposed in Section 3.4.3.

Figure 3

 3 Figure 3.6: VTCR displacement magnitudes in case of different domain sub-divisions. The example is described in Section 3.4.3.

Figure 3

 3 Figure 3.7: Geometric dimensions and boundary conditions of the curved shallow shell described in Section 3.4.4.

  (a) Particular solution rays. (b) Out-of-plane propagative rays.

Figure 3

 3 Figure 3.8: VTCR convergence process. It is a multi-step approach described in Section 3.4.2. The specific problem treated is in Section 3.4.4.

  FEM |w| displacement. (d) VTCR |N x | stress resultant. (e) reference |N x | stress resultant. (f) FEM |N x | stress resultant. (g) VTCR |N y | stress resultant. (h) reference |N y | stress resultant (i) FEM |N y | stress resultant. (j) VTCR |N xy | stress resultant (k) reference |N xy | stress resultant. (l) FEM |N xy | stress resultant.

Figure 3

 3 Figure 3.9: Comparison of |w| and |N| of the vibrational problem explained in Section 3.4.4. It is based on absolute values to highlight the differences.

  Figure 3.10: Boundary conditions and geometric dimensions of the complex frame structure described in Section 3.4.5.

Figure 3 .

 3 Figure 3.11: VTCR convergence process. It is a multi-step approach where ray number types are increased one-by-one up to convergence. It is described in Section 3.4.2. The specific problem treated is in Section 3.4.4.

Figure 3 .

 3 Figure 3.12: Displacement magnitude comparison of the vibrational problem of Section 3.4.5.

Figure 3 .

 3 Figure 3.13: Performance study of the numerical example described in Section 3.4.4.

Figure 3

 3 Figure 3.14: Effect of symmetric distribution policy during discretization refinement example explained in Section 3.6.

Figure 3

 3 Figure3.15: Quasi-symmetric ray distribution for the first eight rays as described in Section 3.6.

Figure 3 .

 3 Figure 3.16: Convergence and conditioning number versus DoFs of a typical VTCR test. The figure is further discussed in Section 3.7.

  (a) Amplitude portrait with 100 DoFs. (b) Amplitude portrait with 200 DoFs. (c) Displacement magnitude with 100 DoFs (|w 100 |).

Figure 3 . 17 :

 317 Figure 3.17: Amplitude portraits and displacement magnitude plots of the same problem with different DoF numbers. Convergence criteria are met and both final solutions are indistinguishable. Yet, the amplitude potrait is instable as DoF number changes.

  4.3 presents multiple in-plane ray types and domain sub-divisions. For the sake of clarity and simplicity the study is restricted to the one-domain version. The considered discretizations are at 100 and 200 DoFs. They are both well beyond convergence (70 DoFs). Figures 3.18 and 3.19 present the results for the normal and shear stress-carrying propagative waves.

  gmres, 100 DoFs. (d) backslash, 100 DoFs. (e) pinv, 200 DoFs. (f) lsqr, 200 DoFs. (g) gmres, 200 DoFs. (h) backslash, 200 DoFs.

Figure 3 .

 3 Figure 3.18: Amplitude portraits of the nornal-stress-carrying propagative in-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs.

  gmres, 200 DoFs. (h) backslash, 200 DoFs.

Figure 3 .

 3 Figure 3.19: Amplitude portraits of the shear stress-carrying propagative in-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs.

Figure 3 .

 3 Figure 3.20: Boundary conditions and sub-domain division of the rectangular plate presented in Section 3.4.3.

Figure 3 .

 3 Figure 3.21: Amplitude portraits of the four sub-domains of the out-of-plane waves computed with pinv and lsqr at 100 and 200 DoFs.

Figure 3 .

 3 Figure 3.22: Amplitude portraits of the out-of-plane waves computed with gmres and backslash at 100 and 200 DoFs.

  (a) pinv, 100 DoFs. (b) lsqr, 100 DoFs. (c) gmres, 100 DoFs. (d) backslash, 100 DoFs. (e) pinv, 200 DoFs. (f) lsqr, 200 DoFs. (g) gmres, 200 DoFs. (h) backslash, 200 DoFs.

Figure 3 .

 3 Figure 3.23: Amplitude portraits of out-of-plane waves computed with pinv, lsqr, gmres, and backslash at 100 and 200 DoFs.

  p-iteration, 4 sub-domains. (d) p-iteration, 8 sub-domains.

Figure 3 .

 3 Figure 3.25: ph-refiment with the sub-division method depicted in Figure 3.25. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. With the increasing of sub-domains, convergence is reached with less rays per sub-domain but more total rays.

Figure 3

 3 Figure 3.26: Sub-domain projected areas tilting effect. It is described in detail in Section 3.7.1.

  (a) p-iteration, 1 domain. (b) p-iteration, 2 sub-domains. (c) p-iteration, 4 sub-domains. (d) p-iteration, 8 sub-domains.

Figure 3 .

 3 Figure 3.27: ph-refiment with the sub-division method depicted in Figure 3.27. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. Even if the VTCR solution deteriorates as the number of sub-divisions increases, each p-refinement converges to a different solution. It is illustrated in Figure 3.28. Such behavior is explained in Section 3.7.1.

  (a) quasi-analytic reference solution |w|. (b) VTCR solution, 2 sub-domains |w|. (c) VTCR solution, 4 sub-domains |w|. (d) VTCR solution, 8 sub-domains |w|.

Figure 3 .

 3 Figure 3.28: Final VTCR displacement solutions generated using the sub-division illustrated in Figure 3.27. The domain is divided in 1, 2, 4, 8 sub-domains and the number of rays is increased up to convergence. Even if the VTCR solution deteriorates as the subdivisions increase, each p-refinement converges (to a different solution). Such behavior is explained in Section 3.7.1.

Figure 3 .

 3 Figure 3.30: Geometry and boundary conditions of the example illustrated in Section 3.7.2.

Figure 3 .

 3 Figure 3.31: VTCR and FEM solutions of the problem described in Section 3.7.2.

  FEM |w| displacement. (d) reference |N x | stress resultant. (e) VTCR |N x | stress resultant. (f) FEM |N x | stress resultant. (g) reference |N y | stress resultant. (h) VTCR |N y | stress resultant (i) FEM |N y | stress resultant. (j) reference |N xy | stress resultant (k) VTCR |N xy | stress resultant. (l) FEM |N xy | stress resultant.

Figure 4 Figure 4

 44 Figure 4.1: Comparison of |w| and |N| of the vibrational problem explained in Section 3.4.4 and reproposed in Section 4.1.1. It is based on absolute values to highlight the differences. Simulation of impact response using a frequency approach

Figure 4 . 3 :

 43 Figure 4.3: Convergence analysis of the out-of-plane rays of the numerical test in Section 4.1.1.

  (a) VTCR 1 sub-domain |w| displacement. (b) VTCR 4 sub-domains |w| displacement.

Figure 4

 4 Figure 4.4: Comparison between VTCR with 1 sub-domain and 4 sub-domains of the problem exposed in Section 4.1.1.

Figure 4 .

 4 Figure 4.3 further validates our conclusions about p-and ph-refinements illustrated at the end of Section 3.7.1.

Figure 4

 4 Figure 4.5: Boundary and geometric dimensions of the frame structure described in Section 4.1.2.

  (a) First step: propagative out-of-plane ray convergence. (b) Second step: evanescent out-of-plane ray convergence. (c) Third step: propagative in-plane ray convergence (normal and shear stress carrying rays).

Figure 4

 4 Figure 4.6: three-step convergence process of VTCR with in-plane inertia of the numerical example described in Section 3.1.1 at 2000 Hz. This process is explained in general in Section 3.4.2.

  Figure 4.5 while Figure 4.9b is divided in 5 sub-domains as depicted in Figure 4.8. VTCR and FEM performances are tested over a frequency band f ∈ [2000; 4000] Hz with a 100 Hz step. Performance analyses over the frequency band [2000; 4000] Hz are reported in Figure 4.10b and Figure 4.10c. The number of DoF of the FEM reference is chosen so that the rule-of-thumb in Equation (3.48) suggested for mid-frequency problems is met.

  converged VTCR with in-plane inertia. (c) converged VTCR without in-plane inertia.

Figure 4 . 7 :

 47 Figure 4.7: Triple comparison of displacement magnitudes of VTCR with and without inplane inertia and a FEM reference of the numerical example described in Section 4.1.2.

Figure 4

 4 Figure 4.8: 5-sub-domain division of the supported half cylinder described in Section 4.1.2.

Figure 4

 4 Figure 4.9: Comparison between VTCR solution with different domain sub-divisions. This demonstrate that the degradation effect highlighted in Sections 3.7.1 and 3.7.2 is solved. The problem is described in Section 4.1.2.

Figure 4 .

 4 Figure 4.10: Performance analyzes over a frequency band[2000; 4000] Hz of the numerical example described in Section 4.1.2. VTCR and FEM are compared considering time and memory consumptions. DoFs required are reported for the sake of completeness, a direct comparison is unmeaningful due to theory differences.

  ki (b) Wave distribution in an orthotropic material.

Figure 4 .

 4 Figure 4.11: Changement of ki to address different wave speeds in orthotropic materials.

  (a) First step: propagative out-of-plane ray convergence. (b) Second step: evanescent out-of-plane ray convergence. (c) Third step: propagative in-plane ray convergence (normal and shear stress carrying rays).

Figure 4 .

 4 Figure 4.12: three-step convergence process of VTCR with in-plane inertia of the numerical example described in Section 4.2 at 3700 Hz. This process is explained in general in Section 3.4.2.

Figure 4 .

 4 Figure 4.13: Comparison of displacement magnitude portraits.

  Figure 4.14: Generic structure of a symmetric skin-core-skin orthotropic sandwich shell with fundamental material properties. Section 4.3 describes the homogenization process.

  to 4.20. Figure4.21 summarizes the difference between the reference and the VTCR solutions for various θ a and frequencies. The used indicator is err FEM as discussed in Section 3.4.2. Performances are exposed in Figure4.22.

Figure 4 .

 4 Figure 4.16: Geometry and boundary conditions of the numerical test presented in Section 4.3.2.

  (a) Out-of-plane propagative rays -1000 Hz. (b) Out-of-plane evanescent rays -1000 Hz. (c) In-plane propagative rays -

Figure 4 .

 4 Figure 4.17: Three-step VTCR convergence process at [1000,2000,4000] Hz. It is further described in general in Section 3.4.2. The error threshold is err V TCR ≤ 0.01. Convergence results and the initial guess VTCR are reported inTable 4.8.

  (a) FEM referenceθ a = 0. (b) VTCR solutionθ a = 0. (c) FEM solution referenceθ a = π/16. (d) VTCR solutionθ a = π/16. (e) FEM solution referenceθ a = π/8. (f) VTCR solutionθ a = π/8.

Figure 4 .

 4 Figure 4.18: Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 1000 Hz with various θ a .

  (a) FEM solution referenceθ a = 0. (b) VTCR solutionθ a = 0. (c) FEM solution referenceθ a = π/16. (d) VTCR solutionθ a = π/16. (e) FEM solution referenceθ a = π/8. (f) VTCR solutionθ a = π/8.

Figure 4 .

 4 Figure 4.19: Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 2000 Hz with various θ a .

  (a) FEM solution referenceθ a = 0. (b) VTCR solutionθ a = 0. (c) FEM solution referenceθ a = π/16. (d) VTCR solutionθ a = π/16. (e) FEM solution referenceθ a = π/8. (f) VTCR solutionθ a = π/8.

Figure 4 .

 4 Figure 4.20: Displacement magnitude comparison of the vibrational problem of Section 4.3.2 at 4000 Hz with various θ a .

Figure 4 .

 4 Figure 4.21: err FEM with frequency and θ a of the problem in Section 4.3.2. err FEM is presented in Section 3.4.2.

Figure 4 .

 4 Figure 4.22: Performance study of the numerical example described in Section 4.3.2.

Figure 4 .

 4 Figure 4.23: Geometric dimensions and boundary conditions of the shell sector described in Section 4.3.4.

Figure 4 .

 4 Figure 4.25: FEM solutions and the VTCR approximation of the vibrational problems described in Section 4.3.4 for c h = [0, 1 4 , 1 2 , 3 4 , 1] when the thickness is described by Equation (4.27).

Figure 4 .

 4 Figure 4.26: err FEM between FEM and VTCR solutions for c h = [0, 1 4 , 1 2 , 3 4 , 1] when thickness is described by either Equation (4.26) or Equation (4.27).

Figure 4

 4 Figure 4.27: A Generic system in control theory. The box with the state-variables represents the system. Inputs modify state-variables that control the output. Input, output, and state-variables are called signals. A more detailed description is provided in Section 4.4.

Frequency-domain approaches

  investigate the problem in the frequency-domain. Each time-function is transformed by mean of the Laplace transform in its corresponding Laplace function s indicates the Laplace transform of the function and s = σ c + ıω is a complex frequency number. The Fourier transform Equation (2.1) is a particular case of the Laplace transform with σ c = 0. In this way, the state-space system in Equation (4.31) becomes sx cs (s) = A c x cs (s) + B c u cs (s) y cs (s) = C c x cs (s) + D c u cs (s) (4.47)

  Figure 4.28. This time the right boundary is clamped and the left edge is subject to in-plane and out-of-plane sinusoidal displacements f s which are alike for the sake of simplicity. They are illustrated in Figure 4.28b. In the time domain the boundary sinusoidally vibrates at the prescribed frequency. The real part represents the maximum displacement while the imaginary part is related to the delay.

  u = 0 w, n = 0 u = [f s , 0, f s ] (a) Boundary conditions.(b) Boundary condition f s along the left edge.

Figure 4 .

 4 Figure 4.28: Boundary conditions of the numerical test presented in Section 4.4.3. Conditions on u and w along the left edge are alike for the sake of simplicity. The imaginary part leads the delay.

  Figure 4.29 illustrates the VTCR displacements at various frequencies. Since the displacement magnitude is constant along the left boundary, Figure 4.29 reports u and w real displacements. Finally, Figure 4.30 compares imaginary and real parts of the VTCR boundary conditions with the imposed values. Since VTCR is a Trefftz method, boundary conditions are approximated.

  .33) to (2.35), this condition is e = econcentrated masses ρ b and springs κ b related to and distributed along the boundary, two additional loads arise so that Equation (4.58) becomes e is the dual displacement quantity of e , -ρ b ω 2 e is the force of inertia of the concentrated mass, and -κ b e is the concentrated spring load. Equation (4.59) proves that concentrated masses and springs behave alike. In fact, inertia and spring loads can be combined in a single term µ

  Figure 4.31: Geometry and boundary conditions of the numerical test presented in Section 4.5.1.

  (a) Out-of-plane propagative rays. (b) Out-of-plane evanescent rays.

Figure 4 .

 4 Figure 4.32: VTCR convergence process. It is a multi-step approach described in Section 3.4.2. Since the problem is completely out-of-plane, in-plane propagative rays are unnecessary. The specific problem treated is in Section 4.5.1.

Figure 4 .

 4 Figure 4.33: VTCR and FEM displacement magnitudes of the example exposed in Section 4.5.1. The difference between the two solutions is err FEM ≈ 6%.

  p 21 = p 2p 1 It is also p 21 = p 2relep 1rele = p 2relep 1rele since the boundary is shared between Ω e and Ω e . Changing the integration variable r reli = p 21 ς + p 1rele (4.68) r rel j = p 21 ς + p 1rele (4.69) ds = L p 21 dς (4.70)

Figure 4 .

 4 Figure 4.34: C ∂ i b Ω e matrix. Its definition is provided in Section 4.6.

  proper dispersion equation in Equation (5.1) provides the speed of sound. Shells present different speeds of sound for in-plane c gip and out-of-plane c goop vibrations. Hereafter subscripts ip and oop refer to in-plane and out-of-plane behaviors. Broadly, for homogeneous isotropic plates

Figure 5 . 2 :

 52 Figure 5.2: HSS3+ structure before and after the pyrotechnic explosions.

  explosion at a generic instant t.

Figure 5 . 3 :

 53 Figure 5.3: Figure 5.3a shows the rocket section during the pyrotechnic detachment. Detonations ignite in A and B and end in C and D following the four explosion paths indicated by the red arrows. The explosions are supersonic. Figure 5.3b highlights their Mach cones. The PFL cylinder is detached from the rocket More details are provided in Section 5.2.

Figure 5 . 4 :

 54 Figure 5.4: Equivalent HSS3+ structure before and after the pyrotechnic explosions.

Figure 5

 5 Figure 5.5: A generic frame structure and the generic case of a boundary shared among multiple sub-domains.

Figure 5 . 6 :

 56 Figure 5.6: General triangle function ∧ and its parameters {a ∧ , b ∧ , h ∧ }.

Figure 5 .

 5 Figure 5.7 depicts the explosion and an accelerometer in a certain, distorted 2D metric. The focus is on the in-plane waves since they are the fastest information-carrying agents. The explosion starts at t = 0 in [0, 0] . It emits vibrational waves during its supersonic propagation along the explosion path. A sensor is located in [a s , b s ] . It starts to perceive the explosion when the Mach cone reaches its position.

  Figure 5.8: Sensor positions of the HSS3+ test.

Figures 5 .

 5 Figures 5.9 and 5.10 illustrate in-plane and out-of-plane sensor data in time and frequency domains for a pair of accelerometers near the explosion path and near the satellite respectively. Since real data are confidential, ticks in the acceleration axes are substituted with letters and the time scales are obscured. Sensor responses at the same altitude are similar. Accelerometer signals at intermediate altitudes present a mixed behavior between Figure 5.9 and Figure 5.10.

  (a) In-plane accelerations -time domain. (b) In-plane accelerations -frequency domain. (c) Out-of-plane accelerations -time domain. (d) Out-of-plane accelerations -frequency domain.

Figure 5 . 9 :

 59 Figure 5.9: Time and frequency domain responses of an accelerometer near the explosion.

  (a) In-plane accelerations -time domain. (b) In-plane accelerations -frequency domain. (c) Out-of-plane accelerations -time domain. (d) Out-of-plane accelerations -frequency domain.

Figure 5 .

 5 Figure 5.10: Time and frequency domain responses of an accelerometer near the satellite.

  Figures 5.12 and 5.13 consider the first milliseconds of the two previous signals where the explosion input leads the vibrational response. Conversely, Figures 5.14 and 5.15 focus on a delayed time window where the modal behavior is already settled.Remark. The vibration cluster at 12 kHz that dominates the frequency response in Figure5.12 is heavily attenuated in space (Figure5.13) and time(Figure 5.14). The com-Simulation of impact response using a frequency approach

  (a) Sensor near the explosion -in-plane. (b) Sensor near the explosion -out-of-plane. (c) Sensor near the satellite -in-plane. (d) Sensor near the satellite -out-of-plane.

Figure 5 .

 5 Figure 5.11: Time-frequency plots of a sensor near the explosion and the satellite.

  (a) In-plane accelerations -time domain. (b) In-plane accelerations -frequency domain. (c) Out-of-plane accelerations -time domain. (d) Out-of-plane accelerations -frequency domain.

Figure 5 .

 5 Figure 5.12: Time and frequency domain responses of an accelerometer near the explosion on a small time window at the beginning of the signal.

  (a) In-plane accelerations -time domain. (b) In-plane accelerations -frequency domain. (c) Out-of-plane accelerations -time domain. (d) Out-of-plane accelerations -frequency domain.

Figure 5 .

 5 Figure 5.13: Time and frequency domain responses of an accelerometer near the satellite on a small time window at the beginning of the signal.

  (a) In-plane accelerations -time domain. (b) In-plane accelerations -frequency domain. (c) Out-of-plane accelerations -time domain. (d) Out-of-plane accelerations -frequency domain.

Figure 5 .

 5 Figure 5.15: Time and frequency domain responses of an accelerometer near the satellite on a small time window at the end of the signal.

  Input function -frequency domain.

Figure 5 .

 5 Figure 5.16: Qualitative behavior of the input function in time and frequency domain determined studying the HSS3+ real data.

  Figure 5.11), Klein assumes a Rayleigh model for the damping matrix C RADIOSS C RADIOSS = α RADIOSS K RADIOSS + β RADIOSS M RADIOSS (5.27)

η

  = a η ω -b η (5.28)where a η and b η are two empirical coefficients that depend on the used materials. Usually, b η ∈ [0; 0.7). Broadly, the in-plane k ip and out-of-plane k oop wavenumbers are

  t 1 and t 2 are the two time windows. t 1 = t 2 for the sake of simplicity. After the transient, the frequency responses of two successive time windows are illustrated in Figures 5.17b and 5.17c. Figures 5.17b and 5.17c are similar. Many peaks in Figure 5.17b persist in Figure 5.17c. Let us consider the high a pi of one of these peaks in Figure 5.17b. After a t1 t2 (a) Time domain. (b) Frequency response of the first time window. (c) Frequency response of the second time window.

Figure 5 .

 5 Figure 5.17: Comparison of the frequency responses on different time windows of a sensor near the satellite to determine the damping coefficient η.

  Figure 5.18. The first instants are neglected to study the modal behavior in the interesting frequency range. Figures 5.18 and 5.19 compare real data and VTCR results. VTCR frequency responses are multiplied by exp(-ηωt) to take into account the time dependency. Color-bars are confidential. However, real data plots and their corresponding VTCR results share the same color-bar.Simulation of impact response using a frequency approach

Figure 5 .

 5 Figure 5.19: Out-of-plane time-frequency plot comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle.

  err f = |a V TCR | -|a re f | k a re f k (5.33)where a V TCR and a re f are the accelerations of the VTCR and the real data results respectively and | | is the absolute value. In order to study the influence of the frequency on the error, four frequency bands are investigated:

  and 5.21. The error indicator err f is based on the kinetic energy and computed on four different frequency bands to determine the influence of the frequency on the error.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

Figure 5 .

 5 Figure 5.20: Frequency response comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

Figure 5 .

 5 Figure 5.21: Frequency response comparisons in a semilog scale between VTCR and real data for a sensor near the explosion, the satellite, and in the middle.

Figure 5

 5 Figure 5.22: Generic SDOF of the SRS illustrated in Section 5.4.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

Figure 5 .

 5 Figure 5.23: SRS comparisons between VTCR and real data for a sensor near the explosion, the satellite, and in the middle.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

Figure 5 .

 5 Figure 5.25: SRS comparisons between VTCR simulations where the variable-thickness shells are approximated with Equation (4.25) (reference model) or Equation (5.42) (modified model) for a sensor near the explosion, the satellite, and in the middle.

  (a) Sensor near the explosion -in-plane response. (b) Sensor near the explosion -out-of-plane response. (c) Sensor near the satellite -in-plane response. (d) Sensor near the satellite -out-of-plane response. (e) Sensor in the middle -in-plane response. (f) Sensor in the middle -out-of-plane response.

Figure 5 .

 5 Figure 5.26: SRS comparisons between VTCR simulations where the variable-thickness shells are approximated with Equation (4.25) (reference model) or Equation (5.42) (modified model) for a sensor near the explosion, the satellite, and in the middle. A semilog scale is used to highlight magnitude differences.

  of unity subordinate to the cover {Ω i }. The partition of unity {c i } is said to be of degree m ∈ N 0 if {c i } ⊂ C m (R n e ). The covering sets are called patches. Let {Ω i } be an open cover of Ω ⊂ R n e and let {c i } be a (M,C ∞ ,C G ) partition of unity subordinate to {Ω i }. Let V i ⊂ H 1 (Ω i ⊂ Ω) be given. Than the space

	Definition 2.

Table 3 .

 3 1: Wave identification using the characteristics of k i , ki , and âi = [a xy , a z ] .

  . It addresses boundary and corner conditions in a weak sense. Conversely to its definition in other works, it explicitly addresses the general case of boundaries and corners shared among multiple sub-domains (more than two) with applied loads. Moreover, this formulation is easily implementable in a computer program. First, the weak variational problem is defined. After, its key properties are enunciated and demonstrated. Definition 8 (The weak variational problem). Find the solution sets {u es , N es , M es } ∈ D e Existence and uniqueness properties Theorem 3. The weak variational problem in Equation (3.21) is equivalent to the general problem described in Chapter 2 in weak sense and the solution of Equation (3.21) exists and is unique. Proof. The solution of the vibrational problem illustrated in Chapter 2 satisfies Equation (3.21). Hence, if the solution of Equation (3.21) exists and is unique, it is the solution of the vibrational problem of Chapter 2. The uniqueness property is proved by contradiction. Let us assume that two different solutions D s1 = {u s1 , N s1 , M s1 } and D s2 = {u s2 , N s2 , M s2 } of the weak form in Equation (3.21) exist over the problem domain Ω. The difference

	such that
			n e e=1 ∑		∂ v e Ω e	δb H
							∂∂ w C e Ω e
	+	n Γe n Γe + 1 Γe	δb H e v es + δq H e w es -δm H e w es , ne ds
	-	1 n Γe + 1 Γe e b -δv H n C e n C e + 1 δq H C e w C es	C e	+	1 n C e + 1	δw H C C e
	-	n ∑ e =1, e =e	1 n Γe + 1 Γe	δb H e v e s + δq H e w e s + δm H e w e s , ne ds
			+	1 n Γe + 1 Γe	δv H e b e s + δw H e q e s + δw e , H ne m e s ds
			-	1 n C e + 1	δq H C e w C e s + δw H C e q C e s C e	= 0
					∀ {δu

e (v esv e ) ds -∂ b e Ω e δv H e b esb e ds + ∂ w e Ω e δq H e (w esw e ) ds -∂ q e Ω e δw H e (q esq e ) ds -∂ w e , ne Ω e δm H e (w es , ne -w e , ne ) ds + ∂ m e Ω e δw e , nH e (m esm e ) ds

+ δw H C e (q C esq C e ) ∂∂ q C e Ω eδq H C e (w C esw C e ) esb e + δw H e (q esq e )δw e , H ne (m esm e ) ds e (q C esq C e ) e , δN e , δM e } ∈ δD e (3.21)

where n is the number of sub-domains and δD e = D e0 is the test function space being the VTCR a Galerkin method.

Simulation of impact response using a frequency approach 3.2.1

  solution of the difference problem. Its equilibrium equations are homogeneous and all boundary, corner, and coupling conditions are zero since both solutions solve the very same vibrational problem. For this reason, δD d = D d0 = D d and the weak form becomes

	n e ∑ e=1	∂ v ed Ω e	b H ed v ed ds -	∂ b ed	Ω e	v H ed b ed ds
	+	∂ w ed Ω e	q H ed w ed ds -	∂ q ed Ω e	w H ed q ed ds
	-	∂ w ed , ne	Ω e	m H ed w ed , ne ds +

∂ m ed Ω e w ed , nH e m ed ds

  ∪ ∂b ed Ω e = ∂w ed Ω e ∪ ∂q ed Ω e = ∂w ed , ne Ω e ∪ ∂m ed Ω e = ∂Ω e \Γ e (3.24) ∂∂q C ed Ω e ∪ ∂∂w C ed Ω e = ∂∂Ω e L exg = 0 since every external load of the difference problem is null. The dissipation energy E Dd is positive by definition. Hence, N d and M d are null. Injecting this result in the homogeneous equilibrium Equations (2.12) and (2.13) implies that u d is null everywhere. For this reason

						(3.25)
	Substituting Equation (2.58) in Equation (3.23) leads to
				E Dd = 0		(3.26)
						22)
	Without loss of generality Equation (3.22) can be restricted to focus only its imagi-
	nary part since it should be true for every test function. Considering that ℑ H	=
	H can be rearranged such that being and are two generic complex vectors or scalars, Equation (3.22) -ℑ
	ℑ	n e ∑ e=1	∂Ω e	b H ed v ed + q H ed w ed -m H ed w ed , ne ds -q H C e w C e ∂∂Ω e	= 0,	(3.23)
			∀ {u ed , N ed , M ed } ∈ D ed

where the problem is supposed well posed. In particular, boundary and corner conditions satisfy Simulation of impact response using a frequency approach ∂v ed Ω e where

  33)where L C [ ] and L C [ ] depends on a and k of the introduced function and on n1 , t1 , n2 , and t2 of the corner considered. The linear form takes into account boundary and corner conditions. It is created in a very similar way. In particular, constant boundary conditions and corner conditions are calculated as illustrated in Equation (3.33) and Equation (3.32) where the shape function part is replaced with its imposed boundary or corner value. Since variable boundary conditions are a theoretical improvement, they are explained in Section 4.4.2.

Table 3 .

 3 Table 3.2. 2: Characteristics of the workstation used for the academic tests.

	Simulation of impact response using a frequency approach

Table 3 .

 3 .6 depicts the VTCR displacement magnitude portraits in case of different domain sub-divisions. 3: Quantities of interest of the example in Section 3.4.3.

	h	thickness	3 mm
	f	frequency	5000 Hz
	E Young modulus 200 GPa
	ν Poisson's ratio	0.3
	ρ	density	7800 Kg/m 3
	η damping factor 0.01

Simulation of impact response using a frequency approach

Table 3 .

 3 4: Ray number and type for initial guess VTCR and converged VTCR. Converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.3 where the error threshold is err V TCR ≤ 0.01.

		DoFs Time [s] RAM [Mb]
	VTCR	70	0.8	2.73
	FEM	15912	1	15
	Table 3.5: DoFs, time, and memory consumption comparisons between VTCR and FEM
	of the numerical example described in Section 3.4.3.

  .9 where |w|, |N xx |, |N xy |, and |N xy | are reported. It is based on the absolute value of the quantities to highlight the differences. Since the extraction of u and v from the Airy's stress function of the reference solution is

		DoFs Time [s] RAM [Mb]
	VTCR	906	4	20
	FEM	342504	18	182
	Table 3.8: DoFs, time, and memory consumption comparisons between VTCR and FEM
	of the numerical example described in Section 3.4.4.
	Simulation of impact response using a frequency approach

Table 3

 3

		thickness	1 mm
	R radius of curvature	2.5 m
	f	frequency	500 Hz
	E	Young modulus	200 GPa
	ν	Poisson's ratio	0.33
	ρ	density	7850 Kg/m 3
	η	damping factor	0.01

.9: Quantities of interest of the example in Section 3.4.5.

Simulation of impact response using a frequency approach

Table 4 .

 4 Table4.1 reports its results. 1: Ray number and type for initial guess VTCR and converged VTCR. converged VTCR is the result of a multi-step convergence process explained in general in Section 3.4.2 and illustrated for this specific case in Section 3.4.4 where the error threshold is err V TCR ≤ 0.01.The triple comparison among FEM, VTCR, and the quasi-analytic reference is proposed in Figure4.1.Performances among old TAPYROSS, its new version with the performance improvement introduced in Section 4.6, and a FEM reference are proposed in Table4.2.

	Types	VTCR initial guess Converged VTCR
	Propagative particular solution	1681	841
	Propagative out-of-plane	200	65
	Evanescent out-of-plane	0	0
	Propagative in-plane		0	0
		DoFs	Time [s] RAM [Mb]
	new VTCR	65	0.3	9.2
	old VTCR	65	4	20
	FEM	342,504	18	182
	Table 4.2: DoFs, time, and memory consumption comparisons between new and old
	VTCR and FEM of the numerical example described in Section 3.4.4 and reanalyzed
	in Section 4.1.1.

Table 4

 4

	1]

.9: Quantities of interest of the example in Section 4.3.4.

  F is related to the Fourier series, L s is the Lamé parameter relative to s, and d s is the boundary length. The series in Equation (4.52) is truncated for computational purposes. Equation (4.52) in matrix form is

	(s) = ∑	Fi e ık Fi s		(4.52)
		i∈Z
	Fi = k Fi =	1 d s 2πi d s	d s /2 -d s /2	1 L s	e -ık Fi s ds	(4.53) (4.54)
	L s = r , s r, s			(4.55)
	where the subscript = (e ık F s ) F	(4.56)
	Simulation of impact response using a frequency approach

Table 4 .

 4 The present test validates the theory introduced in Section 4.4.2. A plate similar to Section 4.5.1 is subject to non-uniform boundary conditions at different frequencies. Since 10: Quantities of interest of the example in Section 4.4.3. the focus is on the imposed boundary conditions, the VTCR boundary values are compared with the imposed ones.Plate geometry and material properties are alike the example in Section 4.5.1 and are reported in Table4.10 for the sake of clarity. However, boundary conditions are different. They are reported in

	4.4.3 Numerical test: a plate subject to a non-uniform boundary con-
	dition
	Simulation of impact response using a frequency approach

Table 5 .

 5 Simulation of impact response using a frequency approach a bpyro [N/m] a qpyro [N/m] 2: Matched explosion input and damping coefficients for the HSS3+ test. The matching processes are illustrated in Sections 5.3.4 and 5.3.5.Computations are performed on a cluster. Its characteristics are reported in Table5.3. Table5.4 exposes results of the VTCR convergence analysis where the error threshold is err V TCR ≤ 0.01. VTCR performances and DoFs are illustrated in Table5.5.

	1230	154	a η 0.005 0.1653 b η	≈ η(around 4 kHz) 0.001
		Cores	100
		Clock Frequency 2.4 GHz
		RAM	80	Gb

Table 5 .

 5 3: Characteristics of the cluster used to compute the VTCR solution in Section 5.4.

	Types	Initial guess VTCR Converged VTCR
	Propagative particular solution	0	0
	Propagative out-of-plane	512	128
	Evanescent out-of-plane	129	65
	Propagative in-plane	128	64

Table 5 .

 5 6 reports the err f for the six accelerometers computed on four frequency bands. Table 5.6: Errors between VTCR and real data accelerations presented in Figures 5.20

	f 1 [%] f 2 [%] f 3 [%]

Simulation of impact response using a frequency approach

See Section 1.1.1.1.

See[Farhat et al. 2003a].Simulation of impact response using a frequency approach

introduced in Section 1.2.2 and described in detail in Section 3.1.Simulation of impact response using a frequency approach

For the interested reader, detailed descriptions of the VTCR versions as well as their pros & cons are provided in Section 1.2.Simulation of impact response using a frequency approach

See Section 4.6 for details about integral computations.Simulation of impact response using a frequency approach

Some angle difficulties prevent to lower this number to 1 sub-domain. Even if a clever algorithm can overcome these problems, its development is out of the aims of the present work.

∇

is the laplacian operator. The term ∇ 2 ∇ 2 w is generated developing ∇ • (∇ • M).Simulation of impact response using a frequency approach

For the sake of simplicity, we suppose that the local orthonormal curvilinear coordinate system α, β, ẑ define also the three symmetry planes of the orthotropic material. In case this assumption is false, this dissertation is still valid but the mathematics becomes more complex. However, the present assumption is always met in the HSS3+ test.

This brief introduction is valid for both time-and frequency-domain representations. Such system versions are examined later on in this Section. Since many letters are already used, every variable used in control theory is denoted by the subscript c to make the notation consistent throughout the whole document without confusing the reader.Simulation of impact response using a frequency approach

Section 3.1.1 provides the definition of straight boundary and polynomial sub-domain.Simulation of impact response using a frequency approach

In case of concentrated masses and springs a term (L [δu i ]) H L [u j ]ds (where is a displacement quantity) arises in the weak form integrals. Hence, the complete 3D matrix would be even bigger (n t × n s × 12) while C ∂Ω e and Q e e∂ Ω e ∩∂ Ω e remain the same.Simulation of impact response using a frequency approach

When they are positioned on a conic structure, they are not rotated of the conic angle.Simulation of impact response using a frequency approach

For the sake of simplicity, only one fixed time window is used.Simulation of impact response using a frequency approach

Simulation of impact response using a frequency approach 5.3), the presented results are for each frequency step. Simulation of impact response using a frequency approach Simulation of impact response using a frequency approach