# Nonlinear vibrations of a ballistic vehicle with friction joints subjected to vibroacoustic excitation - experiments, modeling and simulations

S. Talik, M. Claeys, J.-J. Sinou, H. Valle Canas, J.-P. Lambelin

► To cite this version:

S. Talik, M. Claeys, J.-J. Sinou, H. Valle Canas, J.-P. Lambelin. Nonlinear vibrations of a ballistic vehicle with friction joints subjected to vibroacoustic excitation - experiments, modeling and simulations. Journal of Sound and Vibration, 2024, 588, pp.118526. 10.1016/j.jsv.2024.118526. hal-04642979

HAL Id: hal-04642979

https://hal.science/hal-04642979v1

Submitted on 10 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Nonlinear vibrations of a ballistic vehicle with friction joints subjected to vibroacoustic excitation - experiments, modeling and simulations

S. Talik<sup>a,b</sup>, M. Claeys<sup>a</sup>, J.-J. Sinou<sup>b,c,*</sup>, H. Valle Canas<sup>a</sup>, J.-P. Lambelin<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Nonlinear vibrations

Multipoint correlated random excitation

Industrial structure

## ABSTRACT

This study presents the experiments, modeling and numerical simulations of an industrial aeronautical structure subjected to vibroacoustic excitation. In flight, especially during atmospheric re-entry, an aeronautical structure is exposed to major pressure fluctuations on its external surface. These fluctuations are described as a dynamic excitation which generates structural vibrations of the outer surface and of the components inside the aeronautical structure. Due to some of these internal components such as joints, the assembly may exhibit a nonlinear vibration response. The simulation of these nonlinear vibrations requires reliable modeling of the wall pressure fluctuations and a nonlinear vibration simulation method adapted to the nonlinear modeling of the structure. In this study, a modeling and simulation method is developed to compute the nonlinear vibration response of an industrial assembly to such a surface, random, correlated dynamic excitation. More specifically, numerical simulations are performed by proposing an extension of the well-known Harmonic Balance Method for nonlinear mechanical systems subjected to complex vibroacoustic excitation. The method is validated using a dedicated ground experiment. A metallic industrial assembly representing a ballistic vehicle and including friction joints is used. This structure is subjected to controlled vibroacoustic excitation: diffuse acoustic loading in a reverberant chamber. The structure exhibits a nonlinear vibration response due to friction. The finite element model is validated through an experimental modal analysis with an electrodynamic shaker. The vibroacoustic modeling of the excitation is then validated through a test-simulation comparison using diffuse acoustic field testing at low excitation level. Then the global nonlinear simulation process is validated using diffuse acoustic field testing at increasing excitation levels. Test and simulation results exhibit the same nonlinear behavior: increase of dissipation and softening effect at the main resonances. This work thus represents a new step towards the use of nonlinear vibration simulation methods with real industrial structures and real-life loading.

## 1. Introduction

Many studies in the literature predict the nonlinear behavior of mechanical structures such as academic beams [1–5], jointed structures [6–11] and more or less complex industrial systems [12–15]. In general, studies focus on mechanical systems subjected

<footnote><sup>a</sup> CEA, DAM, CESTA, F-33114 Le Barp, France</footnote>

<footnote><sup>b</sup> Laboratoire de Tribologie et Dynamique des Systèmes UMR CNRS 5513, École Centrale de Lyon, France</footnote>

<footnote><sup>c</sup> Institut Universitaire de France, 75005 Paris, France</footnote>

<footnote>* Corresponding author at: Laboratoire de Tribologie et Dynamique des Systèmes UMR CNRS 5513, École Centrale de Lyon, France.
E-mail address: jean-jacques.sinou@ec-lyon.fr (J.-J. Sinou).</footnote>
to harmonic or random excitation at a single point, the main objective being to measure and identify the non-linear dynamics of a system, or to propose the construction of numerical models capable of reproducing the vibratory behavior observed in simple and controlled experimental tests. It is also worth noting that nowadays modern industrial structures are becoming increasingly complex and require more and more guarantees, both in terms of integrity and operation, from both researchers and engineers. These structures, such as those used in aeronautics or civil engineering, may display a nonlinear behavior that must be taken into account in these guarantee stages before they are put into operation. Therefore, it is important to be able to predict the nonlinear dynamic behavior of a large-scale industrial structure to a so-called real excitation, i.e., encountered during the life cycle of the structure. These excitations can be variable: harmonic, random (or multi-harmonic), shocks, etc. [16]. However, in the particular case of aeronautical structures during an atmospheric re-entry phase, the source of excitation is a correlated multipoint random excitation. Even though some studies have been conducted to examine the nonlinear response of industrial structures to single point excitations, with the industrial tools used today, it is not possible to predict the nonlinear behavior of a large-scale structure (for example with a very large number of degrees of freedom in the finite element model of this structure) to these correlated multipoint random excitations. In a previous publication, the authors showed the importance of considering the correlation of a multipoint random excitation in the vibrational response of a mechanical structure [17]. Different models of a multipoint random excitation exist and can be adapted to many applications, while remaining within the framework of a linear vibratory response [18–20]. However, to the best of our knowledge, no recent study has appeared in the literature on the prediction of a nonlinear vibratory response of an industrial system to a correlated multipoint random excitation. One of the main reasons for this is that current industrial tools do not allow to converge in a suitable time to a nonlinear response of a large-scale structure to a large-scale excitation. The modeling of multipoint correlated random excitations is not adapted to the tools used in industry. Recently, a study was carried out to develop a method to model this type of excitation on a simple academic clamped-clamped beam [21]. In the study proposed, one of the objectives is to adapt this nonlinear strategy for a system comprising a large number of degrees of freedom and subject to a vibroacoustic excitation. The system studied is an industrial model, the Hypervelocity Ballistic 2 (HB2) structure, whose nonlinear vibratory behavior now comes from a friction between two sub-assemblies. Several points represent scientific challenges. First of all, it is necessary to be able to correctly model and implement a vibroacoustic excitation (i.e., a multipoint random excitation) and adapt it to a system with many degrees of freedom. One key point is to considerably reduce the number of stochastic excitation terms, without losing precision in the representation of the real excitation, by projecting the multipoint random excitation on the reduced modal basis of the finite element model of the industrial structure. Secondly it is necessary to adapt the formulation of the non-linear Harmonic Balance Method to solve a nonlinear system subjected to multipoint random excitations. In the following, the effectiveness of producing a reduced model of the vibroacoustic excitation while maintaining an accurate prediction of the nonlinear response of the industrial HB2 structure will be demonstrated. This work can thus be regarded as a new step towards the use of vibration simulation nonlinear methods with real industrial structures and real-life loading. In this context, one of the contributions of this study is to discuss the relevance of both a mechanical modeling strategy and development of specific nonlinear approaches to predict the nonlinear response of a large-scale industrial structure subjected to a real source of vibroacoustic excitation such as correlated multipoint random excitations. In addition, one of the aims of this study is to demonstrate the feasibility of implementing a complete engineering approach based both on experimental tests to understand the dynamic behavior of complex industrial structures subjected to real excitations, and to identify the associated mechanical behavior in order to develop a relevant modeling and numerical simulation strategy.

This paper is organized as follows: first, a complete description is given of the HB2 structure and the proposed experiments for studying a correlated multipoint random excitation, which is an experimental vibroacoustic excitation. The experimental nonlinear vibration response is then observed for increasing levels of the vibroacoustic excitation. The modeling of the structure is then described. The finite element model of the HB2 structure is created, then reduced using a reduction method. A Coulomb friction model is introduced to consider friction in the joints, which is the cause of the nonlinear behavior of the structure. In the next part of this paper, the original modeling of the correlated multipoint random excitation is presented. It uses a projection on the modal basis of the HB2 structure in order to reduce the number of excitation terms to a few degrees of freedom. The final part of this paper concerns the numerical simulation method and the numerical results that are compared to the experimental results. This complete study achieves the objective of predicting the nonlinear vibratory response of a complex industrial structure to a correlated multipoint random excitation as well as the validation of a global numerical strategy combining modeling and simulation based on approximations for modeling complex real excitations, and the consistency of associated non-linear numerical techniques versus vibroacoustic excitations.

## 2. HB2 structure and experiments

### 2.1. *Structure description and diffuse field acoustic testing*

This section details the essential elements of the industrial mechanical structure under consideration, along with the setup used for the experiments with a vibroacoustic excitation. The HB2 structure comes from the external shape of the envelope that has been defined in the scientific literature for aerodynamics studies purposes [22,23]. The structure, representative of a ballistic re-entry vehicle, is an industrial assembly by CEA/CESTA. Fig. 1 shows a picture and a diagram of the HB2 structure. In Fig. 1(a), the model is placed in the acoustic chamber at CEA/CESTA. The reference frame (X, Y, Z) indicated on the image Fig. 1(b) is the global reference frame used for all studies on the HB2 model. The HB2 structure is composed of an envelope, a bottom and a ballast, all made from an aluminum alloy. The outer shell is 1.6 mm thick and is reinforced by circular stiffeners. The flat bottom, on the left
<figure><img src="image_4.png" /><figcaption>Fig. 1. Photo (a) and sectional view of the finite element model (b) of the HB2 structure.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Photo (a) of the internal package with its instrumentation, and diagram (b) of the instrumentation plan.</figcaption></figure>

in Fig. 1(b), is distinct from the rest of the envelope in order to allow insertion of the internal package. The bottom is connected to the envelope by bolted joints.

During assembly, the ballast is inserted into the envelope from the bottom, 23 blades representing 23 friction joints in Fig. 1(b) are compressed during insertion such that the contact is never broken, even during experiments. The HB2 structure has a fundamental transverse mode that stresses the contact joints between the 23 blades and the external envelope and that initiates friction when the structure is subjected to high excitation levels. Most of the mass of the internal package (approximately 25 kg) is concentrated in a solid axisymmetric ballast. A crown of 23 blades serves to center the package in the envelope. A hollow space is left between two blades to allow the passage of instrumentation cables during the experiments. When the package is inserted into the envelope, the blades are deformed and the upper edge of each blade is moved 1 mm inwards. The blade-envelope contact is permanent and the contact force is very low. The value of this force was estimated at 0.1 N using a static finite element simulation of the flexion of the blade with imposed displacement. The Coulomb slip threshold of the blades against the inner wall of the envelope is thus easily reached. At the back, a 1.5 mm thick ferrule is welded to the ballast. The attachment between the internal package and the aeroshell is provided by two flanges assembled to each other by 12 bolted connections.

This industrial structure thus appears as a complex assembly, in which 23 frictional links cause a nonlinear vibratory behavior according to the excitation. To analyze its vibratory behavior, the HB2 structure is instrumented with 21 accelerometers: 15 three-dimensional accelerometers A1 to A15 (with reference Endevco 66M6) and 6 one-dimensional accelerometers A21 to A26 (with reference PCB 352A6). For the sake of readability, only a portion of these accelerometers is photographed in Fig. 2(a) and schematically shown in Fig. 2(b). The set of accelerometers located inside the model are connected to the acquisition setup by means of cables visible in Fig. 1(a) passing through one of the six holes located on the bottom. The origin of the reference frame (X, Y, Z) is located near Accelerometer A4 in the center of the bottom.

The experimental setup used is the reverberant acoustic chamber photographed in Figs. 1(a) and 3. The acoustic sources are two arrays of 12 loudspeakers. The position of the mid-frequency loudspeakers and the natural reverberation properties of the chamber ensure a homogeneous diffuse acoustic field in the room at any point situated more than two meters away from the excitation assemblies. The HB2 structure is suspended at a height of 1.50 m above the ground by means of a winch in order to avoid any acoustic scattering phenomena due to the ground and to ensure free-free mechanical conditions. This vibroacoustic excitation is
<figure><img src="image_6.png" /><figcaption>Fig. 3. Photo of the experimental setup.</figcaption></figure>

recorded by six microphones placed arbitrarily in the acoustic chamber (see Fig. 3). The pressure level measured by each microphone is sent to the acquisition system using conditioners, then the acquisition system calculates the arithmetic mean of the six pressure levels and uses this as the control signal to set the excitation. In this study, the vibroacoustic excitation considered is a white noise, i.e., with a constant Power Spectral Density (PSD), in the frequency band [100–1000] Hz. The excitation is defined with a global RMS level associated with a Sound Pressure Level in decibels (dB SPL). During experiments, the HB2 structure was subjected to increasing levels of excitation from 90 dB to 130 dB in 10 dB steps.

## 2.2. Experimental nonlinear response

This section focuses on the vibration behavior of the HB2 structure solicited by vibroacoustic excitation. For all accelerometers presented in this study, numerical and experimental results are post-processed in order to study the same quantity $a_{RMS}$ which represents the ratio of the RMS acceleration over the RMS acoustic pressure per frequency band. Both numerical and experimental results were post-processed on forty-eighth octave bands, such as

$$a_{RMS}(f_c) = \sqrt{\frac{\int_{f_c^2 - \frac{1}{2N}}^{f_c^2 + \frac{1}{2N}} \text{PSD}_{\text{acceleration}}(f) \, df}{\int_{f_c^2 - \frac{1}{2N}}^{f_c^2 + \frac{1}{2N}} \text{PSD}_{\text{pressure}}(f) \, df}} \quad (1)$$

where $N = 48$ and $f_c$ respectively represent the size of the octave band and the middle frequency of each considered octave band.

The measured quantity $a_{RMS}(f_c)$ is linked to the energy level (in acceleration) of the system divided by the imposed acoustic level (in pressure). With only linear vibroacoustic phenomena, these quantities remain the same whatever the acoustic excitation level. Reading these data as a function of the acoustic level thus allows to correctly observe the nonlinear vibratory phenomena associated with the friction of the blades of the HB2 model.

Fig. 4 shows the processing of the real PSD measurement given by Accelerometer A1 situated at the front end of the envelope for a vibroacoustic excitation level of 90 dB. Due to acoustic resonances present within the reverberant chamber at CEA/CESTA, the superposition of measurements for increasing acoustic levels of excitation does not allow observation of a correct nonlinear behavior, whereas the post-processed data $a_{RMS}$ enables such an observation.

From now on, processed acceleration measurements are analyzed for a limited number of accelerometers. Only the relative accelerations at five acoustic excitation levels (from 90 dB to 130 dB in 10 dB steps) for five accelerometers are presented and analyzed. The five accelerometers, schematized in Fig. 2(b), chosen for the analysis of the vibration measurements are
<figure><img src="image_7.png" /><figcaption>Fig. 4. True experimental data (blue) and post-processed experimental data (red) corresponding to the PSD acceleration given by Accelerometer A1 for a vibroacoustic level of 90 dB. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

* the three-dimensional accelerometer A1 on the external envelope and more precisely on the nose. The measurements associated with this accelerometer are given by dotted curves in Fig. 5;

* the three-dimensional accelerometer A4 at the center of the bottom. The measurements associated with this accelerometer are given by dotted curves in Fig. 6;

* the three-dimensional accelerometer A7 on the front end of the internal package. The measurements associated with this accelerometer are given by dotted curves in Fig. 7;

* the three-dimensional accelerometer A10 on the rear end of the internal package. The measurements associated with this accelerometer are given by dotted curves in Fig. 8;

* the one-dimensional accelerometer A21 on one of the 23 blades of the HB2 structure. The measurements associated with this accelerometer are given by dotted curves in Fig. 9.

For the sake of brevity, only the transverse component (according to the global Z direction) of each acceleration is presented. A zoom into the frequency band [100–200] Hz, where the first vibration mode of the HB2 structure was identified during a modal test (see Section 3), is given for each accelerometer. The acceleration measured along the Z-axis is the most interesting to observe.

Concerning Accelerometer A1, over most of the frequency band, more specifically between 200 Hz and 800 Hz, the curves are all merged, there is no nonlinear vibratory phenomenon observable on this frequency domain. On the other hand, in the frequency band [100–200] Hz, where the first transverse mode of the structure is identified, the dotted curves are distinct. Fig. 5(b) shows very clearly an experimental nonlinear behavior: there is a decrease in the resonance frequency as well as a decrease in the resonance amplitude when the level of acoustic excitation increases. As the excitation level increases, the blades are stressed and rub against the inner part of the aeroshell. This friction has a dual effect on the structure: in stiffness and in dissipation. This friction translates the release of a bonding constraint by the appearance of slippage. It thus leads to a decrease in local stiffness. The local attenuation of stiffnesses leads to a decrease in the resonance frequency of the structure. Moreover, friction is a dissipation phenomenon that contributes to increase the damping of the model. The same nonlinear phenomenon can be observed in Fig. 5(a) for the second transverse vibration mode identified, when the system remains within the framework of a linear behavior, around 925 Hz.

Concerning Accelerometer A4, the acceleration measured along the Z-axis reveals a nonlinear behavior in the two frequency bands [100–200] Hz and [800–1000] Hz: there is a decrease in the resonance frequency along with a decrease in the resonance amplitude as the level of acoustic excitation increases. The level measured for the first transverse mode, located in the frequency band [100–200] Hz and given in Fig. 6(b), remains however much lower than the one measured for Accelerometer A1. This can be explained by the fact that the transverse mode has a lower modal amplitude at the center of the bottom of the structure than on the nose of the envelope. The strong increase in acceleration level that is observed from 300 Hz is due to the presence of vibration modes that are modes of the bottom. These vibroacoustic excitation experiments show that there is little noticeable impact of blade friction on the vibratory behavior of the bottom modes of the HB2 structure.

Concerning the measurements of the accelerometers A7 and A10 located at the front zone and rear zones of the internal package of the HB2 structure, the vibration responses observed along the Z-axis are similar to those described for Accelerometer A1 for the two frequency bands [100–200] Hz and [800–1000] Hz. The nonlinear signature of the HB2 structure, characterized by the friction of the blades against the envelope according to acoustic level, is thus observed on the whole internal package. Concerning the vibration modes of the internal package, identified around 700 Hz and observed in Figs. 7(a) and 8(a), there is no noticeable change in the vibratory behavior of the structure according to acoustic level. These vibroacoustic excitation experiments show that there is little noticeable impact of blade friction on the vibratory behavior of the internal package modes of the HB2 structure.
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Experimental transverse acceleration at Accelerometer A1 on the envelope for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Table 1**</caption>

Values of the ratios of the RMS levels of the acceleration measured by an accelerometer along the Z-axis to the RMS level of the acoustic level (frequency band [100–1000] Hz).

<table><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of experimental RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.434</td><td>1.342</td><td>1.154</td><td>1.023</td><td>0.840</td></tr><tr><td>A4Z</td><td>0.574</td><td>0.558</td><td>0.520</td><td>0.493</td><td>0.457</td></tr><tr><td>A7Z</td><td>0.124</td><td>0.120</td><td>0.113</td><td>0.107</td><td>0.095</td></tr><tr><td>A10Z</td><td>0.189</td><td>0.187</td><td>0.183</td><td>0.177</td><td>0.167</td></tr><tr><td>A21Z</td><td>1.322</td><td>1.307</td><td>1.151</td><td>1.111</td><td>1.064</td></tr></tbody></table>

Finally, the experimental vibratory behavior of the HB2 structure at the location of the blades can be observed via Accelerometer A21 in Fig. 9. This one-dimensional accelerometer measures the acceleration of one of the metal blades along the global Z-axis. It is possible to observe in Fig. 9(a) the two transverse vibration modes respectively located in the two frequency bands [100–200] Hz and [800–1000] Hz. For these two vibration modes, the same experimental nonlinear vibrational phenomenon is observed: there is a decrease in the resonance frequency for the transverse modes as well as a decrease in the resonance amplitude when the level of vibroacoustic excitation increases. The principal mode observable on this accelerometer is a blade mode located between 300 Hz and 400 Hz. There is no noticeable change in the acceleration measured in the frequency band [300–400] Hz according to acoustic level. Another mode is identified in the frequency band [700–800] Hz and this one displays a nonlinear behavior: attenuation of the resonance amplitude.

To summarize of all these experimental results and for the interested reader, Tables 1–3 group together for the studied accelerometers all the ratios of the RMS levels of the acceleration measured along the Z-axis with the RMS level of the pressure according to global vibroacoustic level. The RMS levels are calculated in the frequency band [100–1000] Hz (see Table 1), in the frequency band [100–200] Hz (see Table 2) where the first transverse mode of the HB2 structure is located, and in the frequency band [800–1000] Hz (see Table 3) where the second transverse mode of the HB2 model is identified. Thanks to the data indicated in Tables 1–3, a new quantity is calculated: the coefficient of dissipation of RMS level compared to the RMS level with the referenced acoustic level of 90 dB. This coefficient is defined as the ratio of RMS quotient in a frequency band (i.e., either [100–1000] Hz, or [100–200] Hz or [800–1000] Hz) to the RMS quotient for the 90 dB acoustic level in the same frequency band. With a linear dynamical system, this coefficient would remain equal to 1. These coefficients of experimental dissipation are given for each accelerometer along the Z-axis in Tables 4–6. Their values will be compared with the value of the numerical results given in Section 5.2.2 in order to validate the test-simulation comparisons. The values given in all these tables support the observations on the experimental vibration response of the HB2 structure for the different accelerometers according to acoustic level. For example, Table 5 shows that the experimental response at 130 dB is reduced by a factor of approximately 1/0.36 ≈ 3 in the frequency band [100–200] Hz at the front end of the envelope. Indeed, there is a dissipation of the structure as the blades are stressed to rub against the inner wall of the aeroshell. These vibratory behaviors are observed on all the accelerometers present during the tests: on the envelope (A1), on the bottom (A4), on the internal package (A7, A10) and on one of the 23 blades of the model (A21).

### 3. Modeling of the HB2 structure

Now that the nonlinear behavior of the HB2 structure has been observed experimentally, the paper focuses on the prediction of the nonlinear response of the Finite Element (FE) model of this industrial assembly subjected to vibroacoustic excitation defined
<figure><img src="image_9.png" /><figcaption>Fig. 6. Experimental transverse acceleration at Accelerometer A4 on the bottom for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Experimental transverse acceleration at Accelerometer A7 on the front end of the internal package for different level of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Experimental transverse acceleration at Accelerometer A10 on the rear end of the internal package for different level of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

as a multipoint random excitation. The modeling and simulation method is based on four key points: first a finite element model of the structure is reduced using a dynamic substructuring method. Next, a friction model is introduced for each nonlinear joint in
<figure><img src="image_12.png" /><figcaption>Fig. 9. Experimental transverse acceleration at Accelerometer A21 on one of the 23 blades for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Table 2

