# Nonlinear vibrations of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations - experiments, modeling and simulations

S Talik, Jean-Jacques Sinou, M Claeys, J-P Lambelin

► To cite this version:

S Talik, Jean-Jacques 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. Communications in Nonlinear Science and Numerical Simulation, 2022, 110, pp.106328. 10.1016/j.cnsns.2022.106328. hal-03604341

HAL Id: hal-03604341

https://hal.science/hal-03604341

Submitted on 10 Mar 2022

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>

Research paper

# Nonlinear vibrations of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations - experiments, modeling and simulations

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

Check for
updates

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

## ARTICLE INFO

**Article history:**

Received 15 February 2021

Received in revised form 24 January 2022

Accepted 5 February 2022

Available online 11 February 2022

**Keywords:**

Nonlinear vibrations

Broadband random excitations

Numerical simulation

Nonlinear modeling

## ABSTRACT

This study presents experiments and numerical simulations of a beam with non-ideal boundary conditions and subjected to two broadband correlated or uncorrelated random excitations. Modeling of a beam with non-ideal boundary conditions and additional static pretension, as well as the modeling of multipoint correlated or uncorrelated random excitations, are developed. Based on numerical simulations via an extension of the Harmonic Balance Method, comparisons between experiments and numerical simulations are performed for various broadband bi-point correlated random excitations. In order to achieve such an objective, experiments are performed for four types of correlation between the two broadband random excitations: correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations. Experimental results demonstrate that the different configurations of broadband random excitations drastically affect, not only the amplification or attenuation of some symmetric or anti-symmetric modes of the beam in the vicinity of their primary resonance, but also the appearance of harmonics and a combination of harmonic components. Good correlations are observed between experiments and numerical simulations, thus validating the proposed modeling and the computational strategy for the prediction of the nonlinear vibrational phenomena of the beam system subjected to two correlated or uncorrelated broadband random excitations.

To be noted that the experimental data set including the input and output measurements of the four configurations are provided as supplementary files on Talik et al. [1].

© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ([http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)).

## 1. Introduction

The integration of nonlinearities and the prediction of nonlinear vibrations in mechanical systems are today common topics among researchers and engineers. They allow them to investigate more complex problems in which nonlinearity leads to a variety of complex phenomena. This is particularly the case for the vibration analysis of beams that remains

<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>
an important issue in structural engineering applications in the field of civil engineering or aerospace. The nonlinear vibrations of these structures often lead to unanticipated vibrational behavior that results in severe structural damage. In addition the boundary conditions present in real systems do not correspond to the ideal conditions assumed in the usual mathematical modeling. These small deviations from those ideal conditions can lead to changes in the nonlinear behavior of such systems.

A large number of papers have been devoted to the analysis of beam vibrations. The choice of types of beams (i.e., linear or nonlinear, uniform and non-uniform as well as planar or non-planar) induces a more or less complex modeling [2]. This can also lead to the observation of complex phenomena that involve not only the occurrence of nonlinear dynamic behavior [3] but also potential transitions to chaos for flexible beams [4–10]. In addition, this has a direct impact on the resolution and analysis methods to be used in order to obtain reliable and validated results. A brief state-of-the-art of recent directions in modeling and vibration analysis of beam can be found in [9]. Another challenging part is to be able to predict the nonlinear dynamic behavior of beams with non-ideal boundary conditions. The proposed study will focus more specifically on this second topic. Such a subject of study has received particular attention and remains a challenge due to its inherent difficulties [11–16]. One of the most popular approach to predict the characteristics of beam systems with non-ideal boundary conditions is the perturbation methods, such as multiple scales [2,17]. Even if one of the advantages of the method of multiple scales is to provide an analytical solution of the nonlinear problem, one restriction of such approach is based on the consideration of "small nonlinearities". Thus numerical methods such as the iteration perturbation method [13], the Harmonic Balance Method [18–20], the shooting method [21–23] or the combination of the asymptotic numerical method with the Harmonic Balance Method [24] have been recently proposed to predict harmonically forced vibrations of beam systems.

Although extensive work has been developed to examine the nonlinear response of beam systems with non-ideal boundary conditions subjected to harmonic, multi-harmonic or parametric excitations [19,25,26], little research has been done to determine the effect of different excitation signals on beam systems and the effectiveness of modeling and numerical methods to accurately reproduce these more realistic vibration cases. Although random excitations are one of the most popular excitation techniques for experimental modal analysis, the modeling of broadband excitations and its implementation in numerical methods remains challenging. Indeed simplifying assumptions on excitations are common in computer-aided mechanical design for engineering applications. They make the problem to be treated easier to solve mathematically and sometimes easier to model and incorporate into the system equations. Recently, some numerical studies have investigated the nonlinear response of a beam under one random excitation [15,27]. Roncen et al. [23] studied both experimentally and theoretically the nonlinear characteristics of a beam subjected to random excitation, but the proposed study focused more on the efficiency of numerical methods, such as the Harmonic Balance Method and the shooting method, for robust prediction of the nonlinear behavior of a beam, subjected to one basic random excitation, in an adequate computing time. To the best of our knowledge, there is little work in the literature offering a complete vision combining experimentation, modeling and numerical simulation to validate a global strategy for constructing nonlinear model allowing to predict the nonlinear dynamic behavior of a system subjected to random multipoint excitations. Indeed the consequences of different multipoint solicitations on the nonlinear dynamic behavior of a mechanical system and the potential relationship between the excitation sources and the evolution of the nonlinear signature of the vibrational response are areas of expertise that remain little explored through numerical simulations during the design process of nonlinear mechanical systems. This is mainly due to the complexity of the nonlinear methods to be developed to achieve such an objective. This situation is particularly detrimental since many efforts have been made to enhance vibration testing [28] and to improve the identification and characterization of nonlinearities from experiments [29,30]. So one of the main scientific challenges is not only to better understand the evolution of the nonlinear dynamic behavior of mechanical systems, but also to be able to correctly predict the nonlinear responses of mechanical systems subjected to various complex solicitations. From this point of view, taking into account multiple correlated or uncorrelated random excitations corresponds to solicitations often encountered in practice for many nonlinear dynamic systems. Having a good confidence in modeling through confrontation between experiments and simulations for such complex excitations is nowadays a crucial point. A second major interest for the engineer is to ultimately have a robust and efficient nonlinear method to predict the nonlinear dynamic behavior and performance of a mechanical system in more realistic scenarios. The following study is intended to meet this objective and to better investigate the nonlinear response that exhibits many interesting characteristics when correlated and uncorrelated random excitations are applied on the mechanical system of interest. One of the objectives of this study is to answer this question by proposing a complete approach allowing the development of numerical techniques and the construction of nonlinear models capable of reproducing nonlinear vibrational responses based on different solicitations. The efficiency of the proposed methodology will be studied on the nonlinear response of the clamped-clamped beam subjected to two correlated or uncorrelated broadband random excitations. This work is continuation of the previous work by Claeys et al. [19] and by Roncen et al. [23], who studied the nonlinear behavior of the beam system to swept sine and a basis random excitation, respectively. One of the main original points of the proposed study in comparison with the latter is to undertake a complete experimental protocol to analyze the nonlinear behavior of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations. In addition, a proposal for an appropriate modeling not only of the beam system but also of a multipoint correlated random excitation is discussed in order to reproduce the experimental results by numerical simulation via an extension of the Harmonic Balance Method.
This paper is organized as follows: First, the experimental setup of the clamped-clamped beam with non-ideal boundary conditions is presented with observations of correlation on temporal input data for four different types of correlation (in phase, opposite phase, quadrature phase and uncorrelated excitations). Then various experiments are proposed to highlight the presence of harmonics or a combination of harmonic components. Secondly, a nonlinear model of the beam with non-ideal boundary conditions and additional static pretension is introduced. Some notions about the modeling of two homogeneous correlated or uncorrelated random excitations are discussed. Then, an extension of the classical Harmonic Balance Method (HBM) is proposed for the prediction of the nonlinear dynamic behavior of a mechanical system subjected to two correlated or uncorrelated broadband random excitations. Finally, the numerical simulations are compared with experiments for the appearance and/or attenuation of some symmetric or anti-symmetric modes of the beam in the vicinity of their primary resonance, as well as the appearance of harmonics or combinations of harmonics.

## 2. Experiments

### 2.1. Experimental setup: the CEA-beam benchmark structure

The experimental clamped-clamped beam under study, called the CEA-beam benchmark structure, has been previously described in detail in [19,23]. The following is a brief reminder of the essential elements of description of the mechanical system under consideration.

The complete mechanical structure is presented in Fig. 1(a). The system is composed of a steel beam of dimensions 470 × 20 × 5 mm³ plus two heavy steel blocks of dimensions 100 × 100 × 85 mm³ each. In order to avoid undesirable nonlinear phenomena at the connections between substructures, the beam and the two blocks have been manufactured from a single bulk piece of steel. Due to the smooth transition between the beam and each block, the beam needs to be modeled with non-ideal boundary conditions (see the previous study by Claeys et al. [19] for more details). The whole structure (i.e., the beam plus the two blocks) is screwed to a heavy steel block.

The mechanical system is instrumented with four three-dimensional accelerometers (A1, A2, A3 and A4), two one-dimensional accelerometers (A5 and A6) and two cell forces (F1 and F2), as depicted in Fig. 1(b). Accelerometers A2, A3 and A4 are present to control and to validate the correct embedding of the whole structure during experiments. Indeed, the accelerations measured on these sensors appeared to be negligible, which validates the embedding hypothesis. The beam is subjected to a bi-point random excitation via two electrodynamic shakers (referenced as S1 and S2 in Fig. 1(b)) and using the Multi-Input-Multi-Output (MIMO) control technology [31].

It should be noted that the position of each shaker is fixed for all of the experiments, and symmetrical about the center of the beam. This choice will make it possible to highlight the contributions of symmetric or anti-symmetric modes of the beam, depending on the given correlation between the two broadband random excitations. For a chosen frequency band, four types of correlation between the two excitations will be tested: in phase, opposite phase, quadrature phase and uncorrelated.

At this stage of the study, we recall that many efforts have been made to enhance vibration testing [28]. Experimental modal analysis and vibration testing, in general, are classic techniques for obtaining the dynamic characteristics of linear engineering structures, for instance, resonant frequencies, mode shapes and modal damping. Such basic prior analysis is well established in the industry, in particular during the ground certification of aeronautical and aerospace structures for which the use of multipoint excitations with specific approaches, such as phase resonance or phase separation methods is classically conducted [31–37]. Extension to nonlinear structures is also becoming more and more studied and developed. As reported in [29,30], experiments for dynamic testing of nonlinear vibrating structure and characterization of nonlinearities from experiments have also been conducted in recent years with a focus on identifying nonlinear normal modes of engineering structures [38–44]. Therefore, we recall that the objective of the experimental studies carried out in this study is not to propose a contribution in this sense but only to provide a sufficiently complete and varied experimental database to allow us to test and validate the relevance of the proposed model and the numerical approach developed.

### 2.2. Preamble and modal testing

The purpose of this section is to provide some classical preliminary experimental results that will be used to validate the numerical model proposed in Section 3. Indeed, since the boundary conditions of the beam are not ideal and since many accelerometers are present on the structure, a first preliminary recalibration of the linear model of the beam is essential. These preliminary results will also provide useful basic information for a better understanding of the complex experimental and numerical results for the nonlinear beam subjected to two broadband correlated or uncorrelated random excitations.

First of all, Table 1 gives the natural frequencies (linear modes) of the beam system in the frequency range of interest [20; 1000] Hz. It should be noted that a classical impact hammer testing (in the z-direction) has been used and each natural frequency is calculated using the PolyMAX optimization algorithm [45] for a set of ten repeatable Frequency Response Functions (FRF) and three different impact points (i.e., at the locations of Accelerometers A1, A5 and A6). Moreover, the shape of each mode has been given by the Simcenter Testlab software (i.e., symmetric or anti-symmetric shape of each
<figure><img src="image_5.png" /><figcaption>Fig. 1. A picture (a) and a diagram (b) of the experimental setup which includes six accelerometers (A1 to A6) and two cell forces (F1, F2) connected to the two electrodynamic shakers S1 and S2.</figcaption></figure>

<table><caption>**Table 1**
Identified modes, corresponding frequency and damping ratio from experimental modal analysis.</caption><thead><tr><th>Mode number</th><th>Designation</th><th>Frequency (Hz)</th><th>Damping ratio (%)</th></tr></thead><tbody><tr><td>1</td><td>symmetric bending mode</td><td>108.9</td><td>0.25</td></tr><tr><td>2</td><td>anti-symmetric bending mode</td><td>307.0</td><td>0.12</td></tr><tr><td>1T</td><td>symmetric transverse mode</td><td>417.8</td><td>0.33</td></tr><tr><td>3</td><td>symmetric bending mode</td><td>594.9</td><td>0.13</td></tr><tr><td>4</td><td>anti-symmetric bending mode</td><td>997.3</td><td>0.17</td></tr></tbody></table>

<table><caption>**Table 2**
Measured mass of each accelerometer present on the beam.</caption><thead><tr><th>Accelerometer</th><th>Measured mass (g)</th></tr></thead><tbody><tr><td>A1</td><td>6.0</td></tr><tr><td>A5</td><td>1.1</td></tr><tr><td>A6</td><td>1.1</td></tr></tbody></table>

normal or transverse mode). The first four bending modes of the linear beam are detected, as well as a vibration transverse mode identified at 417.8 Hz. The identification of this transverse mode reveals a transverse component in the hammer impact. This transverse behavior will not be modeled; only the bending modes of the beam in the z-direction will be modeled. However, this first experiment indicates that when the excitation is not perfectly axial (in the z-direction), transverse modes may appear.

The mass of each accelerometer is given in Table 2. It should be noted that the mass of the accelerometers A3 and A4 (≤ 10 g) is not given, since they are situated on each embedding block (≈ 7000 g), so they do not have any influence on the mechanical behavior of the system.

## 2.3. Experimental results for two broadband random excitations

### 2.3.1. Preamble and validation of the temporal input data

One of the key points of the experiments carried out in this study is to be able to ensure that the chosen correlation between the two input signals is respected and well mastered. The objective of this section is to briefly present the selected experimental protocol and its implementation.

Each shaker is driven with a Power Spectral Density (PSD) and the correlation between the two excitations is defined in the Cross Spectral Density (CSD) with the notions of coherence and phase, according to Eq. (1)

$$S_{12}(f) = \sqrt{\gamma_{12}^2(f)S_{11}(f)S_{22}(f)} \exp(j\phi_{12}(f)) \quad (1)$$

where $S_{11}(f)$ and $S_{22}(f)$ are the PSD of the temporal excitations $F_1(t)$ and $F_2(t)$ (which correspond, respectively, to the force measured at the connection with Shakers 1 and 2, respectively). Here, $f$ corresponds to the frequency; $\phi_{12}(f)$ is the phase between the two temporal excitations $F_1(t)$ and $F_2(t)$; and $\gamma_{12}(f)$ defines the coherence between the two forces. It should be noted that preliminary experiments demonstrated that the experimental value of $\gamma_{12}$ varies between 0.05 and 0.98. Although this preliminary result implies that it is not possible in practice to generate two perfectly uncorrelated or correlated multipoint excitations with the electronic devices of this experimental setup, this emphasizes that the following protocol will be sufficient for generating the following desired broadband random excitations:
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Experimental inputs: (a) two experimental random temporal forces F1 (blue) and F2 (red) generated with the MIMO control technology [31] and (b) zoom on the four excitation configurations (from top to bottom: correlated excitations in phase, correlated excitations in opposite phase, correlated excitations in quadrature phase and uncorrelated excitations).</figcaption></figure>

* correlated excitations in phase: $\phi_{12} = 0^{\circ}$ and $\gamma_{12} = 0.98$;

* correlated excitations in opposite phase: $\phi_{12} = 180^{\circ}$ and $\gamma_{12} = 0.98$;

* correlated excitations in quadrature phase: $\phi_{12} = 90^{\circ}$ and $\gamma_{12} = 0.98$;

* uncorrelated excitations: $\gamma_{12} = 0.05$;

For each experiment that will be conducted, the shape of the band-limited white noise PSD is rectangular, i.e., the level of excitation is constant along the entire bandwidth. The sampling frequency is 3200 Hz with a resolution frequency of 0.098 Hz. The total duration of one experiment is 307.2 s. The output signal is divided into 30 time intervals of 10.24 s each. A periodogram estimate computes an approximation of the PSD for each time interval. Then, the PSD of the output signal is calculated by averaging all estimates [46]. The input PSD is constant over time to ensure that the response is stationary.

In order to illustrate and to validate the proposed protocol for the generation of input signals, Fig. 2(a) shows experimental temporal random inputs for one time interval of 10.24 s for the correlated excitations in phase. Moreover, zooms on experimental temporal random inputs are given in Fig. 2(b) for the four different correlation types (in phase, opposite phase, quadrature phase and uncorrelated). These preliminary results validate the experimental protocol for generating a bi-point correlated random excitation with the MIMO control technology.

The experiments carried out comprise two main parts:

* First of all, the global vibrational behavior of the beam system with non-ideal boundary conditions is investigated. More specifically, the effects of the four configurations of the two chosen broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations) are discussed. All experiments are performed with the same level of excitation ($\approx$ 0.71 N RMS) along the bandwidth [20; 1000] Hz. These first experiments will allow us to check the robustness of the proposed nonlinear model to reproduce the vibrational behavior of the beam system subjected to different broadband random excitations.

* Secondly, more particular attention is devoted to the nonlinear signature of the system and the evolution of the nonlinear contributions and harmonic components according to the different configurations of the broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations). In order to achieve such an objective, all experiments are conducted for the same level of excitation ($\approx$ 0.5 N RMS) along the bandwidth [20; 500] Hz and the nonlinear signature of the system is investigated along the bandwidth [500; 1000] Hz. This experimental procedure allows to highlight more particularly the nonlinear contributions over the frequency range [500; 1000] Hz. For the present study, it will allow us to verify the relevance of the nonlinear model as well as the associated numerical approach to reproduce the nonlinear contributions and their variations in function of the four types of correlation between the two broadband random excitations.

It is recalled that one of the original contributions of the present study is to share the data sets to give the opportunity to researchers to conduct additional analysis. All the data sets are available on [1].

### 2.3.2. Effect of the two broadband correlated or uncorrelated random excitations

The experimental results for Accelerometers A1, A5 and A6 for the same level of excitation ($\approx$ 0.71 N RMS) along the bandwidth [20; 1000] Hz are plotted in Fig. 3 for the three configurations of correlated random excitations and the
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_9.png" /><figcaption>**Fig. 3.** Experimental output PSD for Accelerometers A1 (a), A5 (b) and A6 (c) with RMS level of 0.71 N for four different correlation types along the bandwidth [20; 1000] Hz: correlated excitations in phase (blue), correlated excitations in opposite phase (red), correlated excitations in quadrature phase (green), and uncorrelated excitations (orange)</figcaption></figure>
<table><caption>**Table 3**
Frequency values of experimental output PSD for the four configurations of broadband random excitations with RMS level of 0.71 N along the bandwidth [20; 1000] Hz.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Frequency (Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>f<sub>1</sub></th><td>112.1</td><td>111.6</td><td>111.3</td><td>111.4</td></tr><tr><th>f<sub>2</sub></th><td>300.8</td><td>300.5</td><td>300.1</td><td>300.2</td></tr><tr><th>f<sub>3</sub></th><td>581.8</td><td>583.3</td><td>580.8</td><td>582.2</td></tr><tr><th>f<sub>4</sub></th><td>987.2</td><td>986.3</td><td>986.1</td><td>986.1</td></tr></tbody></table>

<table><caption>**Table 4**
Modal damping ratio values of experimental output PSD for the four configurations of broadband random excitations with RMS level of 0.71 N along the bandwidth [20; 1000] Hz.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Damping ratio (%)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>ξ<sub>1</sub></th><td>0.34</td><td>0.3</td><td>0.31</td><td>0.29</td></tr><tr><th>ξ<sub>2</sub></th><td>0.07</td><td>0.1</td><td>0.07</td><td>0.09</td></tr><tr><th>ξ<sub>3</sub></th><td>0.12</td><td>0.22</td><td>0.22</td><td>0.3</td></tr><tr><th>ξ<sub>4</sub></th><td>0.07</td><td>0.03</td><td>0.03</td><td>0.07</td></tr></tbody></table>

<table><caption>**Table 5**
Peak values of experimental output PSD for Accelerometer A1 with RMS level of 0.71 N for the four configurations of broadband random excitations along the bandwidth [20; 1000] Hz.</caption><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>f<sub>1</sub></th><td>689.54</td><td>3.59</td><td>310.05</td><td>523.51</td></tr><tr><th>f<sub>2</sub></th><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><th>f<sub>3</sub></th><td>287.22</td><td>1.27</td><td>111.55</td><td>111.77</td></tr><tr><th>f<sub>4</sub></th><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

configuration of uncorrelated random excitations previously defined in Section 2.3.1. The associated frequency values and the modal damping ratios of the first four vibration bending modes of the beam are given in Tables 3 and 4. Even if these values remain almost unchanged depending on the type of excitations, it can be noted that they are not exactly the same than those measured during the preliminary modal testing (see Table 1). This is mainly due to the presence of the two shakers, small variations in boundary conditions during experiments and a potential nonlinearities contribution [19,47]. To be noted that the potential large variation of modal damping ratios to the kind of excitation is mainly induced by the nonlinear behavior of the beam. As previously mentioned in [48–52] a growth of damping with the vibration amplitude during nonlinear vibrations can be observed during experiments. In addition the identification procedure has been checked by carrying out repeatability tests, even if small variations were also observed on modal damping ratios.

Then Tables 5–7 summarize the peak amplitudes of the experimental output PSD for Accelerometers A1, A5 and A6, respectively. It should be noted that the four largest amplitude peaks retained in this preliminary analysis correspond to amplitude peaks associated with the first four bending vibration modes of the beam with non-ideal boundary conditions. An attenuation or amplification of the system modes appears very clearly depending on the selected excitation configuration. The results obtained are consistent with a relationship between the shape of each eigenmode (symmetric or anti-symmetric vibration modes) and the type of excitation chosen (coherence and phase between the two excitations). For example, PSD peaks for the symmetric modes (anti-symmetric modes, respectively) are amplified (attenuated, respectively) in the case of correlated excitations in phase. Conversely, PSD peaks for the symmetric modes (anti-symmetric modes, respectively) are attenuated (amplified, respectively) in the case of correlated excitations in opposite phase. Moreover, in the two last cases of correlated excitations in quadrature phase or uncorrelated excitations, all of the symmetric and anti-symmetric modes in the frequency range of interest appear very clearly with non-negligible amplitude levels. It is also observed that anti-symmetric modes have no amplitude at the center of the beam (see Fig. 3(a) for $f_2$ and $f_4$) which corresponds to a classic expected result. Moreover, results for Accelerometers A5 and A6 are very similar because of the symmetry of the beam. Finally, some additional peaks that do not correspond to fundamental frequencies in the frequency range of interest (i.e., $f_1, f_2, f_3$ and $f_4$) are also observed. The appearance of such harmonic components and combinations of harmonics, as well as their evolutions with broadband correlated or uncorrelated random excitations, will be discussed in the next section.

### 2.3.3. Analysis of the nonlinear components

As described in [23], the proposed clamped–clamped beam subjected to broadband random excitations may induce a hardening effect and the enlargement of the response peak in the vicinity of the primary resonance, as well as the presence of secondary peaks resulting from the harmonics and a combination of harmonic components generated by the
<table><caption>**Table 6**

Peak values of experimental output PSD for Accelerometer A5 with RMS level of 0.71 N for the four configurations of broadband random excitations along the bandwidth [20; 1000] Hz.</caption><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>344.67</td><td>1.76</td><td>153.52</td><td>260.13</td></tr><tr><td>f<sub>2</sub></td><td>176.25</td><td>24581.10</td><td>15603.50</td><td>10104.3</td></tr><tr><td>f<sub>3</sub></td><td>115.91</td><td>0.76</td><td>48.22</td><td>46.61</td></tr><tr><td>f<sub>4</sub></td><td>0.98</td><td>72.28</td><td>39.54</td><td>20.72</td></tr></tbody></table>

<table><caption>**Table 7**

Peak values of experimental output PSD for Accelerometer A6 with RMS level of 0.71 N for the four configurations of broadband random excitations along the bandwidth [20; 1000] Hz.</caption><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>326.03</td><td>1.74</td><td>147.56</td><td>248.31</td></tr><tr><td>f<sub>2</sub></td><td>176.66</td><td>24724.50</td><td>15691.70</td><td>10114.9</td></tr><tr><td>f<sub>3</sub></td><td>133.54</td><td>0.63</td><td>48.11</td><td>53.70</td></tr><tr><td>f<sub>4</sub></td><td>0.64</td><td>47.57</td><td>27.76</td><td>13.62</td></tr></tbody></table>

primary resonance. In order to analyze in more detail the nonlinear behavior of the system for the different excitation configurations previously chosen and to highlight the contribution of the harmonics and combination of harmonic components, the following experimental protocol is proposed in addition to the previous tests: experiments are performed for one level of excitation (≈ 0.5 N RMS) along the bandwidth [20; 500] Hz, where only the first symmetric mode $f_1$ and the first anti-symmetric mode $f_2$ of the clamped-clamped beam are present, whereas a small level of excitation of approximately 0.0035 N RMS along the bandwidth [500; 1000] Hz is applied. It should be noted that the frequencies of the symmetric mode $f_3$ and the second anti-symmetric mode $f_4$ belong to this second frequency interval [500; 1000] Hz. Using this experimental protocol, the appearance of additional amplitude peaks (other than frequencies $f_3$ and $f_4$) between [500; 1000] Hz corresponds to an amplification of harmonics or combination of harmonic components of the beam.

Both experimental PSD inputs and outputs for the four excitation configurations (i.e., correlated excitations in phase, opposite phase and quadrature phase and uncorrelated excitations) are plotted in Fig. 4. The blue and red curves for the input PSD correspond to the shakers S1 and S2, respectively. The output PSD are those from Accelerometers A1 (red), A5 (blue) and A6 (green). First of all, the frequency values of the first four vibration bending modes of the beam are given in Table 8 and the associated modal damping ratios of the first two vibration bending modes are provided in Table 9. Comparing these frequency values with those of the previous experiments (see Table 3), a small increase in the frequency value of the four primary resonances with an increase of the excitation level is observed. This fact reflects a hardening effect. Then Tables 10–12 give the peak amplitudes of the experimental output PSD for the four primary resonances, respectively. It can be observed that the effects of the different excitation configurations are in accordance with the previous results discussed in Section 2.3.2. Considering more specifically the experimental input and output for the two correlated excitations in phase depicted in Figs. 4(a) and 4(b), several amplitude peaks that are not present in the experimental input signal along the bandwidth [500; 1000] Hz are observed. These results reveal a nonlinear dynamic behavior of the mechanical beam. As indicated in Fig. 5(a), some of these peaks correspond to harmonics of the fundamental frequency $f_1$ (i.e., $5f_1$, $6f_1$, $7f_1$, $8f_1$ and $9f_1$) and a combination of harmonic components such as $2f_1 + f_2$. The contributions of these harmonics are more important at the center of the beam (around Accelerometer A1). This result is consistent with the fact that the maximum amplitude is situated at the center of the beam. Once again, the results obtained for Accelerometers A5 and A6 are very close, due to the geometrical symmetry of the beam system under study and the symmetrical position of the sensors in relation to the center of the beam.

Thus, Figs. 4(c) and 4(d) show the experimental input and output for the two correlated excitations in opposite phase. Here again it is possible to detect several peaks of harmonic components: the two most important contributions correspond to the second and third harmonic components (i.e., $2f_2$ and $3f_2$) of the second fundamental anti-symmetric mode (see Fig. 5(b)). This result is to be put in relation with the fact that the chosen excitation highlights mainly the participation of the second mode (see also Fig. 4(d)). Moreover, three other minor contributions that correspond to the 5th, 6th and 7th harmonics of the first fundamental symmetric mode are visible (see $5f_1$, $6f_1$ and $7f_1$ in Fig. 5(b) for Accelerometer A1).

Similarly the experimental input and output for the two correlated excitations in quadrature phase (or for the two uncorrelated excitations, respectively) are given in Figs. 4(e) and 4(f) (Figs. 4(g) and 4(h), respectively). It can be noted that experimental measurements for both inputs and outputs are very similar for these two configurations. As shown in Figs. 5(c) and 5(d), the number of nonlinear components along the bandwidth [500; 1000] Hz is greater than in the two previous cases (to be compared with Figs. 5(a) and 5(b)). This is explained by the fact that, for these two excitation configurations, all of the modes (symmetric and anti-symmetric) present along the bandwidth [20; 500] Hz are solicited without favoring the amplification of one mode over the other. Thus, the harmonics of the first and second modes (i.e., $5f_1$,
<table><caption>**Table 8**
Frequency values of experimental output PSD for the four configurations of broadband random excitations with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Frequency (Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>f<sub>1</sub></th><td>110.3</td><td>110.5</td><td>110.3</td><td>110.4</td></tr><tr><th>f<sub>2</sub></th><td>298</td><td>298.1</td><td>298</td><td>298</td></tr><tr><th>f<sub>3</sub></th><td>575.1</td><td>-</td><td>576.5</td><td>577.8</td></tr><tr><th>f<sub>4</sub></th><td>970.1</td><td>970.8</td><td>970.1</td><td>970.7</td></tr></tbody></table>

<table><caption>**Table 9**
Modal damping ratio values of experimental output PSD for the four configurations of broadband random excitations with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Damping ratio (%)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>ξ<sub>1</sub></th><td>0.3</td><td>0.31</td><td>0.17</td><td>0.16</td></tr><tr><th>ξ<sub>2</sub></th><td>0.08</td><td>0.11</td><td>0.08</td><td>0.08</td></tr></tbody></table>

<table><caption>**Table 10**
Peak values of the experimental output PSD for Accelerometer A1 with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz for the four configurations of broadband random excitations.</caption><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><th>f<sub>1</sub></th><td>1.10<sup>3</sup></td><td>7.10<sup>0</sup></td><td>1.10<sup>3</sup></td><td>1.10<sup>3</sup></td></tr><tr><th>f<sub>2</sub></th><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><th>f<sub>3</sub></th><td>6.10<sup>-3</sup></td><td>-</td><td>4.10<sup>-3</sup></td><td>3.10<sup>-3</sup></td></tr><tr><th>f<sub>4</sub></th><td>3.10<sup>-6</sup></td><td>1.10<sup>-6</sup></td><td>3.10<sup>-6</sup></td><td>3.10<sup>-6</sup></td></tr><tr><th>2f<sub>1</sub> + f<sub>2</sub></th><td>2.10<sup>-4</sup></td><td>-</td><td>2.10<sup>-3</sup></td><td>2.10<sup>-3</sup></td></tr><tr><th>5f<sub>1</sub></th><td>6.10<sup>-3</sup></td><td>3.10<sup>-3</sup></td><td>3.10<sup>-3</sup></td><td>3.10<sup>-3</sup></td></tr><tr><th>2f<sub>2</sub></th><td>-</td><td>2.10<sup>-1</sup></td><td>9.10<sup>-2</sup></td><td>4.10<sup>-2</sup></td></tr><tr><th>6f<sub>1</sub></th><td>6.10<sup>-4</sup></td><td>1.10<sup>-4</sup></td><td>3.10<sup>-4</sup></td><td>3.10<sup>-4</sup></td></tr><tr><th>f<sub>1</sub> + 2f<sub>1</sub></th><td>-</td><td>-</td><td>2.10<sup>-4</sup></td><td>1.10<sup>-4</sup></td></tr><tr><th>7f<sub>1</sub></th><td>1.10<sup>-5</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><th>8f<sub>1</sub></th><td>1.10<sup>-5</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><th>3f<sub>2</sub></th><td>-</td><td>5.10<sup>-5</sup></td><td>2.10<sup>-5</sup></td><td>1.10<sup>-5</sup></td></tr><tr><th>9f<sub>1</sub></th><td>4.10<sup>-5</sup></td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

$6f_1, 7f_1, 2f_2$ and $3f_2$), as well as the combination of harmonics (i.e., $2f_1 + f_2$ and $f_1 + 2f_2$), are visible. It can also be noted that the contributions $2f_1 + f_2$ and $f_1 + 2f_2$ are more important for the correlated excitations in quadrature phase and the uncorrelated excitations, contrary to the two other excitation configurations, due to the fact that this choice of excitations does not lead to the attenuation of the participation of the first or second modes. For the interested reader, Tables 10–12 summarize the presence of each peak as well as the associated amplitude value of the experimental output PSD for all of the excitation configurations. All of these results demonstrate that the excitation configuration may induce changes on the nonlinear dynamic behavior of the beam system with the more or less marked appearance of nonlinear components.

Finally, it can be observed that the control of the input PSD is not perfect for the four excitation configurations tested with the presence of unexpected resonance peaks along the bandwidth [500; 1000] Hz in Figs. 4(a), 4(c), 4(e) and 4(g). This necessarily may lead to the presence of these additional peaks in the experimental output. This phenomenon may also be due to the coupling between the control at the two excitation points and the nonlinear response of the beam. For example, it can also be noted that the two shakers are positioned near the maxima of the second mode shape of the beam (i.e., the first anti-symmetric mode). Thus, the coupling between the experimental control and the experimental output is stronger for the configuration of correlated excitations in opposite phase (for which the contribution of the second mode is amplified), than for the configuration of correlated excitations in phase (for which the contribution of the second mode is attenuated). Likewise, the two configurations related to the correlated excitations in quadrature phase and the uncorrelated excitations lead to a solicitation of all modes. This implies the presence of a greater number of unexpected resonance peaks along the bandwidth [500; 1000] Hz than in the two previous cases (i.e., correlated excitations in opposite phase and correlated excitations in phase), but with slightly smaller amplitudes.
<figure><img src="image_10.png" /><figcaption>**Fig. 4.** Experimental input and output PSD with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz; (a, b) correlated excitations in phase; (c, d) correlated excitations in opposite phase; (e, f) correlated excitations in quadrature phase; (g, h) uncorrelated excitations; (a, c, e, g) input PSD from Shaker S1 (blue) and Shaker S2 (red); (b, d, f, h) output from Accelerometers A1 (red), A5 (blue) and A6 (green).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 5.** Experimental output PSD along the bandwidth [500; 1000] Hz for the experiment with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz for Accelerometers A1 (red), A5 (blue) and A6 (green); blue vertical lines: harmonics of the 1st mode; red vertical lines: harmonics of the 2nd mode; black vertical lines: combinations of harmonics; (a) correlated excitations in phase; (b) correlated excitations in opposite phase; (c) correlated excitations in quadrature phase; (d) uncorrelated excitations.</figcaption></figure>

<caption>**Table 11**</caption>

Peak values of the experimental output PSD for Accelerometer A5 with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz for the four configurations of broadband random excitations.

<table><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>6.10<sup>2</sup></td><td>3.10<sup>0</sup></td><td>5.10<sup>2</sup></td><td>5.10<sup>2</sup></td></tr><tr><td>f<sub>2</sub></td><td>1.10<sup>2</sup></td><td>1.10<sup>4</sup></td><td>1.10<sup>4</sup></td><td>6.10<sup>3</sup></td></tr><tr><td>f<sub>3</sub></td><td>3.10<sup>-3</sup></td><td>-</td><td>1.10<sup>-3</sup></td><td>1.10<sup>-3</sup></td></tr><tr><td>f<sub>4</sub></td><td>1.10<sup>-5</sup></td><td>6.10<sup>-4</sup></td><td>6.10<sup>-4</sup></td><td>1.10<sup>-3</sup></td></tr><tr><td>2f<sub>1</sub> + f<sub>2</sub></td><td>4.10<sup>-5</sup></td><td>-</td><td>3.10<sup>-4</sup></td><td>3.10<sup>-4</sup></td></tr><tr><td>5f<sub>1</sub></td><td>1.10<sup>-3</sup></td><td>7.10<sup>-4</sup></td><td>7.10<sup>-4</sup></td><td>7.10<sup>-4</sup></td></tr><tr><td>2f<sub>2</sub></td><td>-</td><td>1.10<sup>-1</sup></td><td>6.10<sup>-2</sup></td><td>3.10<sup>-2</sup></td></tr><tr><td>6f<sub>1</sub></td><td>1.10<sup>-4</sup></td><td>1.10<sup>-4</sup></td><td>1.10<sup>-4</sup></td><td>1.10<sup>-4</sup></td></tr><tr><td>f<sub>1</sub> + 2f<sub>1</sub></td><td>-</td><td>-</td><td>4.10<sup>-4</sup></td><td>2.10<sup>-4</sup></td></tr><tr><td>7f<sub>1</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>8f<sub>1</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3f<sub>2</sub></td><td>-</td><td>8.10<sup>-4</sup></td><td>3.10<sup>-4</sup></td><td>2.10<sup>-4</sup></td></tr><tr><td>9f<sub>1</sub></td><td>1.10<sup>-5</sup></td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<caption>**Table 12**</caption>

Peak values of the experimental output PSD for Accelerometer A6 with RMS level of 0.5 N along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz for the four configurations of broadband random excitations.

<table><thead><tr><th rowspan="2">Frequency</th><th colspan="4">PSD ((m s<sup>-2</sup>)<sup>2</sup>/Hz)</th></tr><tr><th>In phase</th><th>Opposite phase</th><th>Quadrature phase</th><th>Uncorrelated</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>6.10<sup>2</sup></td><td>3.10<sup>0</sup></td><td>5.10<sup>2</sup></td><td>5.10<sup>2</sup></td></tr><tr><td>f<sub>2</sub></td><td>1.10<sup>2</sup></td><td>1.10<sup>4</sup></td><td>1.10<sup>4</sup></td><td>6.10<sup>3</sup></td></tr><tr><td>f<sub>3</sub></td><td>3.10<sup>-3</sup></td><td>-</td><td>1.10<sup>-3</sup></td><td>1.10<sup>-3</sup></td></tr><tr><td>f<sub>4</sub></td><td>1.10<sup>-5</sup></td><td>3.10<sup>-4</sup></td><td>6.10<sup>-4</sup></td><td>3.10<sup>-4</sup></td></tr><tr><td>2f<sub>1</sub> + f<sub>2</sub></td><td>5.10<sup>-5</sup></td><td>-</td><td>5.10<sup>-4</sup></td><td>5.10<sup>-4</sup></td></tr><tr><td>5f<sub>1</sub></td><td>1.10<sup>-3</sup></td><td>9.10<sup>-4</sup></td><td>9.10<sup>-4</sup></td><td>9.10<sup>-4</sup></td></tr><tr><td>2f<sub>2</sub></td><td>-</td><td>1.10<sup>-1</sup></td><td>6.10<sup>-2</sup></td><td>3.10<sup>-2</sup></td></tr><tr><td>6f<sub>1</sub></td><td>8.10<sup>-5</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><td>f<sub>1</sub> + 2f<sub>1</sub></td><td>-</td><td>-</td><td>4.10<sup>-4</sup></td><td>2.10<sup>-4</sup></td></tr><tr><td>7f<sub>1</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>8f<sub>1</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3f<sub>2</sub></td><td>-</td><td>1.10<sup>-3</sup></td><td>3.10<sup>-4</sup></td><td>2.10<sup>-4</sup></td></tr><tr><td>9f<sub>1</sub></td><td>6.10<sup>-6</sup></td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

# 3. Modeling

## 3.1. Modeling of a multipoint correlated random excitation

This section details the modeling of a bi-point correlated random excitation. The main idea is to approximate a real excitation by a random Gaussian excitation [53,54] using the central limit theorem. This can be done by simulating a series (of finite length *p*) of cosine functions with weighted amplitudes and evenly spaced frequencies. One of the major advantages will be to be able to directly use this modeling in the numerical method (i.e., the Harmonic Balance Method [23]) implemented in the rest of the study.

For the interested reader, Shinozuka et al. [53] demonstrate that such an approach excitation has the following properties:

* the average of the simulated excitation is zero;

* the autocorrelation function of the simulated excitation tends to the autocorrelation function of the real excitation. As a consequence, using the Wiener-Khintchine relation, the PSD of the simulated excitation tends to the PSD of the real excitation.

* the convergence of the autocorrelation function of the simulated excitation is as $\frac{1}{p^2}$ to the autocorrelation function (the PSD function, respectively) of the real excitation.

As previously presented in Section 2.3.1, the modeling of a multipoint correlated random excitation starts with the cross-spectral density matrix $S(\omega)$ (where $\omega$ is the pulsation, which is defined for two temporal excitations $f_1(t)$ and $f_2(t)$ by

$$S(\omega) = \begin{bmatrix} S_{11}(\omega) & S_{12}(\omega) \\ S_{21}(\omega) & S_{22}(\omega) \end{bmatrix} \qquad (2)$$

where $S_{11}(\omega)$ and $S_{22}(\omega)$ are, respectively, the PSD of the temporal excitations $f_1(t)$ and $f_2(t)$. $S_{21}(\omega)$ is the cross-spectral density function between temporal excitations $f_1(t)$ and $f_2(t)$. The correlation type (i.e., in phase, opposite phase, quadrature phase or uncorrelated excitations) between the two excitations is defined by the extra-diagonal terms of $S(\omega)$ according to Eq. (1). $S(\omega)$ is non-negative definitive [53] and therefore, using the unique Cholesky decomposition

$$S(\omega) = H(\omega)H^*(\omega)^T \qquad (3)$$

where the notation $*$ indicates the complex conjugate and $^T$ is the transpose. $H(\omega)$ is a lower triangular matrix given by

$$H(\omega) = \begin{bmatrix} H_{11}(\omega) & 0 \\ H_{21}(\omega) & H_{22}(\omega) \end{bmatrix} \qquad (4)$$

where the $H_{ij}(\omega)$ coefficients can be determined with the principal minors of the matrix $S(\omega)$ [54]. We have

$$H_{11}(\omega) = S_{11}(\omega)^{\frac{1}{2}} \qquad (5)$$

$$H_{21}(\omega) = \frac{S_{21}(\omega)}{S_{11}(\omega)^{\frac{1}{2}}} \qquad (6)$$

$$H_{22}(\omega) = (S_{22}(\omega) - |H_{21}(\omega)|^2)^{\frac{1}{2}} \qquad (7)$$
Then, the two correlated temporal excitations $f_1(t)$ and $f_2(t)$ can be expressed as

$$f_1(t) = \sqrt{2\Delta\omega} \sum_{n=1}^{p} |H_{11}(\omega_n)| \cos(\omega'_n t + \theta_{11}(\omega_n) + \phi_{1n}) \quad (8)$$

$$f_2(t) = \sqrt{2\Delta\omega} \sum_{n=1}^{p} |H_{21}(\omega_n)| \cos(\omega'_n t + \theta_{21}(\omega_n) + \phi_{1n}) + \sqrt{2\Delta\omega} \sum_{n=1}^{p} |H_{21}(\omega_n)| \cos(\omega'_n t + \theta_{22}(\omega_n) + \phi_{2n}) \quad (9)$$

The parameters present in the two previous equations are defined in such a way that:

* $p$ is the finite length of the series of cosine functions of the approximated excitation. It has to be considered as a power of 2 in order to use Fast-Fourier Transforms (as detailed below for the numerical simulation);

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

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

* $\omega'_n$ is defined as $\omega'_n = \omega_n + \delta\omega_n$, where $\delta\omega_n$ is a small random pulsation introduced to avoid the periodicity of the simulated excitation [53]. In this study, in order to use the Harmonic Balance Method (see Section 4), we must have a periodic excitation so $\delta\omega_n = 0$ and $\omega'_n = \omega_n$;

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

$$\theta_{jk}(\omega_n) = \tan^{-1} \left( \frac{\Im(H_{jk}(\omega_n))}{\Re(H_{jk}(\omega_n))} \right) \quad (10)$$

where $\Im$ and $\Re$ are the imaginary and real parts, respectively. Given that $S(\omega)$ is a Hermitian matrix, we have $\theta_{11}(\omega) = \theta_{22}(\omega) = 0$. The dephasing between temporal excitations $f_1(t)$ and $f_2(t)$ comes from the variable $\theta_{21}(\omega_n)$;

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

From a numerical point of view (computational time and storage), the calculation of such a series may be expensive for high values of $p$. In order to avoid this, it is possible to rewrite Eqs. (8) and (9) by using the Fast-Fourier Transform (FFT) [55] as

$$f_1(t_n) = \sqrt{2\Delta\omega} \Re (\text{IFFT}(\mathbf{H}_{11}[n] e^{j\phi_1[n]})) \quad (11)$$

$$f_2(t_n) = \sqrt{2\Delta\omega} \Re (\text{IFFT}(\mathbf{H}_{21}[n] e^{j\phi_1[n]})) + \sqrt{2\Delta\omega} \Re (\text{IFFT}(\mathbf{H}_{22}[n] e^{j\phi_2[n]})) \quad (12)$$

where IFFT denotes the Inverse of the Fast Fourier Transform and $j$ is the pure imaginary number. $\mathbf{x}[n]$ indicates the $n$th term of the vector $\mathbf{x}$. $\mathbf{H}_{11}$, $\mathbf{H}_{21}$, $\mathbf{H}_{22}$, $\phi_1$ and $\phi_2$ are, respectively, the vectors defined by $\mathbf{H}_{11}[n] = H_{11}(\omega_n)$, $\mathbf{H}_{21}[n] = H_{21}(\omega_n)$, $\mathbf{H}_{22}[n] = H_{22}(\omega_n)$, $\phi_1[n] = \phi_{1n}$ and $\phi_2[n] = \phi_{2n}$. $t_n$ is the discretized time $t_n = ndt$ where $dt$ is the temporal sampling step. It should be noted that Poirion et al. [56] add another independent random variable for the weighted amplitudes in order to respect Box and Muller transformation [57] and to ensure that the simulated force is Gaussian for any fixed $p$, whereas in the study by Shinozuka et al. [53], the simulated force is only asymptotically Gaussian: i.e., when $p \to \infty$. In the present study, it has been chosen not to add these independent random variables.

In conclusion, using the experimental cross-spectral density matrix, it is possible to simulate two correlated or uncorrelated random excitations.

## 3.2. Modeling of the beam system with non-ideal boundary conditions and additional static pretension

The modeling was proposed initially by Nayfeh [58]. An extension of this model is proposed in the following.

The dynamical differential equation of the beam system with non-ideal boundary conditions and two broadband random excitations is defined by

$$\rho A \frac{\partial^2 w}{\partial t^2} + \mu \frac{\partial w}{\partial t} + EI \frac{\partial^4 w}{\partial x^4} = F_1(x_{F_1}, t) + F_2(x_{F_2}, t) + T(t) \frac{\partial^2 w}{\partial x^2} \quad (13)$$

where $\rho$ and $E$ are the mass density and the Young modulus of the beam, respectively. $A$ and $I$ are the cross-sectional area and moment of inertia. $\mu$ corresponds to the linear viscous damping coefficient. $w$ defines the transverse displacement in the reference frame of the beam. $F_1(x_{F_1}, t)$ and $F_2(x_{F_2}, t)$ are the first and second temporal forces situated at the abscissa $x_{F_1}$ and $x_{F_2}$, respectively. It can be noted that, thanks to the low thickness-length ratio of the experimental beam, it is legitimate to neglect the inertial and curvature nonlinear terms pointed out by Nayfeh [58]. For the interested reader a more complete modeling of a beam system including the geometric and internal nonlinearities as well as complete methodologies for solving such nonlinear beam system can be found in [5-9].
$T(t)$ corresponds to the tensile force, which is assumed to be invariant along the length of the beam. It is given by

$$T(t) = T_0 + EA \left( \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2 \right) \quad (14)$$

where $T_0$ is a static pretension and $u$ corresponds to the displacement along the beam axis. For the interested reader, the tensile force $T(t)$ can be assumed to be constant along the beam length due to the fact that the dimensionless quantity $\frac{R_g}{L}$ is very small for the beam under study (and consequently the longitudinal inertia is neglected, see [2], pages 446-455 for a detailed explanation). $L$ corresponds to the beam length (i.e., the effective length $L = 2l$) and $R_g$ denotes the radius of gyration given by $R_g = \sqrt{\frac{l}{A}}$ where $l$ is the second moment of area and $A$ the total cross-sectional area. Considering the geometrical parameters of the beam it comes $\frac{R_g}{L} = 0.003$ which leads without ambiguity to validate the hypothesis of a constant tensile force along the beam length.

The non-ideal boundary conditions of the beam model are defined as

$$T(0) = k_{bound}u(0); T(2l) = -k_{bound}u(2l) \quad (15)$$

$$w(0) = 0; w(2l) = 0 \quad (16)$$

$$EI \frac{\partial^2 w(0)}{\partial x^2} = k_{rot} \frac{\partial w(0)}{\partial x}; EI \frac{\partial^2 w(2l)}{\partial x^2} = -k_{rot} \frac{\partial w(2l)}{\partial x} \quad (17)$$

where the springs $k_{bound}$ and $k_{rot}$ are unknown parameters to be updated, due to the non-ideal boundary conditions (i.e., the deformable block, and each end of the beam being tapered [19]). It should be mentioned that there are several approaches for modeling of the non-ideal boundary conditions: the interested reader can referred to the paper of Rezaei and Zamanian for additional comments on this subject [14].

It is pointed out that only half of the beam has been modeled in the previous studies by Claeys et al. [19] and by Roncen et al. [23] due to the fact that the structure and its excitation are symmetric. In the present study, the two random excitations are not symmetric so the modeling of the beam is carried out over its total length $2l$. Another extension compared to the previous models by Claeys et al. [19] and by Roncen et al. [23] is the inclusion of the mass of Accelerometers A1, A5 and A6. Indeed, experimental tests (not presented in the present study for the sake of brevity) have shown a weak but not negligible influence of these added masses. Thus, taking into account these additional masses will allow a better correlation between experimental and numerical results. It should be noted that the effective length $2l$ is no longer an unknown parameter to be updated, thanks to the introduction of the static pretension $T_0$, unlike the previous studies [19,23]. This novelty will be presented and discussed in the rest of the paper.

Considering the boundary conditions given by Eq. (15), as well as by integrating $T(t)$ between $x = 0$ and $x = 2l$, the tensile force of the beam can be defined as a function that depends only on the transverse displacement $w(x, t)$. Indeed, it comes that

$$\int_{0}^{2l} T(t)dx = T(t)2l = T_02l + EA \left( \int_{0}^{2l} \frac{\partial u}{\partial x} dx + \frac{1}{2} \int_{0}^{2l} \left( \frac{\partial w}{\partial x} \right)^2 dx \right) \quad (18)$$

Then applying the boundary conditions defined by Eq. (15) on the longitudinal displacement $u(x, t)$, we have

$$\int_{0}^{2l} \frac{\partial u}{\partial x} dx = u(2l) - u(0) = -2 \frac{T(t)}{k_{bound}}$$

As a result the expression of the tensile force can be given by

$$T(t) = \left( 1 + \frac{EA}{lk_{bound}} \right)^{-1} \left( T_0 + \frac{EA}{4l} \int_{0}^{2l} \left( \frac{\partial w}{\partial x} \right)^2 dx \right) \quad (19)$$

and the dynamical nonlinear problem associated with the boundary conditions and the static pretension can be defined by

$$\frac{\partial^2 w}{\partial t^2} + \frac{\mu}{\rho A} \frac{\partial w}{\partial t} + \frac{EI}{\rho A} \frac{\partial^4 w}{\partial x^4} = \frac{1}{\rho A} F_1(x_{F1}, t) + \frac{1}{\rho A} F_2(x_{F2}, t) + \left(1 + \frac{EA}{lk_{bound}}\right)^{-1} \frac{T_0}{\rho A} \frac{\partial^2 w}{\partial x^2} \\ + \left(1 + \frac{EA}{lk_{bound}}\right)^{-1} \frac{E}{4l\rho} \int_{0}^{2l} \left(\frac{\partial w}{\partial x}\right)^2 dx \frac{\partial^2 w}{\partial x^2} \quad (20)$$

The next step is to obtain a set of discrete differential equations to be able to apply the numerical simulation process based on the Harmonic Balance Method used in this work. For this purpose, Eq. (20) can be projected on the modal basis of its associated homogeneous equation. After calculations, the projection yields the following discrete equation for the $i$th mode

$$\ddot{w}_i + \frac{\mu}{\rho A} \dot{w}_i + \omega_i^2 w_i = \Gamma_{i,1} F_1(x_{F1}, t) + \Gamma_{i,2} F_2(x_{F2}, t) + \sum_{j=1}^{N_m} \sum_{k=1}^{N_m} \sum_{m=1}^{N_m} \Gamma_{ijkm} w_j w_k w_m \quad (21)$$
with

$$ \Gamma_{i,1} = \frac{1}{\rho A} \int_{0}^{2l} \delta(x, x_{F1}) Y_{i}(x) dx = \frac{1}{\rho A} Y_{i}(x_{F1}) \qquad (22) $$

$$ \Gamma_{i,2} = \frac{1}{\rho A} \int_{0}^{2l} \delta(x, x_{F2}) Y_{i}(x) dx = \frac{1}{\rho A} Y_{i}(x_{F2}) \qquad (23) $$

$$ \Gamma_{ijkm} = -\frac{E}{4l\rho} \left(1 + \frac{EA}{lk_{bound}}\right)^{-1} a_j a_k a_m \int_0^{2l} \frac{dY_j(x)}{dx} \frac{dY_k(x)}{dx} dx \int_0^{2l} \frac{dY_i(x)}{dx} \frac{dY_m(x)}{dx} dx \qquad (24) $$

$$ a_i = \left[ \int_{0}^{2l} Y_i^2(x) dx \right]^{-1} \qquad (25) $$

where $\delta(x, a)$ is a Dirac delta function defined at abscissa $a$. To be noted that vibration modes of the beam are coupled in a nonlinear way in Eq. (21). Thus the transverse displacement $w(x, t)$ is defined by the product of two one-dimensional functions $w(x, t) = \sum_{i=1}^{N_m} a_i w_i(t) Y_i(x)$, where $N_m$ is the number of modes $Y_i(x)$ retained in the modal projection. As a reminder, $w_i(t)$ are the unknown functions to be determined (see Section 4 for the calculation of $w_i(t)$ using the Harmonic Balance Method). The modal solution $Y_i(x)$ can be expressed as

$$ Y_i(x) = \alpha_i \sin(\tau_{1,i}x) + \beta_i \cos(\tau_{1,i}x) + \gamma_i \sinh(\tau_{2,i}x) + \delta_i \cosh(\tau_{2,i}x) \qquad (26) $$

with

$$ \tau_{1,i} = \sqrt{\frac{-\alpha + \sqrt{\Delta_i}}{2 \frac{EI}{\rho A}}} \qquad (27) $$

$$ \tau_{2,i} = \sqrt{\frac{\alpha + \sqrt{\Delta_i}}{2 \frac{EI}{\rho A}}} \qquad (28) $$

$$ \Delta_i = \alpha^2 + 4 \frac{EI}{\rho A} \omega_i^2 \geq 0 \qquad (29) $$

$$ \omega_i = \sqrt{\frac{EI}{\rho A} \left(\frac{\lambda_i}{2l}\right)^4 - \alpha \left(\frac{\lambda_i}{2l}\right)^2} \qquad (30) $$

$$ \alpha = \left(1 + \frac{EA}{lk_{bound}}\right)^{-1} \frac{T_0}{\rho A} \qquad (31) $$

where $\lambda_i$ is the modal parameter of the ith vibration mode. Considering the previous non-ideal boundary conditions, the four variables $\alpha_i$, $\beta_i$, $\gamma_i$ and $\delta_i$ are a solution of

$$ \begin{bmatrix} 0 & 1 & 0 & 1 \\ \tau_{1,i} & \eta \tau_{1,i}^2 & \tau_{2,i} & -\eta \tau_{2,i}^2 \\ \sin(\tau_{1,i} 2l) & \cos(\tau_{1,i} 2l) & \sinh(\tau_{2,i} 2l) & \cosh(\tau_{2,i} 2l) \\ f_1(\tau_{1,i}, \eta) & f_2(\tau_{1,i}, \eta) & f_3(\tau_{2,i}, \eta) & f_4(\tau_{2,i}, \eta) \end{bmatrix} \begin{bmatrix} \alpha_i \\ \beta_i \\ \gamma_i \\ \delta_i \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \qquad (32) $$

with

$$ \eta = \frac{EI}{k_{rot}} \qquad (33) $$

$$ f_1(\tau_{1,i}, \eta) = -\eta \tau_{1,i}^2 \sin(\tau_{1,i} 2l) + \tau_{1,i} \cos(\tau_{1,i} 2l) \qquad (34) $$

$$ f_2(\tau_{1,i}, \eta) = -\eta \tau_{1,i}^2 \cos(\tau_{1,i} 2l) - \tau_{1,i} \sin(\tau_{1,i} 2l) \qquad (35) $$

$$ f_3(\tau_{2,i}, \eta) = \eta \tau_{2,i}^2 \sinh(\tau_{2,i} 2l) + \tau_{2,i} \cosh(\tau_{2,i} 2l) \qquad (36) $$

$$ f_4(\tau_{2,i}, \eta) = \eta \tau_{2,i}^2 \cosh(\tau_{2,i} 2l) + \tau_{2,i} \sinh(\tau_{2,i} 2l) \qquad (37) $$

Finally, each modal parameter $\lambda_i$ is determined by seeking the zeros of the determinant of the 4 × 4 matrix defined in Eq. (32), then each modal shape $Y_i$ can be obtained. The first four modal shapes for the problem under study (i.e., with the presence of a static pretension and non-ideal boundary conditions) are plotted in Fig. 6. Compared to previous work by Claeys et al. [19] and by Roncen et al. [23], all of the vibration modes are orthogonal for the problem defined in Eq. (20) with the proposed non-ideal boundary conditions defined in Eqs. (15), (16) and (17).
<figure><img src="image_12.png" /><figcaption>Fig. 6. The first four modal shapes along the length of the beam: 1st mode (blue) 2nd mode (red) 3rd mode (green) and 4th mode (orange).</figcaption></figure>

Finally, equations describing the dynamic behavior of the beam system are

$$ M\ddot{w}(t) + D\dot{w}(t) + Kw(t) = F(t) + F_{NL}(w(t)) \quad (38) $$

where **M** and **K** are the mass and stiffness matrices. The mass matrix **M** includes the mass contribution of the three accelerometers A1, A5 and A6, and is written as

$$ M = M^{\text{beam}} + M^{\text{mA1}} + M^{\text{mA5}} + M^{\text{mA6}} \quad (39) $$

with:

$$ M_{ij}^{\text{beam}} = \int_{0}^{2l} \rho A a_i a_j Y_i(x) Y_j(x) dx; \quad M_{ij}^{\text{mA1}} = m_{A1} a_i a_j Y_i(x_{A1}) Y_j(x_{A1}) \quad (40) $$

$$ M_{ij}^{\text{mA5}} = m_{A5} a_i a_j Y_i(x_{A5}) Y_j(x_{A5}); \quad M_{ij}^{\text{mA6}} = m_{A6} a_i a_j Y_i(x_{A6}) Y_j(x_{A6}) \quad (41) $$

where $m_{A1}$, $m_{A5}$ and $m_{A6}$ are the masses of Accelerometers A1, A5 and A6, respectively (situated at the abscissa $x_{A1}$, $x_{A5}$ and $x_{A6}$, respectively).

Similarly, the stiffness matrix **K** includes the static pretension and the non-ideal boundary conditions. We have

$$ K = K^{\text{beam}} + K^{k0} + K^{k2l} + K^{T0} \quad (42) $$

with

$$ K_{ij}^{\text{beam}} = \int_{0}^{2l} E l a_i a_j \frac{d^2 Y_i(x)}{dx^2} \frac{d^2 Y_j(x)}{dx^2} dx; \quad K_{ij}^{T0} = \int_{0}^{2l} T_0 2 l a_i a_j \frac{d Y_i(x)}{dx} \frac{d Y_j(x)}{dx} dx \quad (43) $$

$$ K_{ij}^{k0} = k_{\text{rot}} a_i a_j \frac{dY_i(0)}{dx} \frac{dY_j(0)}{dx}; \quad K_{ij}^{k2l} = k_{\text{rot}} a_i a_j \frac{dY_i(2l)}{dx} \frac{dY_j(2l)}{dx} \quad (44) $$

**D** corresponds to the damping matrix, which is defined by a modal damping ratio $\xi_i$ per vibration mode

$$ \forall i \in [1, N_m], \quad D_{ii} = 2\xi_i \Omega_i \quad (45) $$

where $\Omega_i$ corresponds to the natural pulsations of the generalized eigenvalue problem ($K - \omega^2 M$) $X = 0$. It is worth noting that $\Omega_i$ is slightly different from $\omega_i$, since additional masses are present in the model.

The generalized excitation vector $F(t)$ is defined as

$$ \forall i \in [1, N_m], \quad F_i(t) = \rho A \Gamma_{i,1} a_{f1}(t) + \rho A \Gamma_{i,2} a_{f2}(t) \quad (46) $$

where $\Gamma_{i,1}$, $\Gamma_{i,2}$ and $a_i$ are given by Eqs. (22), (23) and (25) for each ith vibration mode of the beam. Moreover, the vector $F_{NL}$ represents the nonlinear contributions given in Eq. (21)

$$ \forall i \in [1, N_m], \quad F_{NL,i}(t) = \rho A \sum_{j=1}^{N_m} \sum_{k=1}^{N_m} \sum_{m=1}^{N_m} \Gamma_{ijkm} w_j(t) w_k(t) w_m(t) \quad (47) $$

For the rest of the study, it was chosen to keep only the first five bending modes, i.e., $w(t) = [w_1(t) \ w_2(t) \ w_3(t) \ w_4(t) \ w_5(t)]^T$.
# 4. Simulation results

## 4.1. Preamble on the Harmonic Balance Method

In this section, the proposed numerical method based on the extension of the classical Harmonic Balance Method for random excitations will be presented. For more details, the interested reader is referred to [23,59].

The Harmonic Balance Method seeks the response of the nonlinear system defined in Eq. (38) as a truncated Fourier series (if this solution exists), such as

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

where $p$ corresponds to 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 the solution $\mathbf{w}$ to be determined. As discussed in [23], the Harmonic Balance Method can be adapted to random excitation by choosing the frequency resolution $\Delta f$ as the fundamental frequency of the excitation. In this case, the multipoint random excitation is foreseen as a multipoint deterministic excitation with one fundamental pulsation $\Omega = 2\pi \Delta f$ (see Section 3.1) and the vector force $\mathbf{F}(t)$ can be defined by a finite Fourier series of order $p$, such as

$$ \mathbf{F}(t) = \sum_{n=1}^{p} (\mathbf{S}_{n,\text{excit}} \sin(n\Omega t) + \mathbf{C}_{n,\text{excit}} \cos(n\Omega t)) \quad (49) $$

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

$$ \mathbf{F}_{\text{NL}}(t) = \mathbf{C}_0 + \sum_{n=1}^{p} (\mathbf{S}_n \sin(n\Omega t) + \mathbf{C}_n \cos(n\Omega t)) \quad (50) $$

where $(\mathbf{C}_0, (\mathbf{S}_n, \mathbf{C}_n)_{\forall n \in \llbracket 1, p \rrbracket})$ are the Fourier coefficients of the nonlinear force $\mathbf{F}_{\text{NL}}$.

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})$, the decompositions (48), (49) and (50) are re-injected into Eq. (38). This leads to a set of $N_m \times (2p + 1)$ nonlinear equations given by

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

$$ \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 (52) $$

The coefficients $(\mathbf{C}_0, (\mathbf{S}_n, \mathbf{C}_n)_{\forall n \in \llbracket 1, p \rrbracket})$ depend on the coefficients $(\mathbf{B}_0, (\mathbf{A}_n, \mathbf{B}_n)_{\forall n \in \llbracket 1, p \rrbracket})$. An extension of the classical Alternate Frequency-Time domain method (AFT-method [60]) is used to calculate these Fourier coefficients $(\mathbf{C}_0, (\mathbf{S}_n, \mathbf{C}_n)_{\forall n \in \llbracket 1, p \rrbracket})$ (see [23,59] for more details).

Finally, the nonlinear equations (51) and (52) are solved by minimizing the following relation

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

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 (54) $$

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

$$ \mathbf{B}_{\text{NL}} = [\mathbf{C}_0 \ \mathbf{S}_1 \ \mathbf{C}_1 \ \dots \ \mathbf{S}_p \ \mathbf{C}_p]^T \quad (56) $$

Minimizing Eq. (53) imposes a convergence criterion given by

$$ \frac{\|\mathbf{H}(\mathbf{X}, \Omega)\|_2}{\|\mathbf{B}\|_2} < \epsilon_{\text{HBM}} \quad (57) $$

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

## 4.2. Numerical results and comparison with the experiments

The objective of this section is to compare the simulation results with the experiments presented in Section 2.3. Specific attention will be paid to the understanding and analysis of the nonlinear dynamic behavior of the beam system. In order to achieve such an objective, results are decomposed into two main parts:
<table><caption>**Table 13**
Values of the physical parameters of the beam.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Mass density ρ</td><td>7850 kg m<sup>-3</sup></td></tr><tr><td>Cross-sectional area A</td><td>10<sup>-4</sup> m<sup>2</sup></td></tr><tr><td>Cross-sectional moment of inertia I</td><td>2.08 10<sup>-10</sup> m<sup>4</sup></td></tr><tr><td>Young modulus E</td><td>205 10<sup>9</sup> Pa</td></tr><tr><td>Half-length of the beam l</td><td>0.235 m</td></tr></tbody></table>

<table><caption>**Table 14**
Values of the numerical parameters for Case 1 and Case 2.</caption><thead><tr><td>Parameter</td><td>Case 1</td><td>Case 2</td></tr></thead><tbody><tr><td>Excitation bandwidth</td><td>[20; 1000] Hz</td><td>[20; 500] Hz</td></tr><tr><td>k<sub>rot</sub></td><td>2320 N m</td><td>1834 N m</td></tr><tr><td>k<sub>bound</sub></td><td>8.11 10<sup>7</sup> N m<sup>-1</sup></td><td>8.11 10<sup>7</sup> N m<sup>-1</sup></td></tr><tr><td>T<sub>0</sub></td><td>2809 N</td><td>3231 N</td></tr><tr><td>m<sub>A1</sub></td><td>0.006 kg</td><td>0.006 kg</td></tr><tr><td>m<sub>A5</sub></td><td>0.00171 kg</td><td>0.00171 kg</td></tr><tr><td>m<sub>A6</sub></td><td>0.00171 kg</td><td>0.00171 kg</td></tr></tbody></table>

* Case 1: comparison between numerical results and experiments are investigated for the dynamic behavior of the beam system with a constant PSD level along the bandwidth [20; 1000] Hz. The objective of this first analysis is to verify the ability of the simulations to reproduce the effects related to the four configurations of the two broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations). Numerical simulations are performed with an ideal excitation, as previously described in Section 3.1 (i.e., with a constant PSD level along the bandwidth [20; 1000] Hz).

* Case 2: more particular attention is devoted to the comparison of the nonlinear signature and the evolution of the harmonic components according to the four configurations of broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations). In order to achieve such an objective, two series of numerical simulations are performed: the first ones consider an ideal excitation with a high constant PSD level along the bandwidth [20; 500] Hz and a low constant PSD level along the bandwidth [500; 1000] Hz. The second ones are made by interpolating the measured experimental excitation for each configuration (as given and described in Figs. 4). This choice will allow us to assess the origin of the observed nonlinear effects on the nonlinear dynamic response of the beam system under study.

Tables 13 and 14 give the values of the physical parameters of the beam system under study and the physical values of the updated parameters, respectively. As previously described in [19,23], some variations of the fundamental frequencies of the beam (identified at low excitation) are observed from one experiment to another. Claeys et al. [19] assume that these frequency variations come from the static constraints in the beam induced by the two embedding blocks screwed to the heavy steel block and they proposed to update the eigenfrequency through an effective half-length of the beam. In the present study, readjustment of eigenfrequencies has been performed by updating the physical parameters given in Table 14. The modal damping ratio of the first five bending modes used for the numerical data are: $\xi_1 = 0.40\%$, $\xi_2 = 0.085\%$, $\xi_3 = 0.45\%$, $\xi_4 = 0.15\%$ and $\xi_5 = 1.5\%$. To be noted that the fifth bending mode is outside the frequency range of interest (with a frequency observed experimentally around 1464 Hz).

### 4.2.1. Case 1 — Reproduction of the effects related to the different configurations of the excitations by simulation

As previously explained, comparison between numerical results and experiments are first performed with a constant PSD level along the bandwidth [20; 1000] Hz.

Figs. 7(a-c), 7(b-d), 8(a-c) and 8(b-d) give the numerical results (blue curves) compared to experimental measurements (red curves) for the four configurations of the two chosen broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase, and uncorrelated excitations, respectively). Results are provided for Accelerometers A1 and A5. It should be noted that results for Accelerometer A6 are similar to those for Accelerometer A5 and therefore they are not provided for the sake of brevity.

It is observed that the numerical simulations are in good agreement with the experimental data. It is clearly shown that the effects related to the different configurations of the two broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase and uncorrelated excitations) are well reproduced by the numerical results. Depending on the selected excitation configuration, attenuation and amplification of each mode in the frequency range of interest [20; 1000] Hz are well reproduced. As previously discussed in Section 2.3.2, all of these results can be easily interpreted by considering the relationship between the shape of each eigenmode (symmetric or anti-symmetric vibration modes) and the type of excitation chosen. These first comparisons between experiments and numerical results validate the approach proposed for the modeling of a multipoint correlated random excitation as well as the modeling of the nonlinear beam with non-ideal boundary conditions and additional static pretension.
<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Experimental (red) and numerical (blue) output PSD for Accelerometers (a–b) A1 and (c–d) A5 with RMS level of 0.71 N for correlated excitations (a–c) in phase along the bandwidth [20; 1000] Hz and (b–d) in opposite phase along the bandwidth [20; 1000] Hz.</figcaption></figure>

We can also observe that the numerical simulation does not reproduce one experimental frequency peak of low amplitude (in the vicinity of 700 Hz). The appearance of this frequency peak is only due to an out-of-plane parasite excitation, which causes the excitation of a transverse mode of the beam. It results from experimental biases inherent in the experimental setup, such as the imperfect positioning of the two shakers.

### 4.2.2. Case 2 – Comparison of the nonlinear signature and harmonic components

In this section, the capability of the numerical simulation to reproduce the experimental nonlinear signature and the evolution of the nonlinear contributions and harmonic components according to the different configurations of broadband random excitations is investigated.

As previously explained in Section 2.3.3, an original experimental protocol has been proposed to analyze the nonlinear behavior of the beam system. The protocol is based on a high level of random excitation along the bandwidth [20; 500] Hz (where only the first symmetric mode $f_1$ and the first anti-symmetric mode $f_2$ are present) combined with a small level of excitation along the bandwidth [500; 1000] Hz. Therefore, the vibratory signature along the bandwidth [500; 1000] Hz brings out a complex behavior, which results from the appearance and the combination of harmonics.

Figs. 9, 10, 11 and 12 illustrate the numerical results (blue and green curves) compared with the experimental measurements (red curves) for the four configurations of the two chosen broadband random excitations. Zooms on the frequency band [500; 1000] Hz with identification of the harmonics and combination of harmonics are provided to better assess the comparison between experiments and numerical simulation. It should be noted that the numerical results defined by the blue lines refer to the use of an idealized excitation, whereas the numerical results given by the green lines refer to the use of an interpolation of the measured experimental excitation. To be noted that one of the assumptions made in this modeling is to consider a linear constant modal damping for all excitation cases. This choice seems relevant
<figure><img src="image_14.png" /><figcaption>**Fig. 8.** Experimental (red) and numerical (blue) output PSD for Accelerometers (a–b) A1 and (c–d) A5 with RMS level of 0.71 N for (a–c) correlated excitations in quadrature phase along the bandwidth [20; 1000] Hz and (b–d) for uncorrelated excitations along the bandwidth [20; 1000] Hz.</figcaption></figure>

and sufficient in the study proposed because of the good comparisons between the experiments and the numerical results for the four configurations of the two chosen broadband random excitations. For the interested reader, nonlinear damping identification and a potential observation of a damping growth with the vibration amplitude during nonlinear vibrations have been carried out in [48–52].

First of all, it is observed that the numerical simulations based on the interpolation of the measured experimental excitation are in good agreement with the experiments. Indeed, the appearance of additional amplitude peaks corresponding to the harmonics or combination of harmonic components are well predicted by the numerical simulation (see, for example, $nf_1$, $nf_2$ or $\pm nf_1 \pm mf_2$ where *n* and *m* are positive integers). Moreover, it can be noted that the evolution (i.e., the appearance, attenuation or amplification) of each harmonic or combination of harmonic components depending on the different configurations of the excitations are well reproduced. Also, it is interesting to note that the frequency peak of the symmetric mode $f_3$ and the second anti-symmetric mode $f_4$ belong to this second frequency interval [500; 1000] Hz. These small contributions for $f_3$ and $f_4$ are only due to the small level of excitation along the bandwidth [500; 1000] Hz. Once again, these comparisons between experiments and numerical results validate the proposed numerical strategy and the modeling of both the nonlinear beam system and the multipoint correlated random excitations.

Secondly, comparing the two series of numerical simulations (i.e., the first ones with an ideal excitation and the second one with an interpolation of the measured experimental excitation) can lead to interesting recommendations to be made on the modeling of the input signal. Indeed, it is clearly shown that considering an ideal excitation (with a high constant PSD level along the bandwidth [20; 1000] Hz, and a low constant PSD level along the bandwidth [500; 1000] Hz) does not allow the nonlinear signature of the beam system to be reproduced exactly, while an interpolation of the measured excitation allows a very good match with the experiments. This is especially true for the appearance of some specific frequency peaks, such as: $2f_1 + f_2$ and $5f_1$ for correlated excitations in phase; $2f_2$ for correlated excitations in
<figure><img src="image_15.png" /><figcaption>**Fig. 9.** Experimental (red) and numerical (blue, green) output PSD for Accelerometers (a–b) A1, (c–d) A5 and (e–f) A6 with RMS level of 0.5 N for correlated excitations in phase along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz — (b, d and f) Zoom on the frequency range [500; 1000] Hz.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 10.** Experimental (red) and numerical (blue, green) output PSD for Accelerometers (a–b) A1, (c–d) A5 and (e–f) A6 with RMS level of 0.5 N for correlated excitations in opposite phase along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz — (b, d and f) Zoom on the frequency range [500; 1000] Hz.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 11.** Experimental (red) and numerical (blue, green) output PSD for Accelerometers (a–b) A1, (c–d) A5 and (e–f) A6 with RMS level of 0.5 N for correlated excitations in quadrature phase along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz — (b, d and f) Zoom on the frequency range [500; 1000] Hz.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Experimental (red) and numerical (blue, green) output PSD for Accelerometers (a–b) A1, (c–d) A5 and (e–f) A6 with RMS level of 0.5 N for uncorrelated excitations along the bandwidth [20; 500] Hz and RMS level of 0.0035 N along the bandwidth [500; 1000] Hz – (b, d and f) Zoom on the frequency range [500; 1000] Hz.</figcaption></figure>
opposite phase: $2f_1 + f_2$, $5f_1$, $2f_2$ and $3f_2$ for correlated excitations in quadrature phase or uncorrelated excitations. These results demonstrate without any ambiguity that a fine modeling of the excitation signal is essential and that the intrinsic nonlinear model of the beam is not sufficient by itself to reproduce the complexity of the experimental results.

It is clearly shown that the effects related to the different configurations of the two broadband random excitations (i.e., correlated excitations in phase, opposite phase and quadrature phase and uncorrelated excitations) are well reproduced by the numerical results. Depending on the selected excitation configuration, the attenuation and amplification of each mode in the frequency range of interest [20; 1000] Hz are well reproduced. As previously discussed in Section 2.3.2, all of these results can be easily interpreted by considering the relationship between the shape of each eigenmode (symmetric or anti-symmetric vibration modes) and the type of excitation chosen. All these comparisons between experimental and numerical results validate the approach proposed for the modeling of a multipoint correlated random excitation, as well as the modeling of the nonlinear beam with non-ideal boundary conditions and additional static pretension.

Finally, some small differences are however noticeable between the simulated results and the experimental tests, in particular for Accelerometer A1. As previously discussed, the main reason is due to experimental biases from the experimental setup. In fact, the punctual external forces for the two shakers may not be perfectly applied in the vertical direction during experiments. This leads to additional modal contributions and interactions between bending modes and other out-of-plane modes, which are not taken into account in the model presented in this study.

## 5. Conclusion

Experiments and numerical simulations of a nonlinear clamped-clamped beam subjected to two broadband correlated, uncorrelated or partially correlated random excitations have been carried out in detail.

Based on experiments, the amplification and the attenuation of some symmetric or anti-symmetric vibration modes of the beam in the vicinity of their primary resonance are clearly shown (depending on the excitation configuration). Therefore, nonlinear phenomena have been identified as the presence of additional frequency peaks resulting from harmonics and a combination of harmonic components generated by the primary resonance of the vibration modes. One of the original contributions is to provide open data of all the experimental tests discussed in Section 2.3 [1].

Thus, a complete numerical strategy based on the formulation of a new modeling of the nonlinear beam with non-ideal boundary conditions and additional static pretension, as well as an extension of the Harmonic Balance Method for random excitations, has been proposed. The efficiency of the proposed methodology is illustrated by performing comparisons between experiments and numerical simulations for different types of random excitations.

Although the validity of the modeling choices and the efficiency of the proposed numerical strategy were unambiguously demonstrated by the high-quality comparisons between the experimental tests and the numerical results for several configurations, there are several directions that should be investigated for future research. Although it is not possible to give an exhaustive list of topics of interest, some studies that seem to have a certain priority and interest may be considered:

* It should be interesting to apply the proposed numerical strategy on a real industrial nonlinear systems subjected to correlated or uncorrelated broadband random excitations. In this context, one of the main challenges for engineers working on the developments of such nonlinear strategy for mechanical structures in a real environment should be the inclusions of modeling errors and uncertainties.

* A directly related point concerns the modeling of more realistic representation of real excitations such as random excitations applied to a surface and spatially correlated. One of the most common example of this kind of excitation source concerns industrial applications with turbulent boundary layer noises. For related nonlinear studies, it is necessary to reduce the number of excitation points while keeping an accurate prediction of the nonlinear response of the system.

* The proposed study focuses on the use of a numerical approach to predict the nonlinear vibrations of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations. Although the proposed strategy has been successfully validated due to the perfect agreement between experimental tests and numerical results, it could be interesting to compare these results with other alternative methodologies. To this end, the authors propose in [1] the data set for experiments discussed in this present study. This database gives the opportunity to researchers to validate analytical and numerical models for the prediction of the nonlinear dynamic behavior of the beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations.

## CRediT authorship contribution statement

S. Talik: Methodology, Software, Validation, Investigation, Writing – original draft, Writing – review & editing, Visualization. J.-J. Sinou: Conceptualization, Methodology, Validation, Writing – original draft, Writing – review & editing, Visualization, Supervision, Project administration. M. Claeys: Conceptualization, Methodology, Software, Validation, Investigation, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition. J.-P. Lambelin: Conceptualization, Methodology, Validation, Investigation, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition.
## Acknowledgments

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

## References

[1] Talik S, Claeys M, Lambelin J-P, Banvillet A, Sinou J-J. Dataset of measurements for the experimental CEA-beam benchmark structure subjected to two correlated or uncorrelated broadband random excitations. Data in Brief 2022.

[2] Nayfeh A, Mook D. Nonlinear oscillations. New York: John Willey And Sons; 1979.

[3] Krysko AV, Awrejcewicz J, Papkova IV, Szymanowska O, Krysko VA. Principal Component Analysis in the Nonlinear Dynamics of Beams: Purification of the Signal from Noise Induced by the Nonlinearity of Beam Vibrations. Adv Math Phys 2017;3038179:1-9.

[4] Zhang W, Wang F, Yao M. Global Bifurcations and Chaotic Dynamics in Nonlinear Nonplanar Oscillations of a Parametrically Excited Cantilever Beam. Nonlinear Dynam 2005;40:251-79.

[5] Awrejcewicz J, Saltykova OA, Zhigalov MV, Hagedorn P, Krysko AV. Analysis of Non-Linear Vibrations of Single-Layered Euler-Bernoulli Beams using Wavelets. Int J Aerosp Lightweight Struct 2011;1(2):203-19.

[6] Awrejcewicz J, Krysko AV, J. M. Saltykova OA, Zhigalov MV. Analysis of regular and chaotic dynamics of the Euler-Bernoulli beams using finite difference and finite element methods. Acta Mech Sinica 2011;27(1):36-43.

[7] Awrejcewicz J, Krysko AV, Soldatov V, Krysko VA. Analysis of the Nonlinear Dynamics of the Timoshenko Flexible Beams Using Wavelets. J Comput Nonlinear Dyn 2012;7:011005.

[8] Awrejcewicz J, Krysko AV, Kutepov IE, Zagniboroda NA, Dobriyan V, Krysko VA. Chaotic dynamics of flexible Euler-Bernoulli beams. Chaos 2013;23:043130.

[9] Awrejcewicz J, Krysko AV, Zagniboroda NA, Dobriyan VV, Krysko VA. On the general theory of chaotic dynamics of flexible curvilinear Euler-Bernoulli beams. Nonlinear Dynam 2015;79:11-29.

[10] Awrejcewicz J, Krysko AV, Papkova IV, Zakharov VM, Erofeev NP, Krylova EY, et al. Chaotic dynamics of flexible beams drive by external white noise. Mech Syst Signal Process 2016;79:225-53.

[11] Tabaddor M. Influence of nonlinear boundary conditions on the single-mode response of a cantilever beam. Int J Solids Struct 2000;37(36):4915-31.

[12] Pakdemirli M, Boyaci H. Non-linear vibrations of a simple-simple beam with a non-ideal support in between. J Sound Vib 2003;268(2):331-41.

[13] Eigoli AK, Ahmadian MT. Nonlinear vibration of beams under nonideal boundary conditions. Acta Mech 2011;218:259-67.

[14] Rezaei M, Zamanian M. A two-dimensional vibration analysis of piezoelectrically actuated microbeam with nonideal boundary conditions. Phys E Low Dimens Syst Nanostructures 2017;85:285-93.

[15] Huang Z, Wang Y, Zhu W, Huang Z. Deterministic and Random Response Evaluation of a Straight Beam with Nonlinear Boundary Conditions. J Vib Eng Technol 2020;82:847-57.

[16] Balasubramanian P, Franchini G, Ferrari G, Painter B, Karazis K, Amabili M. Nonlinear vibrations of beams with bilinear hysteresis ar supports: interpretation of experimental results. J Sound Vib 2021;499.

[17] Nayfeh A. Introduction to perturbation techniques. New York: Wiley; 2011.

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

[19] Claeys M, Sinou J-J, Lambelin J-P, Alcoverro B. Multi-harmonic measurements and numerical simulations of nonlinear vibrations of a beam with non-ideal boundary conditions. Commun Nonlinear Sci Numer Simul 2014;19:4196-212.

[20] Ye S-Q, Mao X-Y, Ding H, Ji J-C, Chen L-Q. Nonlinear vibrations of a slightly curved beam with nonlinear boundary conditions. Int J Mech Sci 2020;168:105294.

[21] Ribeiro P. Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods. Comput Struct 2004;82(17):1413-23.

[22] Ibrahim S, Patel B, Nath Y. Modified shooting approach to the non-linear periodic forced response of isotropic/composite curved beams. Int J Non Linear Mech 2009;44(10):1073-84.

[23] Roncen T, Lambelin J-P, Sinou J-J. Nonlinear vibrations of a beam with non-ideal boundary conditions and stochastic excitations - experiments, modeling and simulations. Commun Nonlinear Sci Numer Simul 2019;74:14-29.

[24] Givois A, Grolet A, Thomas O, Deü J-F. On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models. Nonlinear Dynam 2019;97:1747-81.

[25] Kandil A. Internal resonances among the first three modes of a hinged-hinged beam with cubic and quintic nonlinearities. Int J Non Linear Mech 2020;127:103592.

[26] Sayed M, Mousa A, Mustafa I. Stability and bifurcation analysis of a buckled beam via active control. Appl Math Model 2020;82:649-65.

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

[28] Daborn P, Ind P, Ewins D. Enhanced ground-based vibration testing for aerodynamic environments. Mech Syst Signal Process 2014;49:165-80.

[29] Kerschen G, Worden K, Vakakis AF, Golinval J-C. Past, present and future of nonlinear system identification in structural dynamics. Mech Syst Signal Process 2006;20(3):505-92.

[30] Noel J, Kerschen G. Nonlinear system identification in structural dynamics: 10 more years of progress. Mech Syst Signal Process 2017;83:2-35.

[31] Peeters B, Debille J. Multiple-input-multiple-output random vibration control: Theory and practice. In: Proceedings of the 2002 international conference on noise and vibration engineering. 2002. p. 507-16.

[32] Wright J, Cooper J, Desforges M. Normal-mode force appropriation-theory and application. Mech Syst Signal Process 1999;13:217-40.

[33] Peeters B, Climent H, de Diego R, de Alba J, Ahlquist J, Carreño J, et al. New Excitation Signals for Aircraft Ground Vibration Testing. In: Proceedings of the 26th international modal analysis conference. 2008.

[34] Peeters B, der Auweraer HV, Guillaume P, Leuridan J. The PolyMAX frequency-domain method: a new standard for modal parameter estimation? Shock Vib 2004;11(3-4):395-409.

[35] Göde D, Böswald M, Füllekrug U, Lubrina P. Ground vibration testing of large aircraft: state of the art and future perspectives. In: Proceedings of the 25th international modal analysis conference, Orlando, USA. 2007. p. 1-13.

[36] Grillenbeck A, Dillinger S. Reliability of Experimental Modal Data Determined on Large Spaceflight Structures. Adv Aerosp Appl 2011;351-61.

[37] Govers Y, Böswald M, Lubrina P, Giclais S, Stephan C, Botargues N. AIRBUS A350XWB Ground Vibration Testing: Efficient techniques for customer oriented on-site modal identification. In: Proceedings of the 26th international conference on noise and vibration engineering, Leuven, Belgium. 2014. p. 1-14.
[38] Gibert C. Fitting measured frequency response functions using non-linear modes. Mech Syst Signal Process 2003;17:211-8.

[39] Platten M, Wright J, Dimitriadis G, Cooper J. Identification of multi-degree of freedom non-linear systems using an extended modal space model. Mech Syst Signal Process 2009;23(1):8-29, Special Issue: Non-linear Structural Dynamics.

[40] Kerschen G, Peeters M, Golinval J, Vakakis A. Nonlinear normal modes, Part I: A useful framework for the structural dynamicist. Mech Syst Signal Process 2009;23(1):170-94, Special Issue: Non-linear Structural Dynamics.

[41] Peeters M, Kerschen G, Golinval J. Dynamic testing of nonlinear vibrating structures using nonlinear normal modes. J Sound Vib 2011;330(3):486-509.

[42] Peeters M, Kerschen G, Golinval J. Modal testing of nonlinear vibrating structures based on nonlinear normal modes: Experimental demonstration. Mech Syst Signal Process 2011;25(4):1227-47.

[43] Ehrhardt DA, Allen MS. Measurement of nonlinear normal modes using multi-harmonic stepped force appropriation and free decay. Mech Syst Signal Process 2016;76-77:612-33.

[44] Noël J-P, Renson L, Grappasonni C, Kerschen G. Identification of nonlinear normal modes of engineering structures under broadband forcing. Mech Syst Signal Process 2016;74:95-110, Special Issue in Honor of Professor Simon Braun.

[45] Peeters B. The PolyMAX frequency-domain method: a new standard for modal parameter estimation? Shock Vib 2004;11:395-409.

[46] Welch P. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans Audio Electroacoust 1967.

[47] Roncen T, Sinou J-J, Lambelin J-P. Non-linear vibrations of a beam with non-ideal boundary conditions and uncertainties – Modeling, numerical simulations and experiments. Mech Syst Signal Process 2018;110:165-79.

[48] Alijani F, Amabili M, Balasubramanian P, Carra S, Ferrari G, Garziera R. Damping for large-amplitude vibrations of plates and curved panels, Part 1: Modeling and experiments. Int J Non Linear Mech 2016;85:23-40.

[49] Balasubramanian P, Ferrari G, Amabili M. Identification of the viscoelastic response and nonlinear damping of a rubber plate in nonlinear vibration regime. Mech Syst Signal Process 2018;111:376-98.

[50] Amabili M. Nonlinear damping in nonlinear vibrations of rectangular plates: Derivation from viscoelasticity and experimental validation. J Mech Phys Solids 2018;118:275-92.

[51] Amabili M. Nonlinear damping in large-amplitude vibrations: modelling and experiments. Nonlinear Dynam 2018;93:5-18.

[52] Amabili M. Derivation of nonlinear damping from viscoelasticity in case of nonlinear vibrations. Nonlinear Dynam 2019;97:1785-97.

[53] Shinozuka M, Jan C-M. Digital simulation of random processes and its applications. J Sound Vib 1972;25(1):111-28.

[54] Shinozuka M. Simulation of multivariate and multidimensional random processes. J Acoust Soc Am 1971;49(1B):357-68.

[55] Vetterling W, Flannery B, Press W, Teukolsy S. Numerical Recipes in FORTRAN 77. 1996.

[56] Poirion F, Soize C. Numerical simulation of homogeneous and inhomogeneous Gaussian stochastic vector fields. English ed. La Recherche Aerospatiale; 1989.

[57] Box G, Muller M. A note on the generation of random normal deviates. Tech. rep., Statistical Techniques Research Group. Department of Mathematics. Princeton University; 1958.

[58] Nayfeh A. Nonlinear transverse vibrations of beams with properties that vary along the length. J Acoust Soc Am 1973;53:766-70.

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

[60] Cameron T, Griffin J. An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems. ASME J Appl Mech 1989;56:149-54.