# An efficient approach for predicting the nonlinear vibrations of a beam system subjected to multipoint correlated random excitation

S. Talik, M. Claeys, J.-J. Sinou, J.-P. Lambelin

► To cite this version:

S. Talik, M. Claeys, J.-J. Sinou, J.-P. Lambelin. An efficient approach for predicting the nonlinear vibrations of a beam system subjected to multipoint correlated random excitation. European Journal of Mechanics - A/Solids, 2022, 96, pp.104769. 10.1016/j.euromechsol.2022.104769. hal-03760567

HAL Id: hal-03760567

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

Submitted on 25 Aug 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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# An efficient approach for predicting the nonlinear vibrations of a beam system subjected to multipoint correlated random excitation

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

## ARTICLE INFO

### Keywords:

Vibration

Nonlinear beam

Multipoint correlated random excitation

## ABSTRACT

In linear vibration studies, the statistical description of multipoint random excitations is sufficient to calculate the vibration response of a structure. For nonlinear vibrations, it is necessary to model each excitation point separately, taking into account the correlation between each excitation point. The objective of this paper is to show how to reduce the number of excitation terms while remaining in a formalism compatible with nonlinear vibration studies.

The reduction of the number of stochastic excitation terms can be achieved by Galerkin methods (such as the Karhunen–Loève decomposition). This paper presents an original method which consists of projecting the excitation terms on the eigenmodes of the structure. These two methods are illustrated in the concrete case of a benchmark structure developed by the Commissariat à l'Energie Atomique (CEA), i.e., the mechanical beam system called the CEA-beam benchmark structure, previously studied in Talik et al. (2022), restrained to its first vibration mode and seen as a Duffing oscillator. A random excitation, composed of a consequent number of points of excitation distributed spatially along the structure (more exactly 101 points) and partially correlated, is used to illustrate the effectiveness of the proposed methodology. The proposed method makes it possible to reduce the number of random excitation signals to a single modal excitation term.

## 1. Introduction

Nonlinear vibration studies are becoming increasingly important due to the growing complexity of current mechanical systems. Researchers and engineers are now required to predict and guarantee the nonlinear vibration response of a structure subjected to a specific environment. The nonlinear vibration simulation of structures subjected to single point excitations is now well-known and assimilated by many researchers (Kerschen et al., 2006; Noël and Kerschen, 2017). Before studying complex engineering structures, mainly in the field of civil engineering and aerospace engineering, the potential of original methodologies for the prediction of the nonlinear responses of mechanical systems currently start with case studies of academic structures such as beams. For these structures, the classical results in structural dynamics are well-known and can be fairly easily reproduced, such as vibratory behavior under simple and controlled excitations. These numerous beam system studies allow us to investigate, on one hand, a different numerical approach and, on the other hand, increasingly complex and realistic excitations, while remaining within an acceptable numerical framework. In the case of advanced numerical

strategies, some authors investigate the prediction of the nonlinear vibrations of beams thanks to perturbation methods (Nayfeh and Mook, 1979; Nayfeh, 1993). Other methods, such as the iteration perturbation method (Eigoli and Ahmadian, 2011), shooting method (Ribeiro, 2004; Ibrahim et al., 2009; Roncen et al., 2019a) or Harmonic Balance Method (HBM) (Huang et al., 2011; Claeys et al., 2014; Ye et al., 2020) have been investigated in order to predict the nonlinear vibrations of beam systems. For the study of complex and realistic excitations, to the best of our knowledge, most of the published studies only cover multipoint excitations, which may be a harmonic excitation (Kandil, 2020; Sayed et al., 2020) or random excitation (Fang et al., 1995; Huang et al., 2020). Also, one of the major contributions of this work is to focus on the case of a beam system subjected to a multipoint correlated random excitation and more specifically to propose an efficient strategy for predicting the associated nonlinear vibration behavior.

It is to be noted that the beam system under study, namely the CEA-beam reference structure, has been previously considered for various analyses based on experiments, modeling and numerical simulations.

<footnote><sup>a</sup> CEA, DAM, CESTA, F-33114 Le Burp, 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>
In order to take into account the nonlinear vibrations of this clamped-clamped beam with non-ideal boundary conditions, Claeys et al. (2014) proposed a method of simulating nonlinear vibrations under harmonic or multi-harmonic excitation thanks to the HBM and the shooting method. Roncen et al. (2019a) then adapted the HBM to the one-point random excitation case. Both studies consider a beam model in which non-ideal boundary conditions have been introduced and are representative of the CEA-beam benchmark structure. In Talik et al. (2022), nonlinear vibrations of the CEA-beam subjected to two correlated or uncorrelated broadband random excitations were also investigated. In comparison with these previous works, one of the objectives of the study presented here is therefore to propose a generalization of the numerical strategy to simulate the nonlinear vibration behavior of the CEA-beam system subjected to multipoint correlated random excitation.

One of the major challenges for moving to an even more realistic representation of real excitations is to consider the case of random excitation. In such cases a more rational description of the vibratory responses of mechanical systems requires the use of stochastic nonlinear dynamic model for which specific developments are needed to predict the vibratory responses of mechanical system subjected to random excitation (Bellizzi and Bouc, 1999; Poirel and Price, 2007; Singh et al., 2016; Xu et al., 2016, 2017). For example, turbulent flows are one of the most common source of vibrations in aeronautics. This random excitation source is usually applied over a large surface and spatially correlated. This can lead to difficulties when solving this type of problem in an engineering context and the use of finite element models, both from the point of view of calculation time and the storage size required to solve such a problem. Indeed, applied to a finite element model, there are as many excitation terms as nodes concerned by the excitation. For many finite element models, especially for industrial structures, it can correspond to thousands of excitation terms. It is also worth noting that linear vibrations are well-adapted for multipoint excitations since the statistical description of each excitation term is sufficient to calculate the dynamical response of the structure, but these resolution methods are no longer available in the case of nonlinear structures. Also, the first contribution of this study is to propose the modeling of a multipoint correlated random excitation and a numerical methodology for predicting the nonlinear response of a mechanical system subjected to such complex excitation. One of the advantages of the proposed strategy is that it is based on the use of a numerical method, the HBM, which is easily implemented for finite element models and already used for industrial applications with nonlinearities (Claeys et al., 2016a,b; Roncen et al., 2019c, 2021). The second contribution of this study is to highlight and discuss the possibility of reducting the number of excitation terms remaining in a formalism compatible with nonlinear vibration studies. In the given literature, the Karhunen-Loève (KL) decomposition (Loeve, 1977) is one of the most popular methods and allows researchers to project a random field partially correlated in space onto a reduced number of specific vectors. In the following we propose an alternative method the projection basis is frequency-independent. This characteristic makes it possible to considerably reduce the number of calculations required. As a result the proposed approach provides a promising view of how to simulate the nonlinear vibrational behaviors of complex mechanical systems subjected to multipoint random excitations, whose correlation information depends on the frequency.

The following study is decomposed as follows. In the first part of this paper, the modeling of both the nonlinear beam system and a multipoint random excitation are described. In addition, the adaptation of HBM is presented in the case of to the beam model subjected to a multipoint random excitation with a high number of excitation terms. The proposed numerical results of the nonlinear vibrations of the beam subjected to 101 points' partially correlated random excitation will serve as a reference case to validate the proposed strategy in Section 3. The second part of this paper focuses on the proposed original method

<figure><img src="image_3.png" /><figcaption>Fig. 1. A diagram of the non-ideal clamped-clamped CEA-beam structure solicited by $N_e$ punctual correlated temporal forces.</figcaption></figure>

to reduce the number of excitation terms. This approach consists of projecting the excitation terms on the eigenmodes of the system. The results of the non-linear vibration responses are compared with previous reference results as well as those obtained by applying the KL decomposition (Loeve, 1977) coupled with the Shinozuka decomposition. This application demonstrates the ability of the proposed method to compute the nonlinear response of a system to a multipoint random excitation, with a reduced number of excitation terms, in order to keep a convenient numerical cost.

## 2. Case under study and reference vibrational responses

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

The modeling of the beam system under study is based on the previous analysis of a nonlinear beam with non-ideal boundary conditions introduced in Talik et al. (2022). The dynamical differential equation of this nonlinear beam system subjected to a multipoint correlated random excitation and a tensile force can be defined by (Nayfeh, 1973)

$$ \rho A \frac{\partial^2 w(x,t)}{\partial t^2} + \mu \frac{\partial w(x,t)}{\partial t} + EI \frac{\partial^4 w(x,t)}{\partial x^4} = F(x,t) + T(t) \frac{\partial^2 w(x,t)}{\partial x^2} \quad (1) $$

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(x,t)$ defines the transverse displacement in the reference frame of the beam. $F(x,t)$ is the multipoint correlated random excitation. Due to the low thickness-length ratio $\frac{e}{L}$ (with the thickness $e = 5$ mm and the length of the beam $L = 470$ mm) it can be assumed to neglect the inertial and curvature nonlinear terms (see Nayfeh (1973) for more details). $T(t)$ corresponds to the tensile force. As indicated in Nayfeh and Mook (1979) (see pages 446–455 for a detailed explanation), this tensile force can be assumed to be constant along the beam length due to the fact that the dimensionless quantity $\sqrt{\frac{I}{AL^2}}$ is very small (for the beam under study we have $\sqrt{\frac{I}{AL^2}} = 0.003$). This tensile force $T(t)$ 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 (2) $$

where $T_0$ is a static pretension and $u(x,t)$ corresponds to the displacement along the beam axis.

$N_e$ is defined as a positive integer which represents the number of excitation points present on the beam at abscissa $(x_k)_{\forall k \in [1, N_e]}$ such as

$$ x_k = (k - 1) \frac{L}{N_e - 1} \quad \forall k \in [1, N_e] \quad (3) $$

where $L$ is the length of the beam. In other words, the clamped-clamped beam is subjected to $N_e$ punctual correlated temporal forces evenly distributed along the beam as depicted in Fig. 1.
Therefore, the dynamical differential equation (1) becomes

$$ \rho A \frac{\partial^2 w(x,t)}{\partial t^2} + \mu \frac{\partial w(x,t)}{\partial t} + EI \frac{\partial^4 w(x,t)}{\partial x^4} = \sum_{k=1}^{N_e} f_k(t) \delta(x, x_k) + T(t) \frac{\partial^2 w(x,t)}{\partial x^2} \quad (4) $$

where $\delta(x, x_k)$ is a Dirac delta function which has the following property for any function $f$

$$ \int_{0}^{L} \delta(x, x_k) f(x) dx = \begin{cases} f(x_k) & \text{if } x_k \in [0, L] \\ 0 & \text{otherwise} \end{cases} \quad (5) $$

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

$$ T(t, x = 0) = k_{\text{bound}} u(t, x = 0) = T(t, x = L) = -k_{\text{bound}} u(t, x = L) \quad \forall t \quad (6) $$

$$ w(t, x = 0) = 0 \quad \text{and} \quad w(t, x = L) = 0 \quad \forall t \quad (7) $$

$$ EI \frac{\partial^2 w(t, x = 0)}{\partial x^2} = k_{\text{rot}} \frac{\partial w(t, x = 0)}{\partial x} \quad \text{and} \quad (8) $$

$$ EI \frac{\partial^2 w(t, x = L)}{\partial x^2} = -k_{\text{rot}} \frac{\partial w(t, x = L)}{\partial x} \quad \forall t $$

where the springs $k_{bound}$ and $k_{rot}$ are parameters related to the non-ideal boundary conditions. For more details and explanations on the process of achieving these parameters, the interested reader can refer to Claeys et al. (2014) and Talik et al. (2022). By considering these boundary conditions, as well as by integrating $T(t)$ between $x = 0$ and $x = L$, the tensile force of the beam can be defined as a function that depends only on the transverse displacement $w(x, t)$. It leads to the following expression of the tensile force (Talik et al., 2022).

$$ T(t) = \left(1 + \frac{2EA}{Lk_{\text{bound}}}\right)^{-1} \left(T_0 + \frac{EA}{2L} \int_0^L \left(\frac{\partial w}{\partial x}\right)^2 dx\right) \quad (9) $$

Considering the previous expressions (2) and (4), the dynamical nonlinear problem associated with the boundary conditions and the static pretension can be defined by

$$ \frac{\partial^2 w(x,t)}{\partial t^2} + \frac{\mu}{\rho A} \frac{\partial w(x,t)}{\partial t} + \frac{EI}{\rho A} \frac{\partial^4 w(x,t)}{\partial x^4} = \frac{1}{\rho A} \sum_{k=1}^{N_e} f_k(t) \delta(x, x_k) + \left(1 + \frac{2EA}{Lk_{\text{bound}}}\right)^{-1} \frac{T_0}{\rho A} \frac{\partial^2 w(x,t)}{\partial x^2} + \left(1 + \frac{2EA}{Lk_{\text{bound}}}\right)^{-1} \frac{E}{2L\rho} \int_0^L \left(\frac{\partial w(x,t)}{\partial x}\right)^2 dx \frac{\partial^2 w(x,t)}{\partial x^2} \quad (10) $$

In order to be able to apply the numerical simulation based on the HBM process that is implemented later, Eq. (10) can be projected onto the modal basis of its associated homogeneous equation (Talik et al., 2022). The following explanations explain the steps and calculations that lead to this projection of Eq. (10) into the final form given in Eq. (20). First of all Eq. (10) is projected on the modal basis of its associated homogeneous equation and 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 and $a_i$ are normalization constants defined below by Eq. (23). $w_i(t)$ are the unknown functions to be determined (see Section 2.3.2 for the calculation of $w_i(t)$ using the HBM). The $i$th normal mode of the beam $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) \quad (11) $$

The determination of the four variables $\alpha_i$, $\beta_i$, $\gamma_i$, $\delta_i$ can be performed by considering the non-ideal boundary conditions. This leads to the following system to solve

$$ \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} L) & \cos(\tau_{1,i} L) & \sinh(\tau_{2,i} L) & \cosh(\tau_{2,i} L) \\ 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} \quad (12) $$

<figure><img src="image_4.png" /><figcaption>Fig. 2. The first four modal shapes along the length of the beam: 1st mode (blue), 2nd mode (red), 3rd mode (green) and 4th mode (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

with

$$ f_1(\tau_{1,i}, \eta) = -\eta \tau_{1,i}^2 \sin(\tau_{1,i} L) + \tau_{1,i} \cos(\tau_{1,i} L) \quad (13) $$

$$ f_2(\tau_{1,i}, \eta) = -\eta \tau_{1,i}^2 \cos(\tau_{1,i} L) - \tau_{1,i} \sin(\tau_{1,i} L) \quad (14) $$

$$ f_3(\tau_{2,i}, \eta) = \eta \tau_{2,i}^2 \sinh(\tau_{2,i} L) + \tau_{2,i} \cosh(\tau_{2,i} L) \quad (15) $$

$$ f_4(\tau_{2,i}, \eta) = \eta \tau_{2,i}^2 \cosh(\tau_{2,i} L) + \tau_{2,i} \sinh(\tau_{2,i} L) \quad (16) $$

$$ \eta = \frac{EI}{k_{\text{rot}}} \quad (17) $$

In addition, the expressions for $\tau_{1,i}$ and $\tau_{2,i}$ are

$$ \tau_{1,i} = \sqrt{\frac{-\alpha + \sqrt{\Delta_i}}{2 \frac{EI}{\rho A}}} \quad \text{and} \quad \tau_{2,i} = \sqrt{\frac{\alpha + \sqrt{\Delta_i}}{2 \frac{EI}{\rho A}}} \quad (18) $$

with

$$ \begin{gathered} \Delta_i = \alpha^2 + 4 \frac{EI}{\rho A} \Omega_i^2 \ge 0 \quad ; \quad \Omega_i = \sqrt{\frac{EI}{\rho A} \left(\frac{\lambda_i}{L}\right)^4 - \alpha \left(\frac{\lambda_i}{L}\right)^2} \quad ; \\ \alpha = \left(1 + \frac{2EA}{Lk_{\text{bound}}}\right)^{-1} \frac{T_0}{\rho A} \end{gathered} \quad (19) $$

where $\lambda_i$ is the modal parameter of the $i$th vibration mode. For the interested reader, the first four modal shapes of the beam system under study are plotted in Fig. 2.

After calculations, the projection yields the following discrete equation for the $i$th mode

$$ \dot{w}_i(t) + 2\xi_i^2 \Omega_i \dot{w}_i(t) + \Omega_i^2 w_i(t) = \sum_{k=1}^{N_e} \Gamma_{i,k} f_k(t) + \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 (20) $$

with

$$ \Gamma_{i,k} = \frac{1}{\rho A} Y_i(x_k) \quad (21) $$

$$ \Gamma_{ijkm} = -\frac{E}{2L\rho} \left(1 + \frac{2EA}{Lk_{\text{bound}}}\right)^{-1} a_j a_k a_m \int_0^L \frac{dY_j(x)}{dx} \frac{dY_k(x)}{dx} dx \times \int_0^L \frac{dY_i(x)}{dx} \frac{dY_m(x)}{dx} dx \quad (22) $$
$$a_i = \left[ \int_0^L Y_i^2(x) dx \right]^{-1} \qquad (23)$$

where $\xi_i$ is the modal damping of the ith vibration mode. To be noted that Eq. (20) refers to a Galerkin projection made over the space spanned by the transverse linear modes of the beam. The two right terms of Eq. (20) correspond to the punctual correlated temporal forces and the tensile force with the non-ideal boundary conditions of the beam, respectively.

For the rest of the study, equations describing the system dynamic behavior of the nonlinear beam system, with $w(t) = [w_1(t) \dots w_{N_m}(t)]^T$, are given by

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

where, according to Eq. (20), we have

$$M = \begin{bmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & 1 \end{bmatrix} \qquad (25)$$

$$D = \begin{bmatrix} 2\xi_1 \Omega_1 & 0 & \dots & 0 \\ 0 & 2\xi_2 \Omega_2 & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ 0 & \dots & 0 & 2\xi_{N_m} \Omega_{N_m} \end{bmatrix} \qquad (26)$$

$$K = \begin{bmatrix} \Omega_1^2 & 0 & \dots & 0 \\ 0 & \Omega_2^2 & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ 0 & \dots & 0 & \Omega_{N_m}^2 \end{bmatrix} \qquad (27)$$

$$F(t) = \begin{bmatrix} \sum_{k=1}^{N_e} \Gamma_{1,k} f_k(t) \\ \sum_{k=1}^{N_e} \Gamma_{2,k} f_k(t) \\ \vdots \\ \sum_{k=1}^{N_e} \Gamma_{N_m,k} f_k(t) \end{bmatrix} \qquad (28)$$

$$F_{NL}(w(t)) = \begin{bmatrix} \sum_{j=1}^{N_m} \sum_{k=1}^{N_m} \sum_{m=1}^{N_m} \Gamma_{1jkm} w_j(t) w_k(t) w_m(t) \\ \sum_{j=1}^{N_m} \sum_{k=1}^{N_m} \sum_{m=1}^{N_m} \Gamma_{2jkm} w_j(t) w_k(t) w_m(t) \\ \vdots \\ \sum_{j=1}^{N_m} \sum_{k=1}^{N_m} \sum_{m=1}^{N_m} \Gamma_{N_m jkm} w_j(t) w_k(t) w_m(t) \end{bmatrix} \qquad (29)$$

Now that the dynamical nonlinear problem is set up, the expression of each temporal force $f_k(t)$ will be developed in the next section.

## 2.2. Modeling of a multipoint correlated random excitation

Shinozuka's modeling of a multipoint correlated random excitation (Shinozuka and Jan, 1972; Shinozuka, 1971) is now detailed. As in Section 2.1, we define $N_e$ as the positive integer which corresponds to the number of excitation points on the given domain. In this study, the excitation points are the points present on the beam at abscissa $(x_k)_{\forall k \in \{1, N_e\}}$, and the given domain is the whole length [0, L] of the beam. The main idea of this modeling is to approximate a real multipoint excitation by $N_e$ random Gaussian excitations (Shinozuka and Jan, 1972; Shinozuka, 1971; Shinozuka and Deodatis, 1991) using the central limit theorem. This can be done by simulating $N_e$ series (of finite length $p$ where $p$ is a positive integer) 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 HBM Roncen et al., 2019a) used for nonlinear simulations.

For the interested reader, Shinozuka and Jan (1972) demonstrate that such an approached multipoint excitation has the following properties

* the average of each simulated excitation is zero;

* the autocorrelation function of each simulated excitation tends towards the autocorrelation function of the real excitation applied to each excitation point. As a consequence, using the Wiener-Khintchine relation, the Power Spectral Density (PSD) of each simulated excitation tends towards the PSD of the real excitation applied to each excitation point.

* the convergence of the autocorrelation function of the simulated excitation is as $1/p^2$ to the autocorrelation function (the PSD function, respectively) of the real excitation applied to each excitation point.

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 $N_e$ temporal excitations by

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

where $S_{ii}(\omega)$ is the PSD of the temporal excitation $f_i(t)$ (which is the temporal force applied at abscissa $x_i$) and $S_{ij}(\omega)$ is the cross-spectral density function between temporal excitations $f_i(t)$ and $f_j(t)$. Thanks to the Wiener-Khintchine relation, $S(\omega)$ is positive definite (Shinozuka and Jan, 1972; Shinozuka and Deodatis, 1991) and, using the Cholesky decomposition, one may find

$$S(\omega) = H(\omega)H^*(\omega) \qquad (31)$$

where the notation $*$ denotes the conjugate transpose. $H(\omega)$ is a lower triangular matrix given by

$$H(\omega) = \begin{bmatrix} H_{11}(\omega) & 0 & \dots & 0 \\ H_{21}(\omega) & H_{22}(\omega) & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ H_{N_e 1}(\omega) & H_{N_e 2}(\omega) & \dots & H_{N_e N_e}(\omega) \end{bmatrix} \qquad (32)$$

where the $H_{ij}(\omega)$ coefficients can be determined with the principal minors of the matrix $S(\omega)$ (Shinozuka, 1971). Then, the temporal excitation $f_k(t)$ is defined by

$$f_k(t) = \sqrt{2\Delta\omega} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |H_{kl}(\omega_n)| \cos(\omega'_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) \qquad (33)$$

The parameters present in the previous Eq. (33) 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 (Vetterling et al., 1992);

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

* $\omega_n$ corresponds to the frequency discretization (Shinozuka and Jan, 1972). 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 (Vetterling et al., 1992);

* $\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 (Shinozuka and Jan, 1972). In this study, in order to use the HBM (see Section 2.3.4), we must have a periodic excitation so $\delta\omega_n = 0$ and $\omega'_n = \omega_n$;
• $\theta_{kl}(\omega_n)$ is defined as

$$ \theta_{kl}(\omega_n) = \tan^{-1} \left( \frac{\mathfrak{F}(H_{kl}(\omega_n))}{\mathfrak{R}(H_{kl}(\omega_n))} \right) \quad (34) $$

where $\mathfrak{F}$ and $\mathfrak{R}$ are the imaginary and real parts, respectively. Given that $S(\omega)$ is Hermitian, we have $\theta_{kk}(\omega) = 0$. The dephasing between both temporal excitations $f_k(t)$ and $f_l(t)$ comes from the variable $\theta_{kl}(\omega_n)$;

• $\forall k \in [1, N_e]$, $\phi_{kn}$ 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, since $p$ is chosen as a power of 2 and since the frequency discretization is correctly defined (Vetterling et al., 1992), it is possible to rewrite Eq. (33) by using the Fast-Fourier Transform (FFT) such as

$$ f_k(t_n) = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{l=1}^{k} H_{kl}(\omega_n) \exp(j\phi_{ln}) \right) \right) \quad (35) $$

where IFFT denotes the Inverse of the Fast Fourier Transform and $j$ is a complex number. $t_n$ is the discretized time $t_n = n \cdot dt$ where $dt$ is the temporal sampling step.

In conclusion, using the cross-spectral density matrix $S(\omega)$, it is possible to simulate a multipoint correlated random excitation in the time domain. It should be noted that the multipoint random temporal excitation is defined in the canonical basis of the excitation points defined in the domain $\Omega$.

It is to be noted that two other models of a multipoint random excitation, defined in another basis than the canonical basis, will be also given in Section 3. These two approaches (based on the KL technique and the projection of the excitation on the eigenmodes of the structure, respectively) will be introduced in the last part of this present work in order to allow a reduction of the number of excitation terms.

## 2.3. Application and reference vibrational responses

The objective of this part of the paper is to provide reference results on the non-linear vibrations of the beam system subjected to multipoint random excitation. To this purpose, an extension of the HBM adapted to the use of multi-point correlated random excitations is presented, as well as preliminary results on a linear case in order to verify the relevance of the implemented approach. Then the reference results of the nonlinear response of the beam system to a multipoint correlated random excitation are presented. These results will serve as a reference to validate the reduction of the number of excitation terms through the methods that will be proposed in Section 3.

### 2.3.1. Preamble

The modeling from Shinozuka's work (Shinozuka, 1971) is adapted for the clamped-clamped beam presented in Section 2.1 subjected to a multipoint correlated random excitation. The excitation is defined with a level of excitation $S(\omega)$ along a finite frequency band and a correlation function $C_{ij}(\omega)$ defined between two excitation points $x_i$ and $x_j$. The excitation is considered as identical for all excitation points

$$ S_{ii}(\omega) = S_{jj}(\omega) = S(\omega) \quad \forall (i, j) \in [1, N_e]^2 \quad (36) $$

In order to observe the same nonlinear phenomena as observed in Roncen et al. (2019a), the level of excitation $S(\omega)$ is defined along the finite frequency band [20; 500] Hz and chosen to be rectangular with a high level of excitation around the first mode of the beam along the frequency band [80; 130] Hz (2.60 $10^{-6}$ N$^2$/Hz), and with a low level of excitation along the other frequency bands [20; 80] Hz and [130; 500] Hz (1.25 $10^{-10}$ N$^2$/Hz).

<table><caption>**Table 1**

Values of the geometrical and material parameters of the beam.</caption><thead><tr><th>Parameter</th><th>Value</th></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>Length of the beam L</td><td>0.470 m</td></tr><tr><td>Rotational spring k<sub>rot</sub></td><td>1834 N m</td></tr><tr><td>Linear spring k<sub>bound</sub></td><td>8.11 10<sup>7</sup> N m<sup>-1</sup></td></tr><tr><td>Static pretension T<sub>0</sub></td><td>3231 N</td></tr></tbody></table>

The correlation function $C_{ij}(\omega)$ between two excitation points $x_i$ and $x_j$ is geometrical and considered to be constant for the whole frequency band such as

$$ C_{ij}(\omega) = C_{ij} = \exp\left(-\frac{|x_i - x_j|}{d}\right) \quad (37) $$

where both positions $x_i$ and $x_j$ are defined by Eq. (3). $d$ is a fixed correlation distance. This correlation function is inspired by the correlation functions used for the studies of turbulent boundary layer noises (see Efimtsov (1982) with $\omega \to 0$). The cross-spectral density matrix $S(\omega)$ from Eq. (30) is

$$ S(\omega) = S(\omega)C = S(\omega) \begin{bmatrix} C_{11} & C_{12} & \cdots & C_{1N_e} \\ C_{21} & C_{22} & \cdots & C_{2N_e} \\ \vdots & \vdots & \ddots & \vdots \\ C_{N_e 1} & C_{N_e 2} & \cdots & C_{N_e N_e} \end{bmatrix} \quad (38) $$

For the rest of this study, a number of 101 excitation points ($N_e = 101$) is chosen, and the correlation distance $d$ is defined such as $d/L = 0.1$. To emphasize the objective of validating a numerical modeling for multipoint correlated random excitations and in order to simplify the numerical study of the structure, only the first symmetric bending mode of the beam is considered $N_m = 1$ (for which the natural frequency is situated at $f_1 = 113.0$ Hz). The beam may therefore be seen as a Duffing oscillator (Kovacic and Brennan, 2011) subjected to $N_e$ correlated broadband random excitations.

The geometrical and material parameters used for the description of the beam are given in Table 1. It is to be noted that the complete geometry and technical drawing of the benchmark CEA-beam are available in Claeys et al. (2019).

### 2.3.2. Numerical method

The numerical method used to predict the nonlinear vibrations of a beam subjected to a multipoint random excitation is an extension of the HBM adapted to random excitations. For this purpose, it is necessary to adapt this nonlinear method for the use of multipoint correlated random excitations. This point is briefly presented below.

The HBM seeks the response of the nonlinear system defined in Eq. (24) 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 (39) $$

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 [1, p]})$ are the unknown Fourier coefficients of solution $\mathbf{w}(t)$ to be determined. As discussed in Talik et al. (2022) and Roncen et al. (2019a), the HBM can be adapted to random excitations by choosing 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$. To be noted that the frequency step $\Delta f$ is directly linked to the previously defined pulsation step $\Delta\omega$ (see Section 2.2). Then 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 (40) $$
Moreover, it is assumed that nonlinear contributions $F_{NL}(t)$ can be solved in the finite Fourier series of order $p$

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

where $(C_0, (S_n, C_n)_{\forall n \in [1, p]})$ are the Fourier coefficients of the nonlinear force $F_{NL}(t)$. One of the consequences of the choice of $\Delta f$ as the fundamental frequency is that the number of harmonics $p$ retained in the approximate nonlinear solution (39) can be consequent, which can constitute a major disadvantage for computation time and data storage issues. In order to determine the value of the Fourier coefficients $(B_0, (A_n, B_n)_{\forall n \in [1, p]})$, decompositions (39), (40) and (41) are re-injected into Eq. (24). This leads to a set of $N_m(2p+1)$ nonlinear equations given by

$$KB_0 = C_0 \quad (42)$$

$$\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 [1, p] \quad (43)$$

Coefficients $(C_0, (S_n, C_n)_{\forall n \in [1, p]})$ depend on coefficients $(B_0, (A_n, B_n)_{\forall n \in [1, p]})$. An extension of the classical Alternate Frequency-Time domain method (AFT-method Cameron and Griffin, 1989) is used to calculate these Fourier coefficients $(C_0, (S_n, C_n)_{\forall n \in [1, p]})$ (see Roncen et al. (2019a) and Roncen et al. (2019b) for more details).

Finally, the nonlinear Eqs. (42) and (43) are solved by minimizing the following relation

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

with

$$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 [1, p]} \right) \quad (45)$$

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

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

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

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

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

### 2.3.3. Initial validation based on a linear case

First of all, an initial validation of a linear case is undertaken in order to verify the relevance of the proposed approach in a simple case (for modeling the multipoint correlated random excitation presented in Section 2.2 for the mechanical problem presented in Section 2.1). In order to carry out this analysis, the nonlinear terms of Eq. (20) are considered equal to zero, and the associated dynamical linear equation projected on the first mode is given by

$$\dot{\omega}_1(t) + 2\xi_1 \Omega_1 \dot{\omega}_1(t) + \Omega_1^2 \omega_1(t) = \sum_{k=1}^{N_e} \Gamma_{1,k} f_k(t) \quad (49)$$

where the excitation terms $f_k$ are defined in Eq. (33).

The PSD of the numerical solution of Eq. (49) is compared to the traditional solution from linear vibration studies (Newland, 1975) which will therefore serve as a reference solution. This reference is provided by

$$S_{w_1}(\omega) = \sum_{i=1}^{N_e} \sum_{j=1}^{N_e} H_{1j}(\omega) H_{1j}^*(\omega) S_{ij}(\omega) \quad (50)$$

where $S_{w_1}(\omega)$ is the PSD of the modal displacement $w_1(t)$, $S_{ij}(\omega)$ are the terms of the matrix $S(\omega)$ in Eq. (38), and $H_{1k}(\omega)$ are the transfer functions defined by

$$H_{1k}(\omega) = \frac{\Gamma_{1,k}}{-\omega^2 + 2\xi_1 \Omega_1 j\omega + \Omega_1^2} \quad \forall k \in [1, N_e] \quad (51)$$

where $j$ is the unit imaginary number $\sqrt{-1}$, $\Gamma_{1,k}$, $\xi_1$, and $\Omega_1$ are defined in Section 2.1.

The comparison between the reference solution of the acceleration at the center of the beam and the numerical PSD of the solution of Eq. (49) (denoted by $\text{PSD}_{\text{ref}}$ and $\text{PSD}_{\text{num}}$, respectively) is given in Fig. 3. In this example the level of the excitation is chosen to be rectangular with a high level of excitation along the frequency band [80; 130] Hz (2.60 $10^{-6}$ N$^2$/Hz), and with a low level of excitation along the other frequency bands [20; 80] Hz and [130; 500] Hz (1.25 $10^{-10}$ N$^2$/Hz). It is observed that the proposed modeling of the multipoint correlated random excitation gives a result in perfect agreement with the reference solution. In order to quantify more precisely the adequacy of the proposed modeling, the following error $e_n$ is calculated

$$e_n = \frac{\sqrt{\int_{f_{\text{init}}}^{f_{\text{final}}} \text{PSD}_{\text{ref}}(f) df} - \sqrt{\int_{f_{\text{init}}}^{f_{\text{final}}} \text{PSD}_{\text{num}}(f) df}}{\sqrt{\int_{f_{\text{init}}}^{f_{\text{final}}} \text{PSD}_{\text{ref}}(f) df}} \quad (52)$$

where $f_{\text{init}}$ and $f_{\text{final}}$ correspond to the minimum and maximum frequency values in the chosen frequency band of interest (i.e., $f_{\text{init}} = 20$ Hz and $f_{\text{final}} = 500$ Hz in the present case under study). $e_n$ represents the error between both RMS accelerations at the center of the beam. With the red curve as reference one, $e_n$ is estimated at 1.20%, which unambiguously demonstrates the relevance of the proposed modeling for multipoint correlated random excitation (fully described by Shinozuka's work Shinozuka and Jan, 1972; Shinozuka, 1971; Shinozuka and Deodatis, 1991).

### 2.3.4. Nonlinear response of the beam to a multipoint correlated random excitation

In this section, the nonlinear vibrational response of the clamped-clamped beam subjected to a multipoint correlated random excitation is investigated. The numerical simulations are performed by using the previous modeling presented in Section 2.2. Results are given in Fig. 4 for three different levels of excitation. The first ideal input PSD of the excitation is depicted in Fig. 4(a) where level $S(\omega)$ (from Eq. (38)) is rectangular with a high level of excitation along the frequency band [80; 130] Hz (2.60 $10^{-6}$ N$^2$/Hz), and with a low level of excitation along the other frequency bands [20; 80] Hz and [130; 500] Hz (1.25 $10^{-10}$ N$^2$/Hz). The two other levels of excitation are multiples $2S(\omega)$ and $3S(\omega)$ of the first described ideal input PSD.

The acceleration output PSD at the center of the beam is observed in Fig. 4(b). Considering the comparison between the input PSD and the output PSD, there is one additional amplitude peak in the output PSD that is not present in the input PSD of the [130; 500] Hz band. This reveals nonlinear dynamic behavior of the beam subjected to a multipoint correlated random excitation. This additional peak corresponds to the third harmonic of the first vibration mode of beam $3f_1$, which is the result of the nonlinear cubic force $F_{NL}(t)$ restrained on the first vibration mode of the beam.

It is to be noted that these results are in line with those previously obtained in the study of a clamped-clamped beam subjected to one-point random excitations (Roncen et al., 2019a). Of course, they are not identical because the excitation studied is not the same (switch from a one-point random excitation to a more complex multipoint correlated random excitation), but the appearance of its third harmonic is visible in both cases and represents the main characteristic of the beam's nonlinear behavior. Therefore, this preliminary result shows that the proposed modeling of a multipoint correlated random excitation from
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Acceleration output PSD at the center of the beam from two methods: (red) reference solution and (blue) numerical solution; (b) Zoom along the frequency band [80; 130] Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. (a) Input PSD (b) Acceleration output PSD at the center of the beam. The level of the excitation is defined as $S(\omega)$ (blue), $2S(\omega)$ (green) and $3S(\omega)$ (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Section 2.2 can be adapted to the HBM for nonlinear vibration studies. It should be noted that this nonlinear response of the beam of the complete modeling of the multipoint correlated random excitation will be the reference nonlinear response for the rest of the study.

### 3. Advanced process based on the reduction of the number of excitation terms

Despite the use of Fast Fourier Transforms, the modeling of a multipoint correlated random excitation presented in Section 2.2 may be expensive from a numerical point of view if the number of excitation points becomes high. It may then be interesting or even necessary to find an alternative to this issue. In this section, two strategies will be discussed to reduce the number of excitation terms without losing much information on the excitation level or correlation. The first technique will be based on the well-known KL approach coupled with the Shinozuka decomposition. Then an original approach which consists of decomposing the information of the excitation (for both level and correlation) according to the eigenmodes of the mechanical system is proposed. Both techniques are presented in the rest of this section and adapted for predicting the nonlinear vibrations of the clamped-clamped beam subjected to one multipoint correlated random excitation.

#### 3.1. *Karhunen-Loève technique coupled with the Shinozuka decomposition*

The KL method allows the number of random variables to be reduced. It has already been applied to finite element models (for example the well-known contribution of Ghanem and Red-Horse Ghanem and Red-Horse, 1999). In these applications, a field of random variables is applied to a finite element model. For example, the value of a Young's modulus can vary randomly over the set of points in the structure and a statistical description of this random field is known. In particular, there is a spatial correlation such that two points in close proximity are likely to have the same Young's modulus value. The KL method then allows to reduce the number of random variables, rather than having as many random variables as finite element nodes. Somewhat similarly, in our case, it is the dynamic excitation vector that is random. The statistical properties of these signals are described in the spectral domain. For each excitation point, there are as many random variables as there are frequencies in the spectral decomposition of the signal (with the number of variables denoted by *p*). Multiplying by the number of excitation points gives $p \times N_e$ random variables. The KL method makes it possible to reduce the dimension $N_e$. To our knowledge, this use of the KL method, coupled with the Shinozuka decomposition is original.

In the following the KL modeling of a multipoint correlated random excitation is briefly described. For full details of the KL technique, the
<figure><img src="image_7.png" /><figcaption>Fig. 5. The first four KL eigenvectors of the multipoint correlated random excitation considered as a turbulent boundary layer noise (see Eq. (38)): 1st eigenvector $F_1$ (blue), 2nd eigenvector $F_2$ (red), 3rd eigenvector $F_3$ (green) and 4th eigenvector $F_4$ (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

interested reader is referred to Loeve (1977) and Kerschen et al. (2005). This section focuses on the case study described in Section 2.3. Contrary to Shinozuka's modeling developed in Section 2.2, the multipoint correlated random excitation modeled with the KL technique is not defined in the canonical basis of the excitation points. It is defined in the KL basis such as the temporal vector excitation $F(t)$ from Eq. (28) described as follows

$$F(t) = \sum_{k=1}^{N_e} \zeta_k(t) F_k \qquad (53)$$

where the variables $\zeta_k(t)$ and KL eigenvectors $F_k$ will be made explicit in the following. In the present study, the modeling of the multipoint random excitation is characterized by a cross-spectral density matrix $S(\omega)$ such as $S(\omega) = S(\omega)C$ (see Eq. (38)), in which cross-spectral correlation matrix $C$ is constant along the whole frequency band. Due to the Wiener-Khintchine relation, $C$ is Hermitian. Therefore, using the complex spectral theorem, $C$ has $N_e$ linearly independent eigenvectors ($F_k$) associated to $N_e$ real eigenvalues ($\lambda_k$) such as

$$\begin{bmatrix} C_{11} & C_{12} & \dots & C_{1N_e} \\ C_{21} & C_{22} & \dots & C_{2N_e} \\ \vdots & \vdots & \ddots & \vdots \\ C_{N_e 1} & C_{N_e 2} & \dots & C_{N_e N_e} \end{bmatrix} F_k = \lambda_k F_k \qquad (54)$$

It is to be noted that there are as many KL eigenvectors $F_k$ as the number of excitation terms $N_e$ in the definition of the multipoint excitation. The KL eigenvectors are constant along the frequency band since the cross-spectral correlation matrix $C$ is frequency-independent. The KL eigenvectors are chosen to be an orthonormal basis such as

$$F_i^T C F_j = \lambda_i \delta_{ij} \quad \forall (i, j) \in [1, N_e]^2 \qquad (55)$$

where $\delta_{ij}$ is the Kronecker symbol. The first four KL eigenvectors of the multipoint correlated random excitation considered as a turbulent boundary layer noise (see Eq. (38)) are plotted in Fig. 5. So, in the KL eigenvector basis, the cross-spectral density matrix $S(\omega)$ can be rewritten as

$$S^{\text{KL}}(\omega) = S(\omega) \begin{bmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \lambda_{N_e} \end{bmatrix} \qquad (56)$$

In this basis, the cross-spectral density matrix $S^{\text{KL}}(\omega)$ is diagonal with positive diagonal terms. Therefore, the Cholesky decomposition can be used again and the $S^{\text{KL}}(\omega)$ matrix can be decomposed as follows

$$S^{\text{KL}}(\omega) = \sqrt{S(\omega)} \begin{bmatrix} \sqrt{\lambda_1} & 0 & \dots & 0 \\ 0 & \sqrt{\lambda_2} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sqrt{\lambda_{N_e}} \end{bmatrix} \times \sqrt{S(\omega)} \begin{bmatrix} \sqrt{\lambda_1} & 0 & \dots & 0 \\ 0 & \sqrt{\lambda_2} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sqrt{\lambda_{N_e}} \end{bmatrix} \qquad (57)$$

As for Shinozuka's modeling developed in Section 2.2, each temporal force $\zeta_i(t)$ in the KL eigenvector basis is therefore defined by

$$\zeta_k(t) = \sqrt{2\Delta\omega} \left( \sum_{n=1}^{p} \sqrt{S(\omega_n)} \sqrt{\lambda_k} \cos(\omega_n t + \phi_{kn}) \right) \qquad (58)$$

The parameters $\Delta\omega$, $\omega_n$ and $\phi_{kn}$, present in the previous equation (58), are the same as the ones used in Section 2.2. Thanks to Eq. (53), temporal vector excitation $F(t)$ corresponds to the multipoint random excitation.

$$F(t) = \sqrt{2\Delta\omega} \sum_{k=1}^{N_e} \left( \sum_{n=1}^{p} \sqrt{S(\omega_n)} \sqrt{\lambda_k} \cos(\omega_n t + \phi_{kn}) \right) F_k \qquad (59)$$

The information on the multipoint correlated random excitation is present in variables $S(\omega_n)$, $\lambda_k$ and $F_k$ for each excitation point and for the whole frequency band. The KL technique is very interesting because it enables you to go from $N_e$ random variables which are dependent (see Eq. (33)) of each other to $N_e$ independent random variables.

As in Section 2.2, it is possible to rewrite Eq. (59) by using the Fast-Fourier Transform (FFT) (Vetterling et al., 1992). This leads to

$$F(t_n) = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{k=1}^{N_e} \sqrt{S(\omega_n)} \sqrt{\lambda_k} F_k \exp(j\phi_{kn}) \right) \right) \qquad (60)$$

It is to be noted that Eq. (60) represents another modeling of a multipoint correlated random excitation in comparison with Eq. (35). Each modeling uses the cross-spectral density matrix $S(\omega)$ differently, but the random part of the multipoint excitation is common and stems from independent random phases uniformly distributed between 0 and $2\pi$. By considering all the KL eigenvectors of matrix $C$, which is the reference model of this study, the excitation information (level and correlation) is complete. The linear validation given in Section 2.3.3 is the same for this modeling, including all KL eigenvectors.

Once the multipoint random excitation is defined in the KL eigenvectors basis, in order to calculate the nonlinear response of the system, it is necessary to project this excitation on the vibration modes of the beam ($Y_i(x)_{\forall i \in [1, N_m]}$). This projection is written as

$$\langle F(x, t_n), Y_i(x) \rangle = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{k=1}^{N_e} \sqrt{S(\omega_n)} \sqrt{\lambda_k} \langle F_k(x), Y_i(x) \rangle \exp(j\phi_{kn}) \right) \right) \qquad (61)$$

where $F_k(x)$ represents the interpolated KL eigenvector $F_k$.

Therefore, it leads to $N_e$ excitation terms associated with $N_m$ vibration modes of the beam system. In order to reduce the number of excitation terms, it may be interesting to estimate the influence of each KL eigenvector on each vibration mode of the system. In order to evaluate such a relation, the variable $a_{k,i}$ which represents the influence of the $k$th KL eigenvector on the $i$th vibration mode of the beam is defined as follows

$$a_{k,i} = \frac{\lambda_k \langle F_k(x), Y_i(x) \rangle^2}{\sum_{j=1}^{N_e} \lambda_j \langle F_j(x), Y_i(x) \rangle^2} \qquad (62)$$

Table 2 gives the impact of the first 10 KL eigenvectors on the first four vibration modes of the beam presented in Section 2.1. These results
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** (a) Acceleration output PSD at the center of the beam for reference model (solid lines) and the reduced model (dotted lines) with 1 KL eigenvector; (b) Zoom along the frequency band [320; 360] Hz. The level of the excitation is defined as $S(\omega)$ (blue), $2S(\omega)$ (green) and $3S(\omega)$ (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**
Influences of the first 10 KL eigenvectors on the first four vibration modes on the beam.</caption><thead><tr><th>KL eigenvector</th><th>F<sub>k</sub>(x)</th><th>λ<sub>k</sub></th><th>a<sub>k,1</sub> (%)</th><th>a<sub>k,2</sub> (%)</th><th>a<sub>k,3</sub> (%)</th><th>a<sub>k,4</sub> (%)</th></tr></thead><tbody><tr><td>F<sub>1</sub></td><td></td><td>18.91</td><td>95.93</td><td>–</td><td>6.43</td><td>–</td></tr><tr><td>F<sub>2</sub></td><td></td><td>15.78</td><td>–</td><td>91.97</td><td>–</td><td>12.72</td></tr><tr><td>F<sub>3</sub></td><td></td><td>12.25</td><td>3.97</td><td>–</td><td>81.90</td><td>–</td></tr><tr><td>F<sub>4</sub></td><td></td><td>9.24</td><td>–</td><td>7.77</td><td>–</td><td>72.24</td></tr><tr><td>F<sub>5</sub></td><td></td><td>6.96</td><td>0.09</td><td>–</td><td>11.20</td><td>–</td></tr><tr><td>F<sub>6</sub></td><td></td><td>5.31</td><td>–</td><td>0.23</td><td>–</td><td>14.33</td></tr><tr><td>F<sub>7</sub></td><td></td><td>4.13</td><td>0.007</td><td>–</td><td>0.40</td><td>–</td></tr><tr><td>F<sub>8</sub></td><td></td><td>3.27</td><td>–</td><td>0.02</td><td>–</td><td>0.59</td></tr><tr><td>F<sub>9</sub></td><td></td><td>2.64</td><td>0.001</td><td>–</td><td>0.05</td><td>–</td></tr><tr><td>F<sub>10</sub></td><td></td><td>2.17</td><td>–</td><td>0.004</td><td>–</td><td>0.08</td></tr></tbody></table>

show that only the first symmetric KL eigenvectors (the first anti-symmetric KL eigenvectors, respectively) are needed in order to keep the information on the excitation projected on the symmetric vibration modes of the beam (the anti-symmetric vibration modes, respectively). From the values depicted in Table 2, the multipoint correlated random excitation $F(x,t)$ restrained on the first KL eigenvector describes approximately 96% of the reference modal displacement of the beam’s first vibration mode, calculated with the reference multipoint correlated random excitation considered with all the KL eigenvectors. On the other hand, the multipoint correlated random excitation $F(x,t)$ restrained on the second KL eigenvector describes approximately 92% of the reference modal displacement of the beam’s second vibration mode, calculated with the reference multipoint correlated random excitation considered with all the KL eigenvectors. These results unambiguously demonstrate the relevance of the KL approach developed here and its effectiveness in reducing the number of excitation terms. In other words, the objective of reducing the number of excitation terms by restraining the number of KL eigenvectors in Eq. (58) is fully achieved.

In order to demonstrate the interest and effectiveness of the above, the nonlinear vibrations of the clamped-clamped beam, with a restrained number of KL eigenvectors in the expression of the multipoint correlated random excitation $F(x,t)$, is now studied. The multipoint correlated random excitation $F(x,t)$ projected on the vibration modes basis of the beam has the same formalism (see Eq. (79)) as the one used for the study of nonlinear vibrations of the clamped-clamped beam described in Section 2.3.4. The HBM applied to multipoint random excitations presented in Section 2.3.4 is therefore used to investigate the effect of the reduction of the number of excitation terms in order to predict the nonlinear vibrations of the beam subjected to multipoint random excitation.

The simulation of the nonlinear response of the reduced model by considering only one KL eigenvector in the description of the multipoint correlated random excitation is given in Fig. 6. It is compared to the nonlinear response simulated with the reference model of the excitation described in Section 2.2 (i.e., the original full representation of the multipoint correlated random excitation) for the same three different levels of excitation as presented in Section 2.3.4. It is to be noted that the first level of excitation $S(\omega)$ is chosen to be rectangular with a high level of excitation along the frequency band [80; 130] Hz ($2.60 \times 10^{-6} \text{ N}^2/\text{Hz}$), and with a low level of excitation along the other frequency bands [20; 80] Hz and [130; 500] Hz ($1.25 \times 10^{-10} \text{ N}^2/\text{Hz}$). It is quite clear that the reference and reduced models give the same nonlinear response for the beam system subjected to one multipoint correlated random excitation. The numerical error $e_n$ between the two simulation results is calculated for each level of excitation and it is given in Table 3. For this numerical error, the reference solution of the acceleration at the center of the beam and the numerical PSD (denoted by $PSD_{ref}$ and $PSD_{num}$) are respectively the one found with the complete modeling of the multipoint correlated excitation and the one found with the reduced model of the excitation with 1 KL eigenvector. Not surprisingly, the appearance of the third harmonic $3f_1$ is well predicted for both methods as depicted in Fig. 6(b). This demonstrates the feasibility and efficiency of predicting the nonlinear response of the beam system by applying the KL approach for reducing the number of excitation terms. It can be observed on Fig. 6(b) that there is a shift in the natural frequency of the first mode as a function of the excitation level. The peak of the third harmonic is situated at $F = 339.112 \text{ Hz}$ for the first level of excitation $S(\omega)$, then it goes to $F = 339.354 \text{ Hz}$ for the second level of excitation $2S(\omega)$ and reaches $F = 339.840 \text{ Hz}$ for the third level of excitation $3S(\omega)$. This means that the frequency of the beam’s first vibration mode increases with the excitation level. Respectively, the shift in frequency is represented by frequencies $f_1 = 113.037 \text{ Hz}$ (for $S(\omega)$), $f_1 = 113.118 \text{ Hz}$ (for $2S(\omega)$) and $f_1 = 113.280 \text{ Hz}$ (for $3S(\omega)$).

Whereas the reference model given in Section 2.2 needs $N_e$ excitations terms in order to describe the excitation, only one excitation term, corresponding to one KL eigenvector, has to be used in order to study the nonlinear vibrations of the beam restrained on its first vibration mode. In order to study nonlinear vibrations of the beam with more restrained vibration modes, more KL eigenvectors have to be kept but this number should always be more interesting than the number $N_e$ of excitation terms in the reference model, especially if $N_e$ becomes big. It is worth remembering that $N_e$ defines the number of punctual correlated temporal forces distributed along the mechanical
system under study. In the general context of an engineering problem, it is not uncommon to be confronted with excitations distributed over a non-negligible surface of the mechanical system under study and consequently to end up with a non-negligible number $N_e$. It is therefore obvious that the KL technique could also be interesting in a more generic context and for other engineering applications. However, some limitations must be stated: this technique has shown its efficiency if the correlation of the excitation is constant along the frequency band. If the correlation function becomes frequency-dependent, as is the case for a diffuse field acoustic load applied to aeronautical structures, then the KL eigenvectors ($F_k(\omega)_{\forall k \in \{1, N_e\}}$ associated to the KL eigenvalues ($\lambda_k(\omega)_{\forall k \in \{1, N_e\}}$ are also frequency-dependent. Then, the computing of the reduced multipoint correlated random excitation with the KL technique is no longer interesting for large frequency bands (or/and for calculations with a small frequency discretization), where the KL eigenproblem needs to be solved at each frequency step. Adapted to nonlinear simulations with the use of the HBM to random excitation (see Section 2.3.2), the minimization of Eq. (24) needs to be computed for the whole frequency band. So, all the KL eigenvectors, calculated to each frequency step, need to be stored before the minimization of the whole nonlinear resolution. The KL technique does not seem to be interesting for the simulation of nonlinear vibrations of structures subjected to multipoint correlated random excitations for which the correlation description is frequency-dependent.

## 3.2. Projection on the eigenmodes of the structure

In this section, an original technique which enables to reduce the number of excitation terms while remaining in a formalism adapted for the use of the HBM presented in Section 2.3.4 is discussed. This proposed method is based on the projection of the cross-spectral density matrix $S(\omega)$, on the vibration modes basis of the system $[\mathbf{Y}_1, \dots, \mathbf{Y}_{N_m}]$. The vibration modes have to be normalized, with respect to the mass matrix of the system. As for the KL technique presented in Section 3.1, the temporal vector excitation $F(t)$ from Eq. (28) is now defined as

$$F(t) = \sum_{k=1}^{N_m} y_k(t) \mathbf{Y}_k \quad (63)$$

where the variables $y_k(t)$ will be defined later. $\mathbf{Y}_k$ correspond to the vibration modes of the system reduced to the $N_e$ degrees of freedom of the domain $\Omega$. Therefore, the modeling of the multipoint correlated random excitation given here is defined in the vibration modes basis of the system.

The projection of the cross-spectral density matrix $S(\omega)$ on the vibration modes basis of the system $[\mathbf{Y}_1, \dots, \mathbf{Y}_{N_m}]$ leads to an another matrix $S^{\text{proj}}(\omega)$ defined by

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

with

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

$$\mathbf{P} = \begin{bmatrix} \mathbf{Y}_1 & 0 & 0 \\ 0 & \mathbf{Y}_2 & \vdots \\ \vdots & \vdots & \vdots \\ 0 & \dots & 0 \end{bmatrix} \quad (66)$$

Therefore, matrix $S^{\text{proj}}(\omega)$ is described by

$$S^{\text{proj}}(\omega) = \begin{bmatrix} \mathbf{Y}_1^T S(\omega) \mathbf{Y}_1 & \mathbf{Y}_1^T S(\omega) \mathbf{Y}_2 & \dots & \mathbf{Y}_1^T S(\omega) \mathbf{Y}_{N_m} \\ \mathbf{Y}_2^T S(\omega) \mathbf{Y}_1 & \mathbf{Y}_2^T S(\omega) \mathbf{Y}_2 & \dots & \mathbf{Y}_2^T S(\omega) \mathbf{Y}_{N_m} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{Y}_{N_m}^T S(\omega) \mathbf{Y}_1 & \mathbf{Y}_{N_m}^T S(\omega) \mathbf{Y}_2 & \dots & \mathbf{Y}_{N_m}^T S(\omega) \mathbf{Y}_{N_m} \end{bmatrix} \quad (67)$$

At this stage, it is necessary to demonstrate that matrix $S^{\text{proj}}(\omega)$ is positive semi-definite thanks to the following mathematical property

$$A \text{ is positive semi-definite} \Leftrightarrow \mathbf{x}^T A \mathbf{x} \ge 0 \quad \forall \mathbf{x} \neq \mathbf{0} \quad (68)$$

The following discussion responds to this expectation. First, the matrix $\mathcal{S}(\omega)$ can be rewritten

$$\mathcal{S}(\omega) = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} [\mathcal{S}(\omega)] \begin{bmatrix} 1 & 1 & \dots & 1 \end{bmatrix} \quad (69)$$

Now, let us call $\mathbf{x}$ a vector of dimensions $N_m N_e \times 1$ where at least one of its coefficients (each of them is a vector of dimensions $N_e \times 1$) is not equal to zero. $\mathbf{x}$ is defined as

$$\mathbf{x} = \begin{bmatrix} \mathbf{x}_1 \\ \vdots \\ \mathbf{x}_{N_m} \end{bmatrix} \quad (70)$$

Thanks to Eq. (69), the following relation can be written

$$\mathbf{x}^T \mathcal{S}(\omega) \mathbf{x} = \left(\sum_{i=1}^{N_m} \mathbf{x}_i\right)^T \mathcal{S}(\omega) \left(\sum_{i=1}^{N_m} \mathbf{x}_i\right) \quad (71)$$

Since at least one of the vector $\mathbf{x}$ coefficients is not equal to zero, one of the $\mathbf{x}_i$ vectors is not equal to zero, then $(\sum_{i=1}^{N_m} \mathbf{x}_i)$ is not equal to zero. Thanks to the mathematical property given by Eq. (68) and since $S(\omega)$ is positive definite, then $\mathbf{x}^T S(\omega) \mathbf{x} \ge 0$. Therefore, the matrix $\mathcal{S}(\omega)$ is positive semi-definite. Now, let us apply the same method to demonstrate that matrix $S^{\text{proj}}(\omega)$ is also positive semi-definite. Let us call $\mathbf{z}$ a vector of dimensions $N_m N_e \times 1$, where at least one of its coefficients (each of them is a vector of dimensions $N_e \times 1$) is not equal to zero

$$\mathbf{z} = \begin{bmatrix} \mathbf{z}_1 \\ \vdots \\ \mathbf{z}_{N_m} \end{bmatrix} \quad (72)$$

Thanks to Eq. (64), the following formula can be specified

$$\mathbf{z}^T S^{\text{proj}}(\omega) \mathbf{z} = \mathbf{z}^T \mathbf{P}^T \mathcal{S}(\omega) \mathbf{P} \mathbf{z} = (\mathbf{P} \mathbf{z})^T \mathcal{S}(\omega) (\mathbf{P} \mathbf{z}) \quad (73)$$

Since at least one of the coefficients of vector $\mathbf{z}$ is not equal to zero and $\mathbf{P}$ corresponds to a modal matrix (vibrations modes of a non-rigid structure), then one of the coefficients of the vector $(\mathbf{P} \mathbf{z})$ is not equal to zero. Thanks to the mathematical property given by Eq. (68) and since $\mathcal{S}(\omega)$ is positive semi-definite, then $\mathbf{z}^T S^{\text{proj}}(\omega) \mathbf{z} \ge 0$. This demonstrates that matrix $S^{\text{proj}}(\omega)$ is positive semi-definite.

Then, using the Cholesky decomposition as follows

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

As in Section 2.2, $\mathbf{G}(\omega)$ is a lower triangular matrix given by

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

As previously done in Section 2.2 for the complete modeling of a multipoint correlated random excitation, the $G_{ij}(\omega)$ coefficients can be determined with the principal minors of matrix $S^{\text{proj}}(\omega)$. This method allows you to consider a matrix $S^{\text{proj}}(\omega)$ of dimensions $N_m \times N_m$ which is more interesting than considering matrix $S(\omega)$ of dimensions $N_e \times N_e$, especially if $N_e$ becomes big. Therefore, this method defines $N_m$ excitations which are modal and representative of the complete multipoint correlated random excitation. Then the temporal modal excitation $(y_k(t))_{\forall k \in \{1, N_m\}}$ is given by

$$y_k(t) = \sqrt{2\Delta\omega} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |G_{kl}(\omega_n)| \cos(\omega_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) \quad (76)$$
where variables $p, \Delta\omega, \omega_n, \theta_{kl}(\omega_n)$ and $\phi_{ln}$ are the same as the ones used in Section 2.2. Thanks to Eq. (63), the temporal vector excitation $F(t)$ corresponding to the multipoint random excitation can be written from this modeling as

$$F(t) = \sqrt{2\Delta\omega} \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} \sum_{n=1}^{p} |G_{kl}(\omega_n)| \cos(\omega_n t + \theta_{kl}(\omega_n) + \phi_{ln}) \right) \mathbf{Y}_k \quad (77)$$

As $p$ is considered as a power of 2 and as the frequency discretization is correctly chosen, Eq. (77) can be rewritten by using the Fast-Fourier Transform (FFT)

$$F(t_n) = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} G_{kl}(\omega_n) \exp(j\phi_{ln}) \right) \mathbf{Y}_k \right) \right) \quad (78)$$

where IFFT denotes the Inverse of the Fast Fourier Transform and $j$ is the unit imaginary number $\sqrt{-1}$. $t_n$ is the discretized time $t_n = n \cdot dt$ where $dt$ is the temporal sampling step. It should be noted with this modeling choice, that $N_m$ random variables are dependent of each other, but there are much less random variables than within the full description of the multipoint random excitation (see Section 2.2), which could be an interesting contribution if the methodology is to be adapted to more complex systems.

As the Karhunen-Loève method decomposes the random field into independent random variables, we find in Eq. (60) a simple sum over the number of Karhunen-Loève terms ($N_e$). The modal method involves correlated random variables, which generates a double sum $\sum_{k=1}^{N_m} \sum_{l=1}^{k}$ in Eq. (78). Nevertheless, for the solving of the nonlinear problem, the excitation term must be projected on each eigenmode. This projection is detailed by the following equation

$$\langle F(x, t_n), Y_i(x) \rangle = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{k=1}^{N_m} \left( \sum_{l=1}^{k} G_{kl}(\omega_n) \exp(j\phi_{ln}) \right) \right) \langle Y_k(x), Y_i(x) \rangle \right) \quad (79)$$

The orthogonality of the modal basis allows to suppress the summation $\sum_{k=1}^{N_m} \sum_{l=1}^{k}$ in Eq. (79). This orthogonality is relative to the scalar product associated with the mass operator in the general case. In the chosen example, the linear mass is constant along the beam which leads to

$$\langle F(x, t_n), Y_i(x) \rangle = \sqrt{2\Delta\omega} \mathfrak{R} \left( \text{IFFT} \left( \sum_{l=1}^{i} G_{il}(\omega_n) \exp(j\phi_{ln}) \right) \langle Y_i(x), Y_i(x) \rangle \right) \quad (80)$$

With the Karhunen-Loève method, as many projections as there are structural modes must be performed, whereas with the modal method, the excitation is already projected on the basis of the eigenmodes. Therefore, the computational time advantage brought by the decorrelation of the random variables for the Karhunen-Loève method is balanced by the absence of additional computations for the modal projection in the case of the modal method.

In order to assess the validity of the above proposal, the simulation of the nonlinear response of the beam with only one vibration mode in the description of the multipoint correlated random excitation is given in Fig. 7. Comparison with the reference model of the excitation described in Section 2.2 is also provided for the same three different levels of excitation as shown in Section 2.3.4. The first level of the excitation $S(\omega)$ is chosen to be rectangular with a high level of excitation along the frequency band [80; 130] Hz ($2.60 \cdot 10^{-6} \text{ N}^2/\text{Hz}$), and with a low level of excitation along the other frequency bands [20; 80] Hz and [130; 500] Hz ($1.25 \cdot 10^{-10} \text{ N}^2/\text{Hz}$). Numerical error $e_n$ between the two simulation results is calculated for each level of excitation and given in Table 3. For this numerical error, the reference solution of the acceleration at the center of the beam and the numerical PSD (denoted by PSD<sub>ref</sub> and PSD<sub>num</sub>) are respectively the one found with the complete modeling of the multipoint correlated excitation and the one found with the reduced

<table><caption>**Table 3**
Values of the numerical errors between two simulation results (the simulation used for the reference is the one introduced in Section 2.2).</caption><thead><tr><th>Numerical error e<sub>n</sub> (%)</th><th>S(ω)</th><th>2S(ω)</th><th>3S(ω)</th></tr></thead><tbody><tr><td>Karhunen-Loève decomposition</td><td>0.83</td><td>0.23</td><td>0.64</td></tr><tr><td>Proposed approach</td><td>0.18</td><td>1.10</td><td>0.59</td></tr></tbody></table>

model of the excitation with one vibration mode in the projection of the multipoint correlated random excitation. As previously seen for the KL method, the appearance of the third harmonic $3f_1$ is well predicted as depicted in Fig. 7(b). This shows without any ambiguity the value of the original technique proposed in order to reduce the number of excitation terms, while allowing a correct prediction of the nonlinear signature of the beam system. As described in Section 3.1, there is the same shift in frequency for the first vibration mode of the beam (see Fig. 7(b)) as a function of the level of excitation.

Calculations made in this study were based on a CentOS Linux 7 with a Processor Intel Core i7-9700 K CPU, 8 cores, 3.60 GHz and CPU 32 Go RAM. The average calculation time for computing the temporal vector excitation $F(t)$ is given for each modeling in Table 4. It is clearly shown that the two techniques presented in this study are much more interesting than the complete modeling of the multipoint correlated random excitation. It is to be noted that each modeling of the excitation gives the same calculation time for the resolution of the nonlinear problem. Indeed, each modeling proposes approximately the same level of excitation. The evolution of the nonlinear response will therefore depend on the nonlinear function defined in the mechanical problem. For the problem proposed in this study, the nonlinearity is geometrical and involves each beam vibration mode (see Eq. (22)). Therefore, no difference has been noted in the research for the nonlinear response of the beam, restrained on its first vibration mode, with the three different models of the multipoint correlated excitation.

Whereas the reference model given in Section 2.2 needs $N_e$ excitations terms in order to describe the excitation, only one excitation term, corresponding to one vibration mode, has to be used in order to predict the nonlinear response of the beam system restrained on its first vibration mode.

For the interested reader it is also important to keep in mind that, in the case of a multipoint correlated random excitation for which the correlation function is frequency independent, the proposed method is indeed as efficient as the KL decomposition method. To be noted that the calculation of the projection on the structural eigenmodes of the system is not expensive since the definition of the correlation function and the definition of the structural eigenmodes are analytical and are known as a function of the abscissa. The projection is therefore simple to calculate even if there are more structural eigenmodes to be considered. Thus, even if the correlation information depends on the frequency, the projection basis of this modeling (the vibration modes basis) is not frequency-dependent. Therefore, this proposed original method seems to be more convenient to use than the KL decomposition, in order to study nonlinear vibrations of mechanical systems subjected to correlated excitations for which the correlation depends on the frequency. Indeed, for each frequency in the spectral decomposition, the KL method has to solve and store the eigenvalues and eigenvectors of the spectral density matrix (which is a very expensive step from a numerical point of view if the number of excitation terms is large). On the contrary, the proposed technique will require computing the matrix $S^{\text{proj}}(\omega)$ only once, which can provide a significant benefit when solving finite element problems for industrial applications.

Even if it is obvious when reading and comparing the two formalisms, this latter benefit will need to be validated and verified for more complex mechanical systems. This point is outside the scope of the present study and should be the subject of future works.
<figure><img src="image_9.png" /><figcaption>Fig. 7. (a) Acceleration output PSD at the center of the beam for the reference model (solid lines) and the proposed reduced model (crosses) with 1 vibration mode in the projection of the multipoint correlated random excitation; (b) Zoom along the frequency band [320; 360] Hz. The level of the excitation is defined as S(ω) (blue), 2S(ω) (green) and 3S(ω) (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 4**
Average calculation time to compute the temporal vector excitation F(t) for each modeling (the average is made with 100 calculations).</caption><thead><tr><td>Strategy</td><td>Average calculation time (s)</td></tr></thead><tbody><tr><td>Complete modeling (Eq. (35))</td><td>10.55</td></tr><tr><td>Karhunen-Loève decomposition (Eq. (60))</td><td>0.08</td></tr><tr><td>Proposed approach (Eq. (78))</td><td>0.07</td></tr></tbody></table>

## 4. Conclusion

This paper investigates the nonlinear simulations on a clamped-clamped beam solicited by multipoint random excitation. More specifically one of the main contributions of the proposed study is to reduce the number of excitation points in order to reduce calculation time and data storage, while maintaining an accurate prediction of the system's nonlinear response.

Therefore, the KL decomposition method is initially applied in order to reduce the number of excitation terms. With one vibration mode, the reduction of the number of excitation terms is very interesting and allows you to find the same nonlinear signature as for complete modeling of multipoint random excitation from Shinozuka's work (Shinozuka and Jan, 1972; Shinozuka, 1971; Shinozuka and Deodatis, 1991).

Moreover, a second original method is proposed in order to reduce the number of excitation terms while remaining in a formalism of the HBM applied to random excitations. It is to be noted that this original method is also compatible with nonlinear methods based on the temporal integration of the nonlinear equation. This strategy is based on the projection of excitation on the vibration modes basis of the system. Thanks to this method, the projection basis is frequency-independent, so this method can be easily extended to more complex excitations for which the correlation information depends on the frequency. This proposed original numerical method appears to be an alternative to the KL decomposition method, in order to study and simulate the nonlinear vibratory behavior of a complex mechanical system subjected to random multipoint excitations. This will be the subject of future works.

### CRediT authorship contribution statement

**S. Talik:** Methodology, Software, Validation, Investigation, Writing – original draft, Writing – review & editing, Visualization. **M. Claeys:** Conceptualization, Methodology, Validation, Investigation, Resources,

Writing – original draft, Writing – review & editing, Project administration, Funding acquisition. **J.-J. Sinou:** Conceptualization, Methodology, Validation, Investigation, Writing – original draft, Writing – review & editing, Supervision, Project administration. **J.-P. Lambelin:** Conceptualization, Methodology, Validation, Investigation, Resources, Writing – review & editing, Project administration, Funding acquisition.

### Declaration of competing interest

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

### Data availability

The authors are unable or have chosen not to specify which data has been used.

### Acknowledgments

J.-J. Sinou acknowledges the support of the Institut Universitaire de France.

### References

Bellizzi, S., Bouc, R., 1999. Analysis of multi-degree of freedom strongly non-linear mechanical systems with random input: part I: non-linear modes and stochastic averaging. Probab. Eng. Mech. 14 (3), 229-244.

Cameron, T., Griffin, J., 1989. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. J. Appl. Mech. 56 (1), 149-154.

Claeys, M., Lambelin, J.-P., Chantereau, Y., Sinou, J.-J., 2019. Dataset of multi-harmonic measurements for the experimental CEA-beam benchmark structure. Data Brief 27, 104563.

Claeys, M., Sinou, J.-J., Lambelin, J.-P., Alcoverro, B., 2014. Multi-harmonic measurements and numerical simulations of nonlinear vibrations of a beam with non-ideal boundary conditions. Commun. Nonlinear Sci. Numer. Simul. 19 (12), 4196-4212.

Claeys, M., Sinou, J.-J., Lambelin, J.-P., Todeschini, R., 2016a. Experiments and numerical simulations of nonlinear vibration responses of an assembly with friction joints - Application on a test structure named "Harmony". Mech. Syst. Signal Process. 70-71, 1097-1116.

Claeys, M., Sinou, J.-J., Lambelin, J.-P., Todeschini, R., 2016b. Modal interactions due to friction in the nonlinear vibration response of the "Harmony" test structure: experiments and simulations. J. Sound Vib. 376, 131-148.

Efimsov, B.M., 1982. Characteristics of the field of turbulent wall pressure fluctuations at large Reynolds numbers. Sov. Phys. Acoust. 28, 289-292.

Eigoli, A., Ahmadian, M., 2011. Nonlinear vibration of beams under nonideal boundary conditions. Acta Mech. 218, 259-267.
Fang, J., Elishakoff, I., Caimi, R., 1995. Nonlinear response of a beam under stationary random excitation by improved stochastic linearization method. Appl. Math. Model. 19 (2), 106-111.

Ghanem, R., Red-Horse, J., 1999. Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach. Physica D 133, 137-144.

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

Huang, Z., Wang, Y., Zhu, W., Huang, Z., 2020. Deterministic and random response evaluation of a straight beam with nonlinear boundary conditions. J. Vib. Eng. Technol. 8.

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

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

Kerschen, G., Golinval, J.-C., Vakakis, A.F., Bergman, L.A., 2005. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview. Nonlinear Dynam. 41, 147-169.

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

Kovacic, I., Brennan, M., 2011. The Duffing Equation. John Wiley & Sons.

Loeve, M., 1977. Probability Theory. Springer-Verlag.

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

Nayfeh, A., 1993. Introduction to Perturbation Techniques. In: Wiley Classics Library, J. Wiley, New York.

Nayfeh, A., Mook, D., 1979. Nonlinear Oscillations. John Wiley and Sons, New York.

Newland, D.E., 1975. An Introduction to Random Vibrations and Spectral Analysis.

Noël, J., Kerschen, G., 2017. Nonlinear system identification in structural dynamics: 10 more years of progress. Mech. Syst. Signal Process. 83, 2-35.

Poirel, D., Price, S., 2007. Bifurcation characteristics of a two-dimensional structurally non-linear airfoil in turbulent flow. Comput. Struct. 48 (4), 423-435.

Ribeiro, P., 2004. Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods. Comput. Struct. 82 (17), 1413-1423.

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

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

Roncen, T., Sinou, J.-J., Lambelin, J.-P., 2019c. Experiments and simulations of an industrial assembly with different types of nonlinear joints subjected to harmonic vibrations. J. Sound Vib. 376, 458-478.

Roncen, T., Sinou, J.-J., Lambelin, J.-P., 2021. Experiments and simulations of the structure Harmony-Gamma subjected to broadband random vibrations. Mech. Syst. Signal Process. 159, 107849.

Sayed, M., Mousa, A., Mustafa, I., 2020. Stability and bifurcation analysis of a buckled beam via active control. Appl. Math. Model. 82, 649-665.

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

Shinozuka, M., Deodatis, G., 1991. Simulation of stochastic processes by spectral representation. Appl. Mech. Rev. 44 (4), 191-204.

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

Singh, P., Yeong, H., Zhang, H., Rapti, Z., 2016. Stochastic stability and dynamics of a two-dimensional structurally nonlinear airfoil in turbulent flow. Meccanica 51 (11), 2665-2688.

Talik, S., Sinou, J.-J., Claeys, M., Lambelin, J.-P., 2022. Nonlinear vibrations of a beam with non-ideal boundary conditions and subjected to two correlated or uncorrelated broadband random excitations - experiments, modeling and simulations. Commun. Nonlinear Sci. Numer. Simul. 110, 106328.

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

Xu, Y., Li, Y., Xu, 2016. A method to stochastic dynamical systems with strong nonlinearity and fractional damping. Nonlinear Dynam. 83, 2311-2321.

Xu, Y., Liu, Q., Guo, G., Xu, C., Liu, D., 2017. Dynamical responses of airfoil models with harmonic excitation under uncertain disturbance. Nonlinear Dynam. 89, 1579-1590.

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