Values of the ratios of the RMS levels of the acceleration measured by an accelerometer along the Z-axis to the RMS level of the acoustic level (frequency band [100–200] Hz).

<table><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of experimental RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>0.492</td><td>0.446</td><td>0.345</td><td>0.290</td><td>0.212</td></tr><tr><td>A4Z</td><td>0.014</td><td>0.012</td><td>0.011</td><td>0.010</td><td>0.007</td></tr><tr><td>A7Z</td><td>0.121</td><td>0.111</td><td>0.083</td><td>0.068</td><td>0.046</td></tr><tr><td>A10Z</td><td>0.029</td><td>0.026</td><td>0.021</td><td>0.018</td><td>0.014</td></tr><tr><td>A21Z</td><td>0.104</td><td>0.095</td><td>0.076</td><td>0.067</td><td>0.055</td></tr></tbody></table>

Table 3

Values of the ratios of the RMS levels of the acceleration measured by an accelerometer along the Z-axis to the RMS level of the acoustic level (frequency band [800–1000] Hz).

<table><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of experimental RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>2.973</td><td>2.774</td><td>2.369</td><td>2.083</td><td>1.687</td></tr><tr><td>A4Z</td><td>0.823</td><td>0.770</td><td>0.664</td><td>0.598</td><td>0.512</td></tr><tr><td>A7Z</td><td>0.130</td><td>0.121</td><td>0.103</td><td>0.088</td><td>0.071</td></tr><tr><td>A10Z</td><td>0.204</td><td>0.198</td><td>0.184</td><td>0.174</td><td>0.170</td></tr><tr><td>A21Z</td><td>0.157</td><td>0.144</td><td>0.122</td><td>0.109</td><td>0.093</td></tr></tbody></table>

Table 4

Values of the coefficients of dissipation of the RMS ratio of the acceleration measured by an accelerometer along the Z-axis for each acoustic level (frequency band [100–1000] Hz).

<table><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of experimental dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.94</td><td>0.80</td><td>0.71</td><td>0.59</td></tr><tr><td>A4Z</td><td>1.0</td><td>0.97</td><td>0.91</td><td>0.86</td><td>0.79</td></tr><tr><td>A7Z</td><td>1.0</td><td>0.97</td><td>0.91</td><td>0.87</td><td>0.77</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.99</td><td>0.97</td><td>0.94</td><td>0.88</td></tr><tr><td>A21Z</td><td>1.0</td><td>0.99</td><td>0.87</td><td>0.84</td><td>0.80</td></tr></tbody></table>

this reduced model. The vibroacoustic excitation is projected on the structural modes and the vibration response is simulated in the reduced basis using an adapted Harmonic Balance Method to predict the nonlinear response of the industrial assembly subjected to controlled vibroacoustic excitation. This section details the FE model of the HB2 structure, the reduction method used and the friction modeling introduced in this study in order to simulate the nonlinear response.

The HB2 structure is modeled with the finite element software Abaqus. Figs. 1(b) and 14(a) represent the FE model. Each component of the HB2 structure is defined and meshed separately. The average mesh size used is approximately 10 mm. The complete finite element model of the HB2 model has 171,960 nodes and 128,162 elements. Each component is composed of two types of finite elements with shell solid elements (labeled SC8R) to model the thin walls and volume solid elements (labeled C3D8R) to model the mass parts. Bolted joints are modeled using high rigid connectors. For each joint, a contact surface is defined on each
<caption>**Table 5**</caption>

Values of the coefficients of dissipation of the RMS ratio of the acceleration measured by an accelerometer along the Z-axis for each acoustic level (frequency band [100–200] Hz).

<table><caption>**Table 6**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of experimental dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.91</td><td>0.70</td><td>0.59</td><td>0.36</td></tr><tr><td>A4Z</td><td>1.0</td><td>0.84</td><td>0.75</td><td>0.67</td><td>0.52</td></tr><tr><td>A7Z</td><td>1.0</td><td>0.92</td><td>0.68</td><td>0.56</td><td>0.38</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.90</td><td>0.71</td><td>0.61</td><td>0.49</td></tr><tr><td>A21Z</td><td>1.0</td><td>0.91</td><td>0.73</td><td>0.64</td><td>0.53</td></tr></tbody></table>

Values of the coefficients of dissipation of the RMS ratio of the acceleration measured by an accelerometer along the Z-axis for each acoustic level (frequency band [800–1000] Hz).

<table><caption>**Table 7**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of experimental dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.93</td><td>0.80</td><td>0.70</td><td>0.57</td></tr><tr><td>A4Z</td><td>1.0</td><td>0.94</td><td>0.81</td><td>0.73</td><td>0.63</td></tr><tr><td>A7Z</td><td>1.0</td><td>0.93</td><td>0.79</td><td>0.68</td><td>0.54</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.97</td><td>0.90</td><td>0.85</td><td>0.83</td></tr><tr><td>A21Z</td><td>1.0</td><td>0.92</td><td>0.78</td><td>0.70</td><td>0.59</td></tr></tbody></table>

Numerical and experimental mode pairs identified from the numerical modal analysis of the full finite element model of HB2.

<table><thead><tr><th>Pair</th><th>Mode description</th><th>f<sub>exp</sub> (Hz)</th><th>f<sub>num</sub> (Hz)</th><th>ε<sub>f</sub> (%)</th><th>MAC (%)</th></tr></thead><tbody><tr><td>1</td><td>Transverse mode</td><td>130.9</td><td>136.2</td><td>4.1</td><td>98.1</td></tr><tr><td>2</td><td>Transverse mode</td><td>135.4</td><td>139.4</td><td>2.89</td><td>98.0</td></tr><tr><td>3</td><td>Blades mode</td><td>336.2</td><td>365.2</td><td>8.60</td><td>89.7</td></tr><tr><td>4</td><td>Envelope mode (2 lobes)</td><td>344.7</td><td>408.7</td><td>18.57</td><td>94.5</td></tr><tr><td>5</td><td>Envelope mode (3 lobes)</td><td>474.8</td><td>502.9</td><td>5.91</td><td>90.7</td></tr><tr><td>6</td><td>Envelope mode (3 lobes)</td><td>478.4</td><td>504.8</td><td>5.51</td><td>90.3</td></tr><tr><td>7</td><td>Envelope mode (4 lobes)</td><td>713.7</td><td>685.5</td><td>3.96</td><td>81.7</td></tr><tr><td>8</td><td>Transverse mode</td><td>925.0</td><td>905.9</td><td>2.06</td><td>86.2</td></tr><tr><td>9</td><td>Envelope mode (5 lobes)</td><td>1033.6</td><td>992.5</td><td>3.97</td><td>90.9</td></tr></tbody></table>

side of the joint and a kinematic coupling to a dedicated control node is defined. The joint is then defined as a rigid connector between those two control nodes. Contact between each blade and the outer shell is realized with connectors. On the blade part, the reference nodes are denoted “Ref\_blade\_1” (for blade 1 for example). On the aeroshell part, the reference nodes are denoted “Ref\_contact\_zone\_1” (for the envelope nodes associated with blade 1 for example). The connectors representing each link have a very high stiffness (10¹⁰ N m⁻¹) in the directions r\_blade\_1, θ\_blade\_1 and z\_blade\_1 as well as in rotation around z\_blade\_1 (10¹⁰ N m rad⁻¹) for each blade. The rotations around r\_blade\_1 and around θ\_blade\_1 are left free for each blade.

Once the finite element model has been constructed, it must be validated before the vibratory responses can be studied for various dynamic environments. The validation criterion used in this study is the classical comparison of the numerical modal analysis with the experimental modal analysis. The FE model was validated through an experimental modal analysis with an electrodynamic shaker and a laser vibrometer. Photos of this additional experiment are given in Fig. 10. The chosen excitation is a low-level random force excitation applied at mid-height of the envelope. This input signal is defined by a PSD of constant level (0.10 N RMS) along the frequency band [100–1000] Hz. The accelerometer locations and scanning laser vibrometry measurement points are joined to form a united measurement mesh. The Frequency Response Functions (FRF) are measured for each node of this measurement mesh. Finally, the modal frequencies, damping and shapes of the HB2 structure were identified using the PolyMAX [24] optimization algorithm based on all the measured FRFs. Twenty modes of vibration are experimentally identified along the frequency band [100–1000] Hz. Using FEMTOOLS software, nine pairs of experimental and numerical modes were identified on the basis of a MAC criterion greater than 80%. Table 7 indicates the mode pairs identified from both numerical and experimental modal analyses, within the notations f\_exp, f\_num, ε\_f and MAC respectively indicate the experimental frequency, the numerical frequency, the frequency error and the MAC value of the paired modes. As indicated in Table 7, nine pairs of vibration modes were identified in the frequency band of interest [100–1000] Hz, validating the finite element model of the HB2 structure. As only one blade out of 23 is equipped with an accelerometer, only one blade mode is identified between 300 Hz and 400 Hz, whereas 23 blade modes were computed by numerical modal analysis. The mode pair shapes are illustrated in Figs. 11–13. These figures show the wealth of experimental data obtained by using a 3D scanning laser vibrometer, which produces modal deformations very similar to those obtained numerically.

The global nonlinear vibration response of such a model to vibroacoustic excitation is impossible with the classical simulation method due to the large number of degrees of freedom subjected to dynamic air pressure. A reduction method is thus employed to predict the nonlinear response of the HB2 structure on only a few physical degrees of freedom. The reduction method used in this study is a hybrid method developed by Dassault Systèmes and implemented on Abaqus [25]. This method keeps the static modes
<figure><img src="image_13.png" /><figcaption>Fig. 10. Photos of the experimental set-up for 3D scanning laser vibrometry measurements of the HB2 structure excited at one point by a suspension-mounted electrodynamic shaker.</figcaption></figure>

and the eigenmodes with the chosen degrees of freedom of interest. The static modes are the static responses of the structure to a unit displacement imposed on one chosen degree of freedom and all the others are blocked. This method is analogous to the method of Craig and Bampton [26], the difference comes from the dynamic modes that can be computed with hybrid interfaces: for the computation of these modes, the chosen physical degrees of freedom can be blocked or free, whereas in the Craig and Bampton method, they must be blocked. In this study, this reduction method is used on the HB2 structure in order to keep the following 136 degrees of freedom

* 19 × 3 degrees of freedom for all the degrees of freedom of each accelerometer present on the structure;

* 23 degrees of freedom corresponding to the 23 friction joints between the blades and the external envelope in order to implement a nonlinear numerical method;

* 56 degrees of freedom corresponding to the 56 first vibration modes below 1000 Hz computed from the complete FE model of the HB2 structure, without any additional boundary conditions on the other chosen degrees of freedom.

Fig. 14(b) shows the finite element model of one of the 23 blades present on the HB2 structure. For each blade of the crown, a local cylindrical frame is given and reference nodes are defined in order to numerically link the crown to the aeroshell (observed in blue in Fig. 14(b)). In this study, the 23 chosen degrees of freedom corresponding to the 23 friction joints are those in translation along the $\theta_{blade}$ direction in the local cylindrical frame of each blade. The nonlinear forces are applied to these degrees of freedom in the reduced base. The modeling of the nonlinear links is chosen to be necessary and just sufficient in order to achieve the vibration prediction objectives. It thus takes into account these nonlinearities while ensuring an easy implementation of these contributions in the numerical resolution for the prediction of the nonlinear vibratory response (for more details, see the previous work by Claeys et al. [13]).

Concerning the friction modeling between the blades and the envelope, the friction model used is the Jenkins model, which is represented in Fig. 15 with stiffness connector $k_{bound}$ and Coulomb force $F_{Coulomb}$. Note that, in this study, the parameters $k_{bound}$ and $F_{Coulomb}$ are the same for all blades. For each temporal step $t_i$, and for each blade, a virtual force $\dot{f}_{blade}(t_i)$ is evaluated and compared to the Coulomb force $F_{Coulomb}$ [27]

$$ \dot{f}_{blade}(t_i) = f_{blade}(t_{i-1}) + k_{bound}(u_{blade}(t_i) - u_{blade}(t_{i-1})) \quad (2) $$

Two states are defined: if the virtual force at temporal step $t_i$ is greater than the Coulomb force, then the considered blade is in a slipping state. If not, the considered blade is tied to the envelope, such that

$$ f_{blade}(t_i) = \begin{cases} \dot{f}_{blade}(t_i) & \text{if } |\dot{f}_{blade}(t_i)| < F_{Coulomb} \quad (\text{stick}) \\ \text{sgn}(\dot{f}_{blade}(t_i)) F_{Coulomb} & \text{if } |\dot{f}_{blade}(t_i)| \ge F_{Coulomb} \quad (\text{slip}) \end{cases} \quad (3) $$

In the case of the HB2 structure, this Jenkins model is applied to each friction joint in each local frame of each blade. In this study, each translation along each orthoradial direction is kept in the reduced model. This allows us to define a nonlinear force $F_{NL}(x)$, which depends on the reduced degrees of freedom vector $x$, such that

$$ F_{NL}(x) = \sum_{\text{blade}} (k_{\text{bound}} u_{\text{blade}} - f(u_{\text{blade}})) V_{\text{blade}} \quad (4) $$

where $u_{blade} = V_{blade}^T x$ and $V_{blade}$ defines the connector vector of the considered blade [13].
<figure><img src="image_14.png" /><figcaption>Fig. 11. Pairing of experimental (red) and numerical (blue) mode shapes in YX (a-c-d) and YZ (b-d-f) frames: pair 1 (a-b), pair 2 (c-d) and pair 3 (e-f). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Pairing of experimental (red) and numerical (blue) mode shapes in YX (a-c-d) and YZ (b-d-f) frames: pair 4 (a-b), pair 5 (c-d) and pair 6 (e-f). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 13. Pairing of experimental (red) and numerical (blue) mode shapes in YX (a–c–d) and YZ (b–d–f) frames: pair 7 (a–b), pair 8 (c–d) and pair 9 (e–f). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 14. The FE model of the HB2 structure (a) with some of the accelerometers (yellow crosses) and (b) a zoom on the friction joint between the envelope and the ballast defined in the local frame of the considered blade. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

For the interested reader it is worth remembering that various friction models can be used to simulate the restoring forces in mechanical contacts via friction effects. In the present study, the choice of the Jenkins model is the result of previous works [13,15] showing that this one-parameter model is sufficient to provide a good representation of the experimental response of the mechanical system under study. More specifically the hysteresis phenomenon with the possibility of adhesion or sliding phases at the frictional interface, depending of the evolution of the force in the friction joint and the dynamic behavior of the system are taken into account in this model. One of the advantage of the Jenkins model is that the Coulomb threshold is the unique parameter to be identified. In the present study, the identification of this Coulomb threshold is undertaken by matching simulations with experimental results at 130 dB (corresponding to the highest excitation level proposed in this study). An improvement on this friction model could be to include the possibility of partial slip or micro-slip (e.g. using the Dahl model [28], the Iwan model [29,30], the Valanis model [31] and the LuGre friction model [32] for example).

The nonlinear problem associated with the reduced finite element model of the structure HB2 subjected to vibroacoustic excitation $F_{excit}$ is therefore defined by

$$M\ddot{x} + D\dot{x} + Kx = F_{excit} + F_{NL}(x) \quad (5)$$

where M, D and K respectively define the mass, damping and stiffness matrices in the reduced basis of degrees of freedom of the HB2 structure. $F_{NL}(x)$ is defined with the Jenkins friction model as described by Eq. (4). Thanks to the experimental modal analysis, a modal damping per vibration mode is identified. A Rayleigh model is then introduced to define the damping such that $D = \alpha M + \beta K$. Parameters $\alpha$ and $\beta$ are identified from the modal damping obtained experimentally on the transverse modes located respectively at 131 Hz and 925 Hz, which gives the following parameters: $\alpha = 11.9 \text{ s}^{-1}$ and $\beta = 1.0 \text{ s}^{-7}$ s.

From now on, the vibroacoustic excitation vector $F_{excit}$ needs to be computed in the reduced basis with the complete information on the correlation of the random multipoint excitation applied to the outer surface of the aeroshell, and the complete information on the level of this vibroacoustic excitation.

# 4. Modeling of the vibroacoustic excitation

## 4.1. Coupled FEM-BEM-FEA model

Far from the structure, the acoustic field in the reverberant chamber can be modeled as a diffuse field, band-limited ([100 - 1000] Hz) white noise. Yet the pressure fluctuations at the structure surface are strongly influenced by acoustic scattering phenomena. A vibroacoustic model is developed using Wave6 software. The structural Finite Element Model (FEM) previously presented (without the nonlinear friction joints) is coupled to a boundary element model of the external air volume (BEM). The two models are illustrated in Fig. 16. The numerical modal basis of the FEM, which was presented and validated previously, was computed with Abaqus, then imported into Wave6. A modal damping is then defined. Each modal damping coefficient was chosen to be perfectly equivalent to the Rayleigh model presented previously.

In order to model the surrounding air volume, a closed external surface of the structure is defined and meshed as illustrated in the lower image of Fig. 16. The external volume of this surface is defined as the BEM subsystem. A volume diffuse acoustic field excitation is then introduced. The excitation is implemented using the reciprocity relationship between direct field radiation and
<figure><img src="image_18.png" /><figcaption>Fig. 15. A representation (a) and a force-displacement diagram (b) of the Jenkins friction model.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. From top to bottom, sectional views of structural FEM mesh and BEM mesh.</figcaption></figure>

diffuse reverberant loading [33]. Then a strong coupling is defined between the FEM and the BEM models, on the external surface of the envelope.

In a previous study [34], where only the external envelope of HB2 structure was considered, both internal and external air volumes were modeled using a single BEM model. When adding the internal package, the complex geometry of the internal air volume makes its modeling by BEM almost impossible. An acoustic Finite Element Model would be well-suited for this internal air volume. For the purpose of this research study however, the vibroacoustic modeling has been simplified as much as possible. This study focuses on the structural non-linear behavior. The vibroacoustic model is only introduced in order to compute an accurate dynamic excitation applied to this model. To that end, the BEM model of the external air volume is necessary because acoustic scattering has a major effect on this dynamic excitation. Modeling the internal air volume enables us to compute internal acoustic resonances. It is assumed that these resonances have little effect on the dynamic response of the structure. The modal projection method presented in next section is restricted to structural resonances, even though an extension to acoustic resonances could be developed. Thus, in this work, the internal air volume was not modeled.

The vibroacoustic response of the system is computed for all structural modes at each frequency spanning the frequency band [100–1000] Hz. This response can then be computed for any FEM node by modal superposition. Fig. 17 presents the PSD
<figure><img src="image_20.png" /><figcaption>Fig. 17. Post-processing of the experimental data (blue), of the simulated output data with the complete model of the vibroacoustic excitation (green) and of the simulated output data with the proposed reduced vibroacoustic excitation (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

acceleration simulated at the location of Accelerometer A1 on the envelope (green curve) compared to the measured data of this sensor (blue curve). The corresponding vibroacoustic excitation level is 90 dB. The qualitative comparison of the simulation results with the experiments is excellent. For all accelerometer locations, in Fig. 2(b), the simulated results are close to the experimental data, with a maximum 6 dB difference on the RMS value (in the [100–1000] Hz frequency band). This vibroacoustic modeling, where the friction nonlinear behavior is neglected, is preliminary to the present study. It is the underlying linear model that is relevant at low excitation levels, as illustrated in Fig. 17. These vibroacoustic simulations are performed including all the structural modes, yet as it will be useful in the next section, it is also possible to perform these vibroacoustic simulations by selecting only one or two vibration modes of interest.

## 4.2. Modal projection of the vibroacoustic excitation

This section covers an original method used to define the vibroacoustic excitation in a reduced finite element model. The method for modeling a multipoint correlated random excitation for a academic clamped-clamped beam was developed in a dedicated article [21]. In the following, the modeling of a vibroacoustic excitation in a reduced finite element model will be discussed. One originality is to propose an approximate and controlled modeling of vibroacoustic excitations in order to ensure coherence with the Harmonic Balance Method presented in Section 5.1.

Considering the structural FEM, the surface acoustic loading is a stationary multipoint random excitation, applied to each node of the model connected to the air. Such loading can be characterized using a cross spectral density matrix $S(\omega)$ (where $\omega$ denotes the pulsation), which is defined for the $N_e$ points submitted to the loading by:

$$ S(\omega) = \begin{bmatrix} S_{11}(\omega) & S_{12}(\omega) & \dots & S_{1N_e}(\omega) \\ S_{21}(\omega) & S_{22}(\omega) & \dots & S_{2N_e}(\omega) \\ \vdots & \vdots & \ddots & \vdots \\ S_{N_e,1}(\omega) & S_{N_e,2}(\omega) & \dots & S_{N_e N_e}(\omega) \end{bmatrix} \qquad (6) $$

where the function $S_{ij}(\omega)$ represents the Power Spectral Density (PSD) of the temporal excitation $f_j(t)$, defined to the degree of freedom $i$. The function $S_{ij}(\omega)$ represents the Cross Spectral Density (CSD) between the temporal excitations $f_j(t)$ and $f_j(t)$, which corresponds to the multipoint temporal excitation defined at the degree of freedom $j \neq i$. Thanks to the Wiener-Khintchine relation, the $S(\omega)$ matrix is defined positive [19].

The main idea behind this method is to project the spectral density matrix $S(\omega)$ on the basis of vibration modes $[Y_1^{c.m.}, \dots, Y_{N_m}^{c.m.}]$ (where $N_m$ corresponds to the number of vibration modes) of the complete finite element model of the HB2 structure. These vibration modes must be normalized with respect to the mass matrix of the structure. The vector $F_{excit}^{c.m.}(t)$, representing the acoustic excitation, is written in the vibration mode basis of the HB2 structure

$$ F_{excit}^{c.m.}(t) = \sum_{k=1}^{N_m} y_k(t) Y_k^{c.m.} \qquad (7) $$

The modal variables $y_k(t)$ are sought in this section. The vectors $Y_k^{c.m.}$ correspond to the vibration modes of the complete finite element model considered.
The projection of the spectral density matrix $S(\omega)$ onto the vibration mode basis $[\mathbf{Y}_1^{\text{c.m.}}, \dots, \mathbf{Y}_N^{\text{c.m.}}]$ leads to a new matrix $S^{\text{proj}}(\omega)$ defined by

$$S^{\text{proj}}(\omega) = \mathbf{P}^T \mathcal{S}(\omega) \mathbf{P} \qquad (8)$$

where $\mathbf{P}^T$ denotes the transpose operation. The matrices $\mathcal{S}(\omega)$ and $\mathbf{P}$ are described by

$$\mathcal{S}(\omega) = \begin{bmatrix} \mathbf{S}(\omega) & \mathbf{S}(\omega) & \dots & \mathbf{S}(\omega) \\ \mathbf{S}(\omega) & \mathbf{S}(\omega) & \dots & \mathbf{S}(\omega) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{S}(\omega) & \mathbf{S}(\omega) & \dots & \mathbf{S}(\omega) \end{bmatrix} \qquad (9)$$

$$\mathbf{P} = \begin{bmatrix} \mathbf{Y}_1^{\text{c.m.}} & 0 & 0 \\ 0 & \mathbf{Y}_2^{\text{c.m.}} & \vdots \\ \vdots & \vdots & 0 \\ 0 & \dots & 0 \end{bmatrix} \mathbf{Y}_N^{\text{c.m.}} \qquad (10)$$

The components of matrix $S^{\text{proj}}$ are

$$S_{kl}^{\text{proj}}(\omega) = (\mathbf{Y}_k^{\text{c.m.}})^T S(\omega) (\mathbf{Y}_l^{\text{c.m.}}) \qquad (11)$$

These components are also the excitation terms used in the modal vibroacoustic simulations that are performed using the Wave6 solver. In these linear simulations, $U_n$, the PSD of the displacement of a node $n$ of the structural FE model, is computed using modal superposition

$$U_n(\omega) = \sum_k \sum_l \mathcal{H}_k^*(\omega) S_{kl}^{\text{proj}}(\omega) \mathcal{H}_l(\omega) \mathbf{Y}_k^{\text{c.m.}}|_n \mathbf{Y}_l^{\text{c.m.}}|_n \qquad (12)$$

where $\mathcal{H}_k$ is the modal transfer function

$$\mathcal{H}_k(\omega) = (\Omega_k^2 - 2j\xi_k \Omega_k \omega - \omega^2)^{-1} \qquad (13)$$

$\Omega_k$ is the proper pulsation of mode $k$ that is known and $\xi_k$ is the modal damping that is chosen. When performing the linear vibroacoustic simulations with only one mode $k$, the double sum in Eq. (12) reduces to a unique term, leading to a direct relation between $S_{kl}^{\text{proj}}$ and $U_n$ (when computed with only mode $k$)

$$S_{kk}^{\text{proj}}(\omega) = \frac{U_n(\omega)}{|\mathcal{H}_k(\omega)|^2 (\mathbf{Y}_k^{\text{c.m.}}|_n)^2} \qquad (14)$$

Then, once all the modal acceptances $S_{kk}^{\text{proj}}$ are computed, the cross-terms $S_{kk}^{\text{proj}}$ can be computed using simulations involving two modes $k$ and $l$. When $U_n$ is computed with only modes $k$ and $l$, Eq. (12) lead to

$$S_{kl}^{\text{proj}}(\omega) = \frac{U_n(\omega) - S_{kk}^{\text{proj}}(\omega) (\mathbf{Y}_k^{\text{c.m.}}|_n)^2 - S_{ll}^{\text{proj}}(\omega) (\mathbf{Y}_l^{\text{c.m.}}|_n)^2}{(\mathcal{H}_k^*(\omega)\mathcal{H}_l(\omega) + \mathcal{H}_l^*(\omega)\mathcal{H}_k(\omega)) (\mathbf{Y}_k^{\text{c.m.}}|_n \mathbf{Y}_l^{\text{c.m.}}|_n)} \qquad (15)$$

Thus, performing several linear computations with each mode $k$ and each mode pair $k$ and $l$, all the terms of the matrix $S^{\text{proj}}(\omega)$ are computed. These terms, which represent the acoustic excitation of the structure, take into account the diffuse acoustic field and acoustic scattering effects thanks to the coupled FEM/BEM model. It can then be shown that the matrix $S^{\text{proj}}(\omega)$ is positive semi-definite [21]. The Cholesky decomposition is used as follows

$$S^{\text{proj}}(\omega) = \mathbf{G}(\omega)\mathbf{G}^*(\omega)^T \qquad (16)$$

where $\mathbf{G}(\omega)$ is a lower triangular matrix

$$\mathbf{G}(\omega) = \begin{bmatrix} G_{11}(\omega) & 0 & \dots & 0 \\ G_{21}(\omega) & G_{22}(\omega) & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ G_{N_m 1}(\omega) & G_{N_m 2}(\omega) & \dots & G_{N_m N_m}(\omega) \end{bmatrix} \qquad (17)$$

The $G_{ij}(\omega)$ coefficients are determined with the principal minors of matrix $S^{\text{proj}}(\omega)$. It should be noted that this method is general for any type of multipoint random excitation. The temporal modal excitation $(y_k(t))_{\forall k \in \llbracket 1, N_m \rrbracket}$ is then given by

$$y_k(t) = \sqrt{2\Delta\omega} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |G_{kl}(\omega_n)| \cos(\omega_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) \qquad (18)$$

The parameters present in the previous Eq. (18) are defined in such a way that
* p is the finite length of the series of cosine functions of the approximated excitation. It should be considered as a power of 2 in order to use Fast-Fourier Transforms (FFT) [35];

* $\Delta\omega$ defines the pulsation step along the frequency band $[\omega_{min}, \omega_{max}]$. It is given by $\Delta\omega = \frac{\omega_{max}-\omega_{min}}{p}$;

* $\omega_n$ corresponds to the frequency discretization [19]. It is given by $\omega_n = \omega_{min} + (n - \frac{1}{2}) \Delta\omega$. For the numerical simulation, FFT algorithm discretization is chosen in order to use the Fast-Fourier Transform technique [35];

* $\theta_{kl}(\omega_n)$ is defined as

$$ \theta_{kl}(\omega_n) = \tan^{-1} \left( \frac{\mathfrak{F}(G_{kl}(\omega_n))}{\mathfrak{R}(G_{kl}(\omega_n))} \right) \qquad (19) $$

where $\mathfrak{F}$ and $\mathfrak{R}$ are the imaginary and real parts, respectively;

* $\phi_{ln}$ are independent random phases uniformly distributed between 0 and $2\pi$. The random part of the acoustic excitation comes from these variables.

The cosine series decomposition is related to the use of the Harmonic Balance Method, presented in Section 5.1. Thanks to Eq. (7), the description of the multipoint random excitation vector in the full basis of the finite element model considered given by the vector $F_{excit}^{c.m.}(t)$ can be written

$$ F_{excit}^{c.m.}(t) = \sqrt{2\Delta\omega} \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |G_{kl}(\omega_n)| \cos(\omega_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) Y_k^{c.m.} \qquad (20) $$

The reduction method used on Abaqus ensures that the vibration modes of the reduced model are exactly equal to those of the complete model projected in the reduced basis. Therefore, the vector of the multipoint random excitation in the reduced basis $F_{excit}(t)$ can be written

$$ F_{excit}(t) = \sqrt{2\Delta\omega} \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |G_{kl}(\omega_n)| \cos(\omega_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) Y_k \qquad (21) $$

where the vectors $Y_k$ correspond to the vibration modes of the reduced finite element model considered.

The vibroacoustic modeling of the excitation was validated in a preliminary study through a test-simulation comparison using diffuse acoustic field testing at low excitation level (90 dB), when the nonlinearities are supposed to be small in the overall response of the structure. This validation is illustrated for one accelerometer in Fig. 17. Three curves are plotted in this figure, the blue curve corresponds to the post-processing of the experimental PSD of acceleration for Accelerometer A1, the green curve and the red curve correspond to the post-processing of the numerical PSD of acceleration of the degree of freedom closest to Accelerometer A1. The green curve comes from the linear vibroacoustic simulation presented in Section 4.1. During this simulation performed using Wave6 software, the complete acoustic excitation defined by Eq. (20) is computed on all the structural nodes connected to air. The modal acceptance terms (i.e., the $G_{ij}(\omega)$ coefficients of the $G(\omega)$ matrix) are calculated over the whole frequency band [100–1000] Hz. The red curve is then calculated using the nonlinear simulation method presented in Section 5.1 where the acoustic loading is computed using Eq. (21). Due to calculation time issues, the modal acceptance terms $G_{kk}$ are only calculated at frequencies corresponding to the resonance frequencies $\Omega_k$ of the vibration modes of the numerical model of the HB2 structure. The cross-modal acceptance terms $G_{kl}$ with $k \neq l$ are computed at the median pulsation between $\Omega_k$ and $\Omega_{k}$. This results in the following approximation

$$ G_{kl}(\omega_n) = G_{kl} \left( \frac{\Omega_k + \Omega_l}{2} \right) \quad \forall n \in [1, p] \qquad (22) $$

Note that information on the vibroacoustic excitation level for all frequencies in the frequency band [0–1000] Hz is not complete, as the cross-modal acceptances terms are identified at specific frequencies. This will result in an underestimation of the excitation, and therefore somewhat underestimated results from the reduced model compared to the results from the full model.

Therefore, Fig. 17 highlights two important facts:

* the comparison between the blue curve and the green curve is very satisfactory, it validates the construction of the FE model of this industrial structure with the construction of each component and the contact between each sub-assembly. The slight difference is interpreted as being due to the acoustic excitation term which is modeled as a homogeneous acoustic diffuse field whereas experimentally it is influenced by the acoustic modes of the reverberant chamber;

* the comparison between the green curve and the red curve is also interesting, it indicates that the modeling of a reduced multipoint random excitation is consistent with the whole model. Both curves are merged at frequencies corresponding to the resonance frequencies of the vibration modes, and the output acceleration of the structure is underestimated at frequencies far from the vicinity of these resonance frequencies. As this study focuses on the prediction of a nonlinear phenomenon which only occurs at resonance frequencies, this result is satisfactory and validates the construction of the reduced vibroacoustic excitation.

It should be noted that the second resonance peak in the frequency band [800–1000] Hz is lower for the green curve (coming from the complete linear vibroacoustic simulation). In fact, this is due to a loss of energy with acoustic emission taken into account in the Wave6 software.
## 5. Numerical simulation

Henceforth, the nonlinear force $F_{NL}(x)$ and the vibroacoustic excitation $F_{excit}$ in the reduced element model of the HB2 structure define the nonlinear problem described by Eq. (5), which can be solved using the Harmonic Balance Method (HBM) adapted to multipoint random excitations [17].

### 5.1. Harmonic balance method adapted to random excitations

For this purpose, this nonlinear method must be adapted for the use of acoustic excitations seen as multipoint correlated random excitations. The HBM seeks the response of the nonlinear system defined in Eq. (5) as a truncated Fourier series (if this solution exists), such that

$$x(t) = \mathbf{B}_0 + \sum_{n=1}^{p} (\mathbf{A}_n \sin(n\Omega t) + \mathbf{B}_n \cos(n\Omega t)) \quad (23)$$

where $p$ is the chosen order of the truncated Fourier series and $(\mathbf{B}_0, (\mathbf{A}_n, \mathbf{B}_n)_{\forall n \in \llbracket 1, p \rrbracket})$ are the unknown Fourier coefficients of solution $x(t)$ to be determined. As discussed in [17,36], the HBM can be adapted to random excitations by choosing the frequency resolution $\Delta f$ as the fundamental frequency of the excitation.

Thus in the case of vibroacoustic excitation, the multipoint random excitation can be foreseen as a multipoint deterministic excitation with one fundamental pulsation $\Omega = 2\pi\Delta f$ and using Eq. (21), vector force $F_{excit}(t)$ can be defined by a finite Fourier series of order $p$, such that

$$F_{\text{excit}}(t) = \sum_{n=1}^{p} (S_{n,\text{excit}} \sin(n\Omega t) + C_{n,\text{excit}} \cos(n\Omega t)) \quad (24)$$

$$S_{n,\text{excit}}(t) = -\sqrt{2\Delta\omega} \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} |G_{kl}(\omega_n)| \sin(\theta_{kl}(\omega_n) + \phi_{ln}) \right) Y_k \quad (25)$$

$$C_{n,\text{excit}}(t) = \sqrt{2\Delta\omega} \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} |G_{kl}(\omega_n)| \cos(\theta_{kl}(\omega_n) + \phi_{ln}) \right) Y_k \quad (26)$$

Moreover, it is assumed that nonlinear contributions $F_{NL}(t)$ can be solved in the finite Fourier series of order $p$

$$F_{\text{NL}}(t) = C_0 + \sum_{n=1}^{p} (S_n \sin(n\Omega t) + C_n \cos(n\Omega t)) \quad (27)$$

where $(C_0, (S_n, C_n)_{\forall n \in \llbracket 1, p \rrbracket})$ are the Fourier coefficients of the nonlinear force $F_{NL}(t)$. One of the consequences of the choice of $\Delta f$ as the fundamental frequency is that the number of harmonics $p$ chosen in the approximate nonlinear solution (23) can be significant, which can constitute a major disadvantage in terms of computation time and data storage issues. In order to determine the value of the Fourier coefficients $(\mathbf{B}_0, (\mathbf{A}_n, \mathbf{B}_n)_{\forall n \in \llbracket 1, p \rrbracket})$, decompositions (23), (24) and (27) are re-injected into Eq. (5). This leads to a set of $N_m(2p+1)$ nonlinear equations given by

$$\mathbf{K}\mathbf{B}_0 = \mathbf{C}_0 \quad (28)$$

$$\begin{bmatrix} \mathbf{K} - (n\Omega)^2 \mathbf{M} & -n\Omega \mathbf{D} \\ n\Omega \mathbf{D} & \mathbf{K} - (n\Omega)^2 \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{A}_n \\ \mathbf{B}_n \end{bmatrix} = \begin{bmatrix} \mathbf{S}_{n,\text{excit}} \\ \mathbf{C}_{n,\text{excit}} \end{bmatrix} + \begin{bmatrix} \mathbf{S}_n \\ \mathbf{C}_n \end{bmatrix} \quad \forall n \in \llbracket 1, p \rrbracket \quad (29)$$

Coefficients $(C_0, (S_n, C_n)_{\forall n \in \llbracket 1, p \rrbracket})$ depend on coefficients $(\mathbf{B}_0, (\mathbf{A}_n, \mathbf{B}_n)_{\forall n \in \llbracket 1, p \rrbracket})$. An extension of the classic Alternate Frequency-Time domain method (AFT-method [37]) is used to calculate these Fourier coefficients $(C_0, (S_n, C_n)_{\forall n \in \llbracket 1, p \rrbracket})$ (see [36,38] for more details) thanks to Eq. (4).

Finally, the nonlinear Eqs. (28) and (29) are solved by minimizing the following function

$$H(\mathbf{X}, \Omega) = \mathbf{A}\mathbf{X} - \mathbf{B} - \mathbf{B}_{\text{NL}}(\mathbf{X}) \quad (30)$$

with

$$\mathbf{A} = \text{Diag} \left( \mathbf{K}, \begin{bmatrix} \mathbf{K} - (n\Omega)^2 \mathbf{M} & -n\Omega \mathbf{D} \\ n\Omega \mathbf{D} & \mathbf{K} - (n\Omega)^2 \mathbf{M} \end{bmatrix}_{\forall n \in \llbracket 1, p \rrbracket} \right) \quad (31)$$

$$\mathbf{B} = [0 \quad \mathbf{S}_{1,\text{excit}} \quad \mathbf{C}_{1,\text{excit}} \quad \dots \quad \mathbf{S}_{p,\text{excit}} \quad \mathbf{C}_{p,\text{excit}}]^T \quad (32)$$

$$\mathbf{B}_{\text{NL}} = [\mathbf{C}_0 \quad \mathbf{S}_1 \quad \mathbf{C}_1 \quad \dots \quad \mathbf{S}_p \quad \mathbf{C}_p]^T \quad (33)$$
<figure><img src="image_21.png" /><figcaption>Fig. 18. Experimental (dotted lines) and numerical (solid lines) transverse acceleration at Accelerometer A1 on the envelope for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Minimizing Eq. (30) imposes a convergence criterion described as

$$ \frac{\|\mathbf{H}(\mathbf{X}, \mathbf{\Omega})\|_2}{\|\mathbf{B}\|_2} < \epsilon_{HBM} \qquad (34) $$

where $\epsilon_{HBM}$ is a chosen numerical precision and $\|\cdot\|_2$ is the quadratic norm. For the rest of the study, $\epsilon_{HBM}$ is chosen to be equal to $10^{-3}$.

Note that the proposed extension of the Harmonic Balance Method to vibroacoustic excitation is based on a redefining of the fundamental pulsation, as well as approximation of a multipoint correlated random excitation by reducing the number of excitation terms while remaining in a formalism of the HBM applied to random excitations.

## 5.2. Numerical results

In this section, the simulation results are compared with the experimental results of the test campaign on the HB2 structure subjected to acoustic excitation. This comparison is performed on the 5 accelerometers whose measurements were presented in Section 2. The parameter of the friction model used to perform this set of simulations is the same for all excitation levels: $F_{Coulomb} = 0.05$ N. In order to appreciate the order of magnitude of this parameter, a static calculation with imposed displacement on an isolated metal blade is performed on Abaqus. The radial force applied to a blade is estimated at 0.1 N. The order of magnitude of the Coulomb threshold thus seems reasonable.

In a first step, the numerical results are analyzed alone in order to identify the numerical nonlinear vibration signature for some vibration modes of the HB2 structure associated with the friction model chosen in Section 3. In a second step, these results are presented for all the degrees of freedom located near the accelerometers in order to be compared to the measurements of the experimental campaign performed with an acoustic excitation. This comparison is performed through a graphical superposition of the experimental and numerical results processed by octave bands (see Section 2.2). The comparison of the numerical and experimental dissipation coefficients by acoustic level also allows us to conclude on the relevance of the nonlinear simulation method developed. As a reminder, these coefficients are introduced in Section 2.2 to quantify the nonlinear phenomenon observed experimentally.

### 5.2.1. Study of the numerical nonlinear vibration signature

The first vibration mode of the HB2 structure in the study is the first transverse mode experimentally identified at approximately 131 Hz (see Table 7). First, the numerical study focuses on the degree of freedom corresponding to the Z component of the closest node from the experimental location of the first accelerometer A1 located on the nose of the envelope of the HB2 structure. The results from the numerical simulations are given by the solid lines in Fig. 18 for the different acoustic levels. Concerning the first transverse mode of the HB2 structure located in the frequency band [100–200] Hz, a distinction of the curves according to the acoustic level is observed numerically (see Fig. 18(b)), which is representative of a nonlinear vibratory signature. This is the same nonlinear vibration signature as that observed experimentally: there is a decrease in amplitude and a decrease in resonance frequency according to acoustic level. This phenomenon is very well-represented for the first transverse mode at Accelerometer A1 along the global Z-axis. The four other degrees of freedom chosen for this numerical study are the ones corresponding respectively to locations closest to Accelerometers A4, A7, A10 and A21 along the Z-direction. The results from the numerical simulations are respectively given by the solid lines in Figs. 19–22 for the different acoustic levels. The same nonlinear vibrational phenomenon is observed on the first transverse mode, i.e., there is a decrease in resonance amplitude and a decrease in resonance frequency with acoustic level.

The second vibration mode studied in this section is the second transverse mode, experimentally identified at approximately 925 Hz (see Table 7). For all degrees of freedom respectively corresponding to Accelerometers A1, A4, A7 and A10 along the
<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Experimental (dotted lines) and numerical (solid lines) transverse acceleration at Accelerometer A4 on the bottom for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Z-direction, neither the resonance amplitude, nor the resonance frequency of the second transverse mode located in the frequency band [800–1000] Hz are strongly modified by the acoustic level: the nonlinear signature is not well-represented numerically for this vibration mode.

For other vibration modes present in the numerical modal analysis of the FE model of the HB2 structure, the analysis of the numerical results varies according to the chosen degree of freedom. For example, for the degree of freedom corresponding to Accelerometer A1, it is possible to observe in Fig. 18(a) that the solid curves are merged between 300 Hz and 850 Hz, where certain resonances are represented. There are three resonances between 300 Hz and 400 Hz and three resonances between 580 Hz and 700 Hz. No nonlinear vibratory phenomena are observable at the respective frequencies of these modes. However, an additional resonance at 830 Hz is observed in Fig. 20(a) for the degree of freedom corresponding to Accelerometer A7. The amplitude of this resonance is completely attenuated as the acoustic level increases. This mode is very pronounced at the front end of the inner package compared to the resonance observed in Figs. 18(a) and 19(a) for Accelerometers A1 and A4. This mode corresponds to the first longitudinal mode, it is not observable for the nodes located on the envelope corresponding to the accelerometers A1 and A4. The fourth accelerometer is Accelerometer A10 located at the rear end of the inner package of the HB2 structure. The results from the numerical simulations given by the solid lines in Fig. 21 for the different acoustic levels are consistent with the numerical results given for Accelerometer A7. The nonlinear signature is observed for the mode previously mentioned at 830 Hz. The second transverse mode located around 925 Hz is weakly resonant along the Z-axis. Indeed, the transverse mode observed on Abaqus has a relatively low modal amplitude at the location of the node associated with Accelerometer A10. Finally, the internal package modes and some blade modes, located between 350 Hz and 850 Hz, resonate very strongly according to the acoustic level for Accelerometer A21 located on one of the 23 blades of the HB2 structure (see Fig. 22(a)). A nonlinear signature is notably observed at 500 Hz and 710 Hz for which the resonance concerned is completely attenuated at 130 dB. These frequencies correspond to coupled blade and envelope modes according to the numerical modal analysis.

To summarize these observations of the numerical results and for the interested reader, Tables 8–10 group together, for the accelerometers studied, all the ratios of the RMS levels of the simulated acceleration along a global axis to the RMS level of the sound level according to global sound level. The RMS levels are calculated in the frequency band [100–1000] Hz (see Table 8), in the frequency band [100–200] Hz (see Table 9) where the first transverse mode of the HB2 model is located, and in the frequency band [800–1000] Hz (see Table 10) where the second transverse mode of the HB2 model is located. Using the data given in these tables, a coefficient of dissipation introduced in Section 2.2 is calculated in the different frequency bands: either [100–1000] Hz (see Table 11), or [100–200] Hz (see Table 12) or [800–1000] Hz (see Table 13). The values shown in all these tables support the observed results on the simulated vibration response of the HB2 structure for these different degrees of freedom according to acoustic level. Indeed, there is a dissipation of the structure as the blades are stressed by the acoustic level. The numerical nonlinear signature of the HB2 structure is only well-predicted by simulation in the frequency band [100–200] Hz. Comparison with the experimental results and a discussion are given in the following section.

### 5.2.2. Experiments-simulations comparison

Now that the analysis of the numerical results is given for different of the HB2 structure vibration modes, the comparison with the experiments is carried out. As a reminder, the superposition of numerical (solid lines) and experimental (dotted lines) results is given in Figs. 18–22.

Concerning Accelerometer A1, the simulated acceleration level is much lower than the one measured between 200 Hz and 800 Hz. This is mainly due to the choice of the calculation of the cross-modal acceptance terms for the acoustic excitation given in Section 4. Indeed, for this study, the choice was to insist on the representation of the acoustic level for the resonance frequencies
<figure><img src="image_23.png" /><figcaption>Fig. 20. Experimental (dotted lines) and numerical (solid lines) transverse acceleration at Accelerometer A7 on the front end of the internal package for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21. Experimental (dotted lines) and numerical (solid lines) transverse acceleration at Accelerometer A10 on the rear end of the internal package for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Experimental (dotted lines) and numerical (solid lines) transverse acceleration at Accelerometer A21 on one of the 23 blades for different levels of vibroacoustic excitation: 90 dB (black), 100 dB (blue), 110 dB (green), 120 dB (orange) and 130 dB (red) - (a) in the frequency band [100–1000] Hz and (b) zoom in the frequency band [100–200] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<caption>**Table 8**</caption>

Values of the ratios of the RMS levels of the simulated acceleration (at the location of the degree of freedom corresponding to an accelerometer) along a global axis to the RMS level of the sound level according to global sound level (frequency band [100–1000] Hz).

<table><caption>**Table 9**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of numerical RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.418</td><td>1.415</td><td>1.319</td><td>1.177</td><td>1.299</td></tr><tr><td>A4Z</td><td>0.651</td><td>0.652</td><td>0.634</td><td>0.596</td><td>0.616</td></tr><tr><td>A7Z</td><td>0.165</td><td>0.166</td><td>0.162</td><td>0.164</td><td>0.173</td></tr><tr><td>A10Z</td><td>0.263</td><td>0.262</td><td>0.249</td><td>0.231</td><td>0.233</td></tr><tr><td>A21Z</td><td>0.695</td><td>0.702</td><td>0.678</td><td>0.682</td><td>0.652</td></tr></tbody></table>

Values of the ratios of the RMS levels of the simulated acceleration (at the location of the degree of freedom corresponding to an accelerometer) along a global axis to the RMS level of the sound level according to global sound level (frequency band [100–200] Hz).

<table><caption>**Table 10**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of numerical RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>0.435</td><td>0.432</td><td>0.151</td><td>0.092</td><td>0.158</td></tr><tr><td>A4Z</td><td>0.013</td><td>0.013</td><td>0.005</td><td>0.003</td><td>0.004</td></tr><tr><td>A7Z</td><td>0.101</td><td>0.100</td><td>0.036</td><td>0.026</td><td>0.044</td></tr><tr><td>A10Z</td><td>0.023</td><td>0.023</td><td>0.008</td><td>0.006</td><td>0.011</td></tr><tr><td>A21Z</td><td>0.103</td><td>0.102</td><td>0.039</td><td>0.018</td><td>0.014</td></tr></tbody></table>

Values of the ratios of the RMS levels of the simulated acceleration (at the location of the degree of freedom corresponding to an accelerometer) along a global axis to the RMS level of the sound level according to global sound level (frequency band [800–1000] Hz).

<table><caption>**Table 11**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Ratio of numerical RMS levels ((m s<sup>-2</sup>)/Pa)</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>2.981</td><td>2.975</td><td>2.785</td><td>2.482</td><td>2.741</td></tr><tr><td>A4Z</td><td>0.831</td><td>0.828</td><td>0.778</td><td>0.710</td><td>0.787</td></tr><tr><td>A7Z</td><td>0.111</td><td>0.111</td><td>0.113</td><td>0.149</td><td>0.183</td></tr><tr><td>A10Z</td><td>0.212</td><td>0.200</td><td>0.150</td><td>0.109</td><td>0.135</td></tr><tr><td>A21Z</td><td>0.839</td><td>0.850</td><td>0.799</td><td>0.829</td><td>0.713</td></tr></tbody></table>

Values of the coefficients of dissipation of the RMS ratio of the simulated acceleration of a degree of freedom near an accelerometer along the Z-axis for each acoustic level (frequency band [100–1000] Hz).

<table><caption>**Table 12**</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of numerical dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.99</td><td>0.93</td><td>0.83</td><td>0.92</td></tr><tr><td>A4Z</td><td>1.0</td><td>1.00</td><td>0.97</td><td>0.92</td><td>0.95</td></tr><tr><td>A7Z</td><td>1.0</td><td>1.00</td><td>0.98</td><td>0.99</td><td>1.05</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.99</td><td>0.95</td><td>0.88</td><td>0.88</td></tr><tr><td>A21Z</td><td>1.0</td><td>1.01</td><td>0.98</td><td>0.98</td><td>0.94</td></tr></tbody></table>

Values of the coefficients of dissipation of the RMS ratio of the simulated acceleration of a degree of freedom near an accelerometer along the Z-axis for each acoustic level (frequency band [100–200] Hz).

<table><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of numerical dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.99</td><td>0.35</td><td>0.21</td><td>0.36</td></tr><tr><td>A4Z</td><td>1.0</td><td>0.99</td><td>0.36</td><td>0.22</td><td>0.30</td></tr><tr><td>A7Z</td><td>1.0</td><td>0.99</td><td>0.36</td><td>0.25</td><td>0.44</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.99</td><td>0.36</td><td>0.26</td><td>0.47</td></tr><tr><td>A21Z</td><td>1.0</td><td>0.99</td><td>0.38</td><td>0.17</td><td>0.14</td></tr></tbody></table>

of the HB2 structure vibration modes, while neglecting the information on the excitation level for frequencies located outside the scope of these resonance frequencies for reasons of calculation time. This explains the fact that the comparison of experimental and numerical acceleration levels is good at the resonance frequencies of the transverse modes of the HB2 model (see Fig. 18(a)), but also the fact that the acceleration level simulated decreases sharply away from these resonances. This is particularly the case
<table><caption>**Table 13**
Values of the coefficients of dissipation of the RMS ratio of the simulated acceleration of a degree of freedom near an accelerometer along the Z-axis for each acoustic level (frequency band [800–1000] Hz).</caption><thead><tr><th rowspan="2">Accelerometer</th><th colspan="5">Coefficient of numerical dissipation</th></tr><tr><th>90 dB</th><th>100 dB</th><th>110 dB</th><th>120 dB</th><th>130 dB</th></tr></thead><tbody><tr><td>A1Z</td><td>1.0</td><td>0.99</td><td>0.93</td><td>0.83</td><td>0.92</td></tr><tr><td>A4Z</td><td>1.0</td><td>0.99</td><td>0.94</td><td>0.85</td><td>0.95</td></tr><tr><td>A7Z</td><td>1.0</td><td>0.99</td><td>1.01</td><td>1.34</td><td>1.64</td></tr><tr><td>A10Z</td><td>1.0</td><td>0.95</td><td>0.71</td><td>0.52</td><td>0.64</td></tr><tr><td>A21Z</td><td>1.0</td><td>1.01</td><td>0.95</td><td>0.98</td><td>0.85</td></tr></tbody></table>

from 180 Hz along the Z-axis for the degree of freedom near Accelerometer A1. The nonlinear experimental and numerical vibration signatures remain coherent: there is a decrease in resonance amplitude and a decrease in resonance frequency for the first transverse mode of the HB2 structure. Numerically, this behavior does not differ between 90 dB and 100 dB, then a very strong dissipation of the resonance amplitude from 100 dB appears and finally a decrease of the resonance frequency is observed from the acoustic level of 120 dB. Experimentally, the two phenomena occur simultaneously according to the acoustic level, which is not represented numerically. This may be due to the choice of friction model used, with an overly low Coulomb threshold coefficient. Between 120 dB and 130 dB the simulated resonance frequency remains low and the modal amplification increases. This simulated evolution is also related to the choice of a Jenkins model, similar trends have been observed and explained by Al Sayed et al. [39]. The Jenkins model has two linear asymptotic states: at very low excitation level, the blades remain stuck to the envelope, at very high excitation level, the blade are almost always sliding on the envelope surface which is the second asymptotic state. The maximal dissipation is reached when stick and slip phases are of equal duration. Experimentally this increase of amplification at high level is not observed, revealing another limit of the Jenkins model in this application. However, the comparison remains very satisfactory in terms of the objective of representing the nonlinear vibratory phenomenon for the first transverse mode along the Z direction. Indeed, as given in Table 12 at the 130 dB acoustic level, a coefficient of numerical dissipation of 0.36 is reached numerically which is exactly the same coefficient of experimental dissipation shown in Table 5. The comparison of the numerical and experimental coefficients of dissipation provides new information on the ability of the model to represent and predict the nonlinear signature of the HB2 structure. Overall, this comparison for Accelerometer A1 is satisfactory.

Concerning the comparison for Accelerometer A4, contrary to Accelerometer A1, the simulated acoustic level is well-represented since the simulated acceleration levels along the Z-axis axes are of the same order of magnitude as the measurements, with the same trends: a progressive increase in acceleration from 300 Hz. Yet, the representation of the nonlinear vibration signature is very satisfactory along the Z-axis for the first transverse mode of the model.

Concerning the comparison for the accelerometers A7 and A10, located on either side of the internal package of the HB2 structure, the analysis of the superposition of the data set is similar to the comparisons for Accelerometers A1 and A4 for the 2 main transverse modes. However, the comparison of the experimental and numerical dissipation coefficients on the frequency band [100–200] Hz, where the first transverse mode of the model is located (see Table 12), is very satisfactory. For Accelerometer A7, a coefficient of numerical dissipation of 0.44 is obtained at 130 dB, very close to the coefficient of experimental dissipation of 0.38. For Accelerometer A10, a coefficient of numerical dissipation of 0.47 is obtained at 130 dB, very close to the coefficient of experimental dissipation of 0.49. These comparisons for Accelerometers A7 and A10 are thus satisfactory.

Concerning Accelerometer A21 located on one of the 23 blades of the HB2 structure, the superposition of all the data is shown in Fig. 22. Contrary to what is analyzed for the other accelerometers in the paper, the first transverse mode along the Z-axis for Accelerometer A21 is completely attenuated at 130 dB. The vibration phenomenon is thus overestimated for this mode as indicated by the coefficients of numerical dissipation between 100 Hz and 200 Hz (see Table 12). Two other vibration modes, corresponding to coupled blade and envelope modes, located numerically at 500 Hz and 710 Hz, are completely attenuated for the numerical results and they are also attenuated for the experimental results. Indeed, the resonances located experimentally at 480 Hz and 720 Hz are also completely attenuated from 130 dB. This vibratory behavior is thus very well-reproduced by the simulation and does not only concern the transverse modes of the HB2 model. Finally, the comparison deteriorates from 750 Hz since the modes of the blades resonate more strongly in the simulation results than experimentally. The acoustic behavior of the internal air volume must have a significant effect on the vibroacoustic response of the blades. Thus, the vibroacoustic modeling of this internal air volume and its coupling with the nonlinear modeling presented in this article would be a promising way to reduce those test-simulation differences.

In conclusion, in light of the elements presented in this section, the nonlinear vibration simulation of the HB2 model solicited by an acoustic excitation is satisfactory on many points: it is possible to account for a nonlinear vibration phenomenon at different locations of the industrial structure according to different degrees of freedom. The numerical results presented are consistent with the experiments with notably very good comparisons in the frequency band [100–200] Hz where the first transverse mode of the HB2 structure is located. However, some other numerical results are not in complete agreement with the measurements of the test campaign performed. In particular, the nonlinear vibratory phenomenon on the second transverse mode of the model is not sufficiently well-represented compared to the experimental data. The limitations in test-simulation comparison must be related to the choices made for the nonlinear joint model: the choice of friction model given in Section 3, the choice of degrees of freedom in translation along the direction θ only and the choice of the Coulomb threshold. However, all these choices enabled to predict a complex nonlinear vibratory phenomenon which is found on several vibration modes (at 139 Hz, 500 Hz and 710 Hz) of an
industrial structure. This meets one of the objectives of this paper, which was to propose a controlled, accurate and sufficient model to reproduce certain nonlinear phenomena. The overall satisfaction with the simulation results and the associated experiments-simulations confrontations comes from the very good correlation on the level of simulated excitation thanks to an original modeling of a multipoint random excitation on a reduced number of degrees of freedom, and on the other hand on the representation of a complex nonlinear vibratory phenomenon on an industrial structure according to the level of this simulated excitation.

## 6. Conclusion

In this paper, the work proposed is an application of an original approach to modeling a multipoint correlated random excitation in an industrial context. This application is dedicated to the use of a simulation method of nonlinear vibrations with an extension of the Harmonic Balance Method. First, the structure studied (the HB2 structure) was designed to enhance nonlinear vibration thanks to friction joints in order to be representative of an aerospace structure submitted to a vibroacoustic load, representative of atmospheric re-entry. In order to experimentally study such a nonlinear behavior, an original setup, a reverberation chamber with loudspeakers arrays, was developed at the CEA/CESTA site. A nonlinear vibration response was then obtained experimentally and may now be confronted to numerical results.

Thus, this paper developed the complete modeling of both the structure and the excitation. The finite element model of the HB2 structure is reduced thanks to a reduction method implemented in the Abaqus software, allowing few degrees of freedom to be kept. The number of degrees of freedom retained is just sufficient and necessary for the observation of the nonlinear response and for the introduction of a nonlinear force. This nonlinear force, modeled in its simplest form with a Jenkins model, is applied to only one degree of freedom per friction joint in the HB2 structure. Once the reduction method used, the modeling of acoustic excitation must be stored in a reduced excitation vector. In order to achieve such an objective, an original method previously published by the authors was used. This consists in projecting the spectral density matrix onto the vibration modes of the HB2 structure. The reduced vibroacoustic excitation is therefore modal and it was validated through a test-simulation comparison using diffuse acoustic field testing at low excitation level. Finally, using an adapted Harmonic Balance Method, the nonlinear response of the HB2 structure was predicted. Numerical results are consistent with experiments in the vicinity of the vibration modes, mainly the first transverse mode, of the assembly.

## CRediT authorship contribution statement

S. Talik: Writing – original draft, Visualization, Validation, Software, Investigation, Data curation. M. Claeys: Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. J.-J. Sinou: Writing – review & editing, Supervision, Methodology, Conceptualization. H. Valle Canas: Writing – review & editing, Resources, Data curation. J.-P. Lambelin: Writing – review & editing, Validation, Supervision, Methodology, Investigation, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The data that has been used is confidential.

## Acknowledgments

The authors gratefully acknowledge the support of the CEA/CESTA, France teams that helped with this study. They acknowledge in particular the experimental service that provided the measures presented in this paper. J.-J. Sinou acknowledges the support of the Institut Universitaire de France, France.

## References

[1] T. Anderson, A. Nayfeh, B. Balachandran, Experimental verification of the importance of the nonlinear curvature in the response of a cantilever beam, J. Vib. Acoust. 118 (1) (1996) 21-27.

[2] J. Fang, I. Elishakoff, R. Caimi, Nonlinear response of a beam under stationary random excitation by improved stochastic linearization method, Appl. Math. Model. 19 (2) (1995) 106-111.

[3] J. Huang, R. Su, Y. Lee, S. Chen, Nonlinear vibration of a curved beam under uniform base harmonic excitation with quadratic and cubic nonlinearities, J. Sound Vib. 330 (21) (2011) 5151-5164.

[4] A. Nayfeh, Nonlinear transverse vibrations of beams with properties that vary along the length, J. Acoust. Soc. Am. 53 (3) (1973) 766-770.

[5] J. Ray, C. Bert, Nonlinear vibrations of a beam with pinned ends, J. Eng. Ind. 91 (4) (1969) 997-1004.

[6] H. Festjens, G. Chevallier, J. Dion, A numerical tool for the design of assembled structures under dynamic loads, Int. J. Mech. Sci. 75 (2013) 170-177.

[7] M. Eriten, M. Kurt, G. Luo, D. McFarland, L. Bergman, A. Vakakis, Nonlinear system identification of frictional effects in a beam with a bolted joint connection, Mech. Syst. Signal Process. 39 (1-2) (2013) 245-264.

[8] H. Festjens, G. Chevallier, J. Dion, Nonlinear model order reduction of jointed structures for dynamic analysis, J. Sound Vib. 333 (7) (2014) 2100-2113.
[9] J. Armand, L. Salles, C. Schwingshackl, D. Süß, K. Willner, On the effects of roughness on the nonlinear dynamics of a bolted joint: A multiscale analysis, Eur. J. Mech. A Solids 70 (2018) 44–57.

[10] R. Lacayo, L. Pesaresi, J. Groß, D. Fochler, J. Armand, L. Salles, C. Schwingshackl, M. Allen, M. Brake, Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver, Mech. Syst. Signal Process. 114 (2019) 413–438.

[11] M. Jin, G. Kosova, M. Cenedese, W. Chen, A. Singh, D. Jana, M. Brake, C. Schwingshackl, S. Nagarajaiah, K. Moore, J.-P. Noël, Measurement and identification of the nonlinear dynamics of a jointed structure using full-field data; part ii - nonlinear system identification, Mech. Syst. Signal Process. 166 (2022) 108402.

[12] T. Detroux, L. Renson, L. Masset, G. Kerschen, The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems, Comput. Methods Appl. Mech. Engrg. 296 (2015) 18–38.

[13] M. Claeys, J.-J. Sinou, J.-P. Lambelin, R. Todeschini, Modal interactions due to friction in the nonlinear vibration response of the harmony test structure: Experiments and simulations, J. Sound Vib. 376 (2016) 131–148.

[14] D. Roettgen, M. Allen, Nonlinear characterization of a bolted, industrial structure using a modal framework, Mech. Syst. Signal Process. 84 (2017) 152–170.

[15] T. Roncen, J.-J. Sinou, J.-P. Lambelin, Experiments and simulations of the structure harmony-gamma subjected to broadband random vibrations, Mech. Syst. Signal Process. 159 (2021) 107849.

[16] G. Kerschen, K. Worden, A. Vakakis, J.-C. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mech. Syst. Signal Process. 20 (3) (2006) 505–592.

[17] S. Talik, J.-J. Sinou, M. Claeys, J.-P. Lambelin, Nonlinear vibrations of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations - experiments, modeling and simulations, Commun. Nonlinear Sci. 110 (2022) 1–27.

[18] M. Shinozuka, Simulation of multivariate and multidimensional random processes, J. Acoust. Soc. Am. 49 (1B) (1971) 357–368.

[19] M. Shinozuka, C.-M. Jan, Digital simulation of random processes and its applications, J. Sound Vib. 25 (1) (1972) 111–128.

[20] F. Poirion, C. Soize, Numerical simulation of homogeneous and inhomogeneous Gaussian stochastic vector fields, La Rec. Aérosp. 1989 (1989).

[21] S. Talik, M. Claeys, J.-J. Sinou, J.-P. Lambelin, An efficient approach for predicting the nonlinear vibrations of a beam system subjected to multipoint correlated random excitation, Eur. J. Mech. A Solids 11 (2022).

[22] J. Gray, Summary Report on Aerodynamic Characteristics of Standard Models HB-1 and HB-2, Tech. Rep. AEDC-TDR-64-137, Arnold Engineering Development Center, 1964.

[23] D. Damjanovic, B. Rasuo, S. Mandic, S. Vukovic, J. Isakovic, Usability of the comparative experimental-numerical supersonic test cases with the HB reference model, in: Proceedings of the 29th Congress of the International Council of the Aeronautical Sciences, 2014.

[24] B. Peeters, H.V. der Auwaerer, P. Guillaume, J. Leuridan, The polymax frequency-domain method: a new standard for modal parameter estimation? Shock Vib. 11 (3–4) (2004) 395–409.

[25] M. Kim, V. Belsky, M. Belyi, Substructure generation using automated multilevel substructuring, 2013, patent US 2013/0124150A1 and EP 2597578A1.

[26] R. Craig, M. Bampton, Coupling of substructure for dynamic analysis, AIAA J. 1 (2) (1968) 1313–1319.

[27] J. Guillen, C. Pierre, An efficient, hybrid, frequency-time domain method for the dynamics of large-scale dry-friction damped structural systems, in: IUTAM Symposium on Unilateral Multibody Contacts, Springer, 1999, pp. 169–178.

[28] P. Dahl, Solid friction damping of mechanical vibrations, J. Guid. Control Dynam. 14 (1976) 1675–1682.

[29] W. Iwan, A distributed-element model for hysteresis and its steady-state dynamic response, J. Appl. Mech. 33 (1966) 893–900.

[30] B. Deaner, M. Allen, M. Starr, D. Segalman, H. Sumali, Application of viscous and iwan modal damping models to experimental measurements from bolted structures, J. Vib. Acoust. 117 (27) (2015) 1–12.

[31] L. Gaul, J. Lenz, Nonlinear dynamics of structures assembled by bolted joints, Acta Mech. 125 (1–4) (1997) 169–181.

[32] L. Gaul, R. Nitsche, The role of friction in mechanical joints, Appl. Mech. Rev. 54 (2001) 93–106.

[33] P.J. Shorter, R.S. Langley, On the reciprocity relationship between direct field radiation and diffuse reverberant loading, J. Acoust. Soc. Am. 117 (2005) 85–95.

[34] M. Claeys, H. Valle-Canas, B. Alcoverro, Full-frequency vibroacoustic modeling of a ballistic re-entry aeroshell and validation through diffuse field acoustic testing, Appl. Sci. 12 (11) (2022).

[35] W. Vetterling, B. Flannery, W. Press, S. Teukolsy, Numerical Recipes in Fortran 77, Cambridge University Press, 1992.

[36] T. Roncen, J.-P. Lambelin, J.-J. Sinou, Nonlinear vibrations of a beam with non-ideal boundary conditions and stochastic excitations - experiments, modeling and simulations, Commun. Nonlinear Sci. 74 (2019) 14–29.

[37] T. Cameron, J. Griffin, An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems, J. Appl. Mech. 56 (1) (1989) 149–154.

[38] T. Roncen, J.-J. Sinou, J.-P. Lambelin, Experiments and nonlinear simulations of a rubber isolator subjected to harmonic and random vibrations, J. Sound Vib. 451 (2019) 71–83.

[39] B.A. Sayed, E. Chatelet, S. Baguet, G. Jacquet-Richardet, Dissipated energy and boundary condition effects associated to dry friction on the dynamics of vibrating structures, Mech. Mach. Theory 46 (4) (2011) 479–491.