# Analysis and Physical Interpretation of the Uncertainty Effect in Structural Dynamics

Louis Jezequel, Hugo de Filippis, Alexy Mercier

► To cite this version:

Louis Jezequel, Hugo de Filippis, Alexy Mercier. Analysis and Physical Interpretation of the Uncertainty Effect in Structural Dynamics. Journal of Engineering Mechanics - ASCE, 2022, 148 (2), 10.1061/(ASCE)EM.1943-7889.0002040. hal-04083030

HAL Id: hal-04083030

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

Submitted on 21 Mar 2025

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
# Analysis and Physical Interpretation of the Uncertainty Effect in Structural Dynamics

Louis Jézéquel<fnref n="1" />; Hugo de Filippis<fnref n="2" />; and Alexy Mercier<fnref n="3" />

**Abstract:** In order to determine the steady-state response of random dynamical systems, the polynomial chaos expansion (PCE) approach appears to be the most suitable, but its convergence is weak around the eigenfrequencies. Indeed, increasing the order of the polynomial creates a side effect, which decreases the convergence of the expansion. Previous studies used convergence accelerating methods to reduce this side effect and the main aim of this study is to present an original approach to obtain a closed-form solution for a class of random dynamical system. The analytical developments are based on the introduction of a new class of modes called eigenmodes of perturbation and on the resolution of periodic systems. This approach gives a physical interpretation of the damping effect observed on the mean response of linear random dynamical systems by referring to the propagation of waves in a fictitious semiperiodic system associated with the eigenmodes of perturbation. The dynamic models used to expose the method are simple but the study is extended to complex structures using modal synthesis. DOI: 10.1061/(ASCE)EM.1943-7889.0002040. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.

## Introduction

In many scientific fields, determining the random response of a dynamic system is a fundamental problem that necessitates specific methods. The Monte Carlo simulation, which has become a reference, is a well-known method despite its main drawback of requiring high computing performance on complex systems. That is why other methods have been developed and used instead of the latter, like the disturbance method, which used the expansion of random variables into Taylor series (Kleiber and Hien 1992); the Neumann method, which used the Neumann series expansion (Yamazaki et al. 1988); and the Karhunen-Loève expansion (Huang et al. 2001).

Moreover, the polynomial chaos expansion (PCE), originally developed by Wiener (1938), regained attention in the early 1990s following the publication of Ghanem and Spanos' book (Ghanem and Spanos 1991), and is widely adopted nowadays. The polynomials used in the PCE are hypergeometrical polynomials according to the Askey scheme (Xiu and Karniadakis 2002). Among the different polynomials that exist, the most commonly used are the Hermite polynomials, which have an orthogonal property with respect the Gaussian probability density functions (PDF), the Legendre polynomials with a uniform PDF, the Laguerre polynomials with a gamma PDF, and the Jacobi polynomials with a beta PDF.

Nevertheless, this method also includes drawbacks, namely using polynomials to expand the response generates undesirable

interferences around the eigenfrequencies, which decrease by increasing the order of the polynomials, but this slows down the convergence (Jacquelin et al. 2015a).

To remedy this issue, convergence accelerating methods have been experimented with, like the Aitken and Shank transformations (Keshavarzzadeh et al. 2014; Jacquelin et al. 2015b), but they still need high-order polynomials. Some researchers have used the Padé approximants method (Brezinski 1996; Baker and Graves-Morris 1996), which expands the dynamical random response on a fraction with polynomials in both the numerator and the denominator. A generalization of the Padé approximants has also been used in a convergence study (Matos 1996; Jacquelin et al. 2017). In this paper, the Chebyshev polynomials are used to obtain a closed-form of steady-state response of a class of random dynamic systems. Note that certain articles have focused on the use of Chebyshev polynomials in different subjects, like the analysis of bifurcation in stochastic systems (Ma et al. 2006; Liu et al. 2014), for solving variable order fractional differential-integral equations (Liu et al. 2016), for estimating geotechnical shear strength parameters (Li et al. 2006), and for a nonparametric identification procedure (Argoul and Jezequel 1989).

After the presentation of the specificity of using Chebyshev polynomials of the second kind in the context of the PCE method, a new class of modes, the so-called eigenmode of perturbation, is presented in the second part of this study. These modes allow expression of the steady-state response of a mechanical structure disturbed by a random stiffness parameter. Thanks to these early developments, the analytical computation of the limit of the polynomial chaos expansion was then obtained. Simple examples are used to illustrate the method and interpret the solution by wave analysis in a periodic fictitious space associated with the random parameter. In the last part, the method proposed is extended to complex structures using modal synthesis.

## Expansion on Chebyshev Polynomials

The system studied is a linear random dynamical system with $N$ degrees of freedom (DOF), excited with a harmonic force vector $F$ of frequency $\omega$. The system is defined by invariant mass and

<footnote n="1">¹University Professor, Laboratoire de Tribologie et Dynamique des Systèmes, École Centrale de Lyon, Écully 69134, France. Email: louis.jezequel@ec-lyon.fr</footnote>

<footnote n="2">²Student, École Nationale des Travaux, Publics de l'État, Vaulx-en-Velin 69120, France. Email: hugo.defilippis@entpe.fr</footnote>

<footnote n="3">³Ph.D. Student, Laboratoire de Tribologie et Dynamique des Systèmes, École Centrale de Lyon, Écully 69134, France (corresponding author). Email: alexy.mercier@ec-lyon.fr</footnote>

<footnote>Note. This manuscript was submitted on December 16, 2020; approved on September 8, 2021; published online on November 30, 2021. Discussion period open until April 30, 2022; separate discussions must be submitted for individual papers. This paper is part of the *Journal of Engineering Mechanics*, © ASCE, ISSN 0733-9399.</footnote>
stiffness matrices, written as **M**, **K**, respectively, and a random stiffness matrix **A** associated with an uncertain parameter x, which is a zero-mean random variable. The sources of uncertainties modeled by the random parameter x and the matrix **A** can come from the elastic parameters of a material constituting a part of the structure or from imprecise boundary conditions. The steady-state response vector **U** is the solution of the following matrix system of N × N dimension:

$$(-\omega^2\mathbf{M} + \mathbf{K} + x\mathbf{A})\mathbf{U}(\omega, x) = \mathbf{F}(\omega) \quad (1)$$

with $x \in [-1, 1]$.

If the notation $\tilde{\mathbf{K}} = \mathbf{K} - \omega^2 \mathbf{M}$ is adopted and if $\omega$ are not explicitly noted in the vectors, the governing equation becomes

$$(\tilde{\mathbf{K}} + x\mathbf{A})\mathbf{U}(x) = \mathbf{F} \quad (2)$$

According to the polynomial chaos expansion, the random steady-state dynamical response of the system is expressed as follows:

$$\mathbf{U}(x) = \sum_{n=0}^{\infty} \mathbf{U}_n P_n(x) \quad (3)$$

In this study, the Chebyshev polynomials of the second kind will be used considering their specific properties compared to the other possible polynomials (Mason 1993). In particular, the Chebyshev polynomials satisfy the following recurrence equation:

$$P_{n+1}(x) - 2xP_n(x) + P_{n-1}(x) = 0 \quad (4)$$

Note that this recurrence equation has constant coefficients, contrary to the Legendre and Laguerre polynomials, which verify the following relationships:

**Legendre :**

$$\forall n \ge 1 \quad (n+1)L_{n+1}(x) - (2n+1)xL_n(x) + nL_{n-1}(x) = 0 \quad (5)$$

**Laguerre :**

$$\forall n \ge 1 \quad (n+1)L_{n+1}(x) + (x-2n-1)L_n(x) + nL_{n-1}(x) = 0 \quad (6)$$

The probability density function associated with the Chebyshev polynomials of the second kind is given by

$$f(x) = \frac{2}{\pi} \sqrt{1 - x^2} \quad \text{on } [-1; 1] \\ \langle P_m(x), P_n(x) \rangle = \int_{-1}^{1} P_m(x) P_n(x) f(x) dx = \begin{cases} 1 & \text{if } m = n \\ 0 & \text{if } m \neq n \end{cases} \quad (7)$$

A shown in Fig. 1, this probability density function seems to be well adopted to describe the uncertainties in the product design and process.

Using Eqs. (3) and (4), the random dynamical response of the system from Eq. (1) at the n order becomes

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Probability density function for the Chebyshev polynomials of the second kind.</figcaption></figure>

$$\begin{bmatrix} \mathbf{K} - \omega^2 \mathbf{M} & 0 & \cdots & 0 \\ 0 & \mathbf{K} - \omega^2 \mathbf{M} & & \vdots \\ \vdots & & \mathbf{K} - \omega^2 \mathbf{M} & 0 \\ 0 & \cdots & 0 & \mathbf{K} - \omega^2 \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{U}_0 \\ \vdots \\ \mathbf{U}_n \end{bmatrix} + \frac{1}{2} \begin{bmatrix} \mathbf{A} & 0 & \cdots & 0 \\ \mathbf{A} & 0 & \cdots & 0 \\ \vdots & \mathbf{A} & 0 & \cdots \\ 0 & \mathbf{A} & 0 & \cdots \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ 0 \\ \vdots \\ 0 \end{bmatrix} \quad (8)$$

which can be associated with the following recurrence equation:

$$\forall n \ge 1 \quad (\mathbf{K} - \omega^2 \mathbf{M})\mathbf{U}_n + \frac{1}{2}\mathbf{A}(\mathbf{U}_{n-1} + \mathbf{U}_{n+1}) = 0 \quad (9)$$

As already mentioned (Jacquelin et al. 2015a), the convergence of PCE near the resonance frequencies is weak. Thus, for low damping systems, the truncation to the order n of PCE induces n response peaks at each resonance frequency. To illustrate this point, the absolute value of the mean response $\mathbf{U}_0$ for the two degrees of freedom (2DOF) system defined in Fig. 2 has been calculated in the case $(P/k) = 0.1$ and $\xi = (c/2\sqrt{km}) = 0.02$. The darker curves shown in Fig. 3 were obtained with n = 3 and the lighter curves are the exact solutions corresponding to n = $\infty$. The analytical approach presented in this study is aimed at obtaining this exact solution in the following form:

$$\forall n \ge 1 \quad \mathbf{U}_n = \alpha_i^n \mathbf{Z}_i \quad (10)$$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Two DOF system with viscous damping.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Absolute value of the first component of $\mathbf{U}_0$ called $|U_0^1|$; and (b) absolute value of the second component of $\mathbf{U}_0$ called $|U_0^2|$.</figcaption></figure>

This particular solution verifies Eq. (9), which becomes

$$ \forall n \ge 1 \quad ((1 + \alpha_i^2)\mathbf{A} + 2\alpha_i \bar{\mathbf{K}})\mathbf{Z}_i = 0 \quad (11) $$

To continue the developments, it is necessary to link the vector $\mathbf{Z}_i$ to a new class of modes, which will be presented in the next section.

## Eigenmodes of Perturbation

### Definition

A new eigenvalue problem is defined as follows:

$$ \forall i = 1, 2, \dots, N \quad \lambda_i \bar{\mathbf{K}} \mathbf{Y}_i = \mathbf{A} \mathbf{Y}_i \quad (12) $$

For the $m$ eigenmodes $\mathbf{Y}_i$ with $\lambda_i \ne 0$ called eigenmodes of perturbation, Eq. (12) is written as

$$ \forall i = 1, 2, \dots, m \quad \bar{\mathbf{K}} \mathbf{Y}_i = -\Theta_i \mathbf{A} \mathbf{Y}_i \quad \text{with} \quad \lambda_i = -\frac{1}{\Theta_i} \quad (13) $$

The shape modes associated with $\lambda_i = 0$ are such that $\mathbf{A}\mathbf{Y}_i = 0$, constituting the kernel of $\mathbf{A}$, which has dimension $N - m$. The eigenvectors are mass normalized

$$ \mathbf{Y}_i^T \mathbf{M} \mathbf{Y}_i = 1 \quad (14) $$

By using the symmetry of the matrix $\bar{\mathbf{K}}$ and $\mathbf{A}$, the following relationships are obtained:

$$ \mathbf{Y}_i^T \mathbf{A} \mathbf{Y}_j = a_i \delta_{ij} \quad (15) $$

$$ \mathbf{Y}_i^T \bar{\mathbf{K}} \mathbf{Y}_j = -\Theta_i a_i \delta_{ij} \quad (16) $$

$$ \mathbf{Y}_i^T \mathbf{K} \mathbf{Y}_i = \omega^2 - \Theta_i a_i = b_i \quad (17) $$

where $\mathbf{Y}_i$, $a_i$, $\Theta_i$ are a function of $\omega$ and $\delta_{ii}$ is the Kronecker function.

In the case of a damped system associated with a viscous damping matrix $\mathbf{C}$, the eigenmodes of perturbation $\mathbf{Y}_i$ and the eigenvalue $\Theta_i$ are complex and the normalizing equation Eq. (14) must be modified

$$ \tilde{\mathbf{Y}}_i^T \mathbf{M} \mathbf{Y}_i = 1 \quad (18) $$

where $\tilde{\mathbf{Y}}_i$ is the complex conjugate of $\mathbf{Y}_i$.

## Steady State Response Expression Using the Eigenmodes of Perturbation

The solution of the steady state dynamic system defined in Eq. (2) for a given value of the disturbance parameter $x$ can be expressed in the base of the eigenmodes of perturbation

$$ \mathbf{U}(x) = \mathbf{U}_A + \sum_{i=1}^{m} q_i \mathbf{Y}_i \quad \text{with} \quad m \le N \quad (19) $$

where $q_i$ is the $i$th modal coordinate of the $i$th eigenmode of perturbation $\mathbf{Y}_i$.

The term $\mathbf{U}_A$ belongs to the kernel of $\mathbf{A}$. Using this previous response expression in Eq. (2), we obtain

$$ \sum_{i=1}^{m} q_i \bar{\mathbf{K}} \mathbf{Y}_i + \sum_{i=1}^{m} q_i x \mathbf{A} \mathbf{Y}_i + \bar{\mathbf{K}} \mathbf{U}_A = \mathbf{F} \quad (20) $$

multiplying this last equation by $\bar{\mathbf{K}}^{-1}$ gives

$$ \mathbf{U}_A = \bar{\mathbf{U}} - \sum_{i=1}^{m} q_i \left(1 - \frac{x}{\Theta_i}\right) \mathbf{Y}_i \quad (21) $$

with $\bar{\mathbf{U}} = \bar{\mathbf{K}}^{-1} \mathbf{F}$ as the nondisturbance steady state solution. Then, by reporting Eq. (21) in Eq. (19)

$$ \mathbf{U}(x) = \bar{\mathbf{U}} + \sum_{i=1}^{m} q_i \frac{x}{\Theta_i} \mathbf{Y}_i \quad (22) $$

By using the orthogonality relationship in Eq. (15), the expression of $q_i$ is also given as follows:

$$ q_i = \frac{\mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}}}{a_i \left(1 - \frac{x}{\Theta_i}\right)} \quad (23) $$

Finally, the general dynamical response of the perturbed system is

$$ \mathbf{U}(x) = \left[ \mathbf{I} + \sum_{i=1}^{m} \frac{x}{a_i(\Theta_i - x)} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \right] \bar{\mathbf{U}} \quad (24) $$

or with $\lambda_i$

$$ \mathbf{U}(x) = \left[ \mathbf{I} - \sum_{i=1}^{m} \frac{x \lambda_i}{a_i (1 + x \lambda_i)} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \right] \bar{\mathbf{U}} \quad (25) $$

The last expressions of the response of the system defined in Eq. (2) can be used directly to compute statistical moments of the response. Indeed, in the context of the Monte Carlo method, they are used to obtain the solution for each $x$ without having to reverse the linear system defined in Eq. (2). But in this study, these relationships will be used to obtain a closed-form solution of this class of random dynamic systems.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Plate with localized stiffness.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Plate with an additional thickness.</figcaption></figure>

## Illustration Examples

Two examples are provided to illustrate the previous section. The plate shown in Figs. 4 and 5 is 60 cm long and 40 cm wide with a thickness of 2 mm, and has the material properties of steel: a Young's modulus, $E$, equal to 210 GPa, a Poisson's ratio $\nu$ of 0.3, and a density $\rho$ of 7,800 kg m$^{-3}$. The discrete model of this plate was obtained by the finite-element method using only 35 rectangular elements. In the case of a localized stiffness at Point P (see Fig. 4), the eigenvalue problem, Eq. (12), has only one nonzero eigenvalue $\Theta_1$ due to the form of the localized stiffness, which gives a diagonal matrix $\mathbf{A}$ with only the nonzero term $k$ at the node of the finite-element model corresponding to Point P. The nondimensional parameter $D$ is introduced

$$D = \frac{k}{k_{eq}} \qquad (26)$$

with $k_{eq}$ = stiffness of the plate on Point P approximated using the finite-element model.

The evolution of the eigenvalue $\Theta_1$ is shown in Fig. 6 for a value of $D$ equal to 0.25. In this case, $\Theta_1 = 0$ corresponds to the resonance frequencies of the nondisturbed system and the asymptotes are the antiresonance frequencies at Point P.

In this very specific case, if we limit the analysis of the displacement $u(P)$ at Point P, Eq. (25) gives

$$\frac{u(P)}{\bar{u}(P)} = \frac{1}{1 - \frac{x}{\Theta_1}} \qquad (27)$$

If we introduce the mechanical impedance $Z(\omega)$ of the plate at Point P, we can write the following equilibrium equations:

$$F = (kx + Z)u(P) \qquad F = Z\bar{u}(P) \qquad (28)$$

By eliminating the applied force at Point P between the two previous equations, we obtain the following relation:

$$\frac{u(p)}{\bar{u}(p)} = \frac{1}{1 + kxZ^{-1}} \qquad (29)$$

By comparing Eqs. (27) and (28), it appears that the unique spectral value $\Theta_1$ is given by

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of $\Theta_1$ for the localized stiffness plate with a value of $D = 0.25$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution of $\Theta_i$ for the plate with an extra thickness.</figcaption></figure>

$$\Theta_1 = -k^{-1}Z(\omega) \qquad (30)$$

In the case of a plate stiffened along one of its sides, shown in Fig. 5, the number $m$ of eigenmodes of perturbation depends on the number of the nodes of the finite-element model. The values of the $\Theta_i$ given in Fig. 7 are five in number because the discrete model of the plate gives a matrix $\mathbf{A}$ whose codimension of the kernel is five.

## Closed-Form Solution of the Random Dynamic System

### Method Exposition

Comparing Eq. (11) with the eigenvalue Eq. (13), we obtain the following relationships:

$$\forall i \quad \mathbf{Z}_i = \mathbf{Y}_i \qquad (31)$$

$$\alpha_i^2 - 2\Theta_i\alpha_i + 1 = 0 \qquad (32)$$

The evolution of $\alpha_i$ can be studied in two different cases.

When $|\Theta_i| > 1$, the second-order Eq. (32) has two real solutions, but only that with the module less than 1 must be taken into account to ensure the convergence of the power series of Eq. (10) (cf. Fig. 8)

$$\begin{cases} \forall \Theta_i > 1 \quad \alpha_i(\Theta_i) = \Theta_i - \sqrt{\Theta_i^2 - 1} \\ \forall \Theta_i < -1 \quad \alpha_i(\Theta_i) = \Theta_i + \sqrt{\Theta_i^2 - 1} \end{cases} \qquad (33)$$

When $|\Theta_i| < 1$, solutions $\alpha_1$ and $\alpha_2$ are conjugate complex roots with a unitary modulus. In order to choose between the two
<figure><img src="image_9.png" /><figcaption>Fig. 8. $|\Theta_i| > 1$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. $|\Theta_i| < 1$.</figcaption></figure>

solutions, it is necessary to introduce a damping term as a damping matrix C in Eq. (8) and choose the complex solution with the modulus less than 1. By passing to the limit by making the matrix C tend toward 0, it can be shown that the solution to be taken into account is that of the upper half circle (cf. Fig. 9)

$$|\Theta_i| < 1 \quad \alpha_i = \cos \beta_i + i \sin \beta_i \quad \text{with} \quad \cos \beta_i = \Theta_i \quad (34)$$

A physical interpretation of this choice will be given in the case of a simple degrees of freedom system by attaching it to propagation waves and the Sommerfield boundary conditions.

The general solution $\mathbf{U}_n$ is the combination of the particular solution defined in Eq. (10)

$$\forall n \ge 1 \quad \mathbf{U}_n = \sum_{i=1}^{m} \alpha_i^n H_i \mathbf{Y}_i + \overline{\mathbf{U}_n} \quad (35)$$

with the unknown $H_i$, which is the participation of each particular solution; and $\overline{\mathbf{U}_n}$ belongs to the kernel of the operator A.

According to the PCE, the response becomes

$$\mathbf{U}(x) = \mathbf{U}_0 + \sum_{n=1}^{\infty} \left( \sum_{i=1}^{m} \alpha_i^n H_i \mathbf{Y}_i + \overline{\mathbf{U}_n} \right) P_n(x) \quad (36)$$

The analytical calculation of the coefficients $H_i$ and the vectors $\overline{\mathbf{U}_n}$ using Eq. (24) is given in the Appendix, and the results are

$$H_i = (1 + \alpha_i^2) \frac{\mathbf{Y}_i \mathbf{A} \bar{\mathbf{U}}}{a_i} \quad (37)$$

$$n \ge 1 \quad \overline{\mathbf{U}_n} = \mathbf{0} \quad (38)$$

The vector $\mathbf{U}_0$ is obtained in the following form:

$$\mathbf{U}_0 = \bar{\mathbf{U}} + \sum_{i=1}^{m} \alpha_i^2 \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \quad (39)$$

If we consider that the probability density of the random variable $x$ is given in Fig. 1, the orthogonality relationships of Eq. (7) show that $\mathbf{U}_0$ corresponds to the statistic mean value of the response $\mathbb{E}(\mathbf{U})$.

The analytical variance expression of the response for each component $\mathbf{U}^k$ of the response vector $\mathbf{U}$ is equal to

$$\mathbb{V}(\mathbf{U}^k) = \mathbb{E}((\mathbf{U}^k - \mathbf{U}_0^k)(\tilde{\mathbf{U}}^k - \tilde{\mathbf{U}}_0^k)) \quad (40)$$

where $\tilde{\mathbf{S}}$ is the complex conjugate notation of term S.

By using Eq. (3), we can write

$$\mathbf{U}^k - \mathbf{U}_0^k = \sum_{n=1}^{\infty} \mathbf{U}_n^k P_n(x) \quad (41)$$

By using Eq. (35), the variance of the component $\mathbf{U}^k$ is expressed as

$$\mathbb{V}(\mathbf{U}^k) = \left( \sum_{n=1}^{\infty} \sum_{i=1}^{m} \alpha_i^n H_i Y_i^k \right) \left( \sum_{l=1}^{\infty} \sum_{j=1}^{m} \tilde{\alpha}_j^l \tilde{H}_j Y_j^k \right) \mathbb{E}(P_n P_l) \quad (42)$$

By applying the orthogonality relation from Eq. (7), we obtain

$$\mathbb{V}(\mathbf{U}^k) = \sum_{i=1}^{m} \sum_{j=1}^{m} H_i \tilde{H}_j \tilde{Y}_j^k Y_i^k \sum_{n=1}^{\infty} (\alpha_i \tilde{\alpha}_j)^n \quad (43)$$

The following identity can be used:

$$\sum_{n=1}^{\infty} (\alpha_i \alpha_j)^n = \sum_{n=1}^{\infty} \frac{\alpha_i \tilde{\alpha}_j}{1 - \alpha_i \tilde{\alpha}_j} \left(1 - \lim_{n \to +\infty} (\alpha_i \tilde{\alpha}_j)^n\right) \quad (44)$$

by taking into account the damping of the mechanical system, $\lim_{n \to +\infty} (\alpha_i \alpha_j)^n = 0$ because $|\alpha_i \alpha_j| < 1$.

Finally, the variance expression becomes

$$\mathbb{V}(\mathbf{U}^k) = \sum_{i=1}^{m} \sum_{j=1}^{m} \frac{\alpha_i \tilde{\alpha}_j}{1 - \alpha_i \tilde{\alpha}_j} (1 + \alpha_i^2)(1 + \tilde{\alpha}_j^2) \tilde{Y}_j^k Y_i^k d_i \tilde{d}_j \quad (45)$$

with

$$d_i = \frac{\mathbf{Y}_i \mathbf{A} \bar{\mathbf{U}}}{a_i} \quad (46)$$

## Interpretation of Results Using Simple Models

In this section, the analytical results determined previously will be applied to simple systems to physically interpret the polynomial chaos expansion proposed. The first random system studied is a single DOF (1DOF) system without damping subjected to the harmonic force shown in Fig. 10.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Single DOF system.</figcaption></figure>

The governing equation of the system is

$$(-\omega^2 m + k + xA)U = F \\ \text{with } u = U \exp(i\omega t) \quad \text{and} \quad f = F \exp(i\omega t) \quad (47)$$

In this simple case, the unique solution $\Theta_1$ is given by the solution of the spectral problem defined in Eq. (13)

$$A\Theta_1 = -k + \omega^2 m \quad (48)$$

Nondimensional terms are introduced: $D = (A/k) < 1$ is the dispersion factor $U_{stat} = (F/k)$, and $\Omega = (\omega/\omega_0)$ is the reduced frequency where $\omega_o = \sqrt{(k/m)}$ is the natural frequency. Eq. (48) then becomes

$$D\Theta_1 = -1 + \Omega^2 \quad (49)$$

The curve giving $D\Theta_1$ as a function of $\Omega$ in the event of $D = 0.2$ highlights a frequency band where the $\alpha_i$ value is imaginary with a module of 1 while its limit values are $\Omega_1 = \sqrt{1 - \bar{D}}$ and $\Omega_2 = \sqrt{1 + \bar{D}}$ (cf. Fig. 11). A physical interpretation can link this zone to a propagation band of a discrete semi-infinite periodic system. In this way, the equation system in Eq. (8) in the case of a single DOF system can be illustrated by the mechanical model shown in Fig. 12. In order to obtain a physical model, coefficient $A$ has been replaced by $-h$, the latter parameter $h$ considered positive is interpreted as the stiffness of a spring. The frequency band [$\Omega_1; \Omega_2$] corresponds to a propagation zone, and the propagation coefficient $\alpha_1$ must have a negative imaginary part in order to comply with the Sommerfeld condition. In Eq. (32), the solution chosen $\alpha_i$ must have a positive imaginary part because coefficient $D$ is considered positive, while in the mechanical system, $h$ is positive and $A$ is negative.

According to Eq. (39), the mean response $U_0$ is obtained as follows:

<figure><img src="image_12.png" /><figcaption>Fig. 12. Semi-infinite periodic system.</figcaption></figure>

$$U_0 = \bar{U} + \frac{1}{a_1} \alpha_1^2 Y_1 Y_1^T A \bar{U} \quad (50)$$

Note that $Y_1 A Y_1^T$ is equal to $a_1$, and the terms $Y_1 Y_1^T A$ can be switched in this case of 1DOF system; therefore the previous equation becomes

$$U_0 = \bar{U}(1 + \alpha_1^2) \quad \text{with} \quad \bar{U} = \frac{F}{k(1 - \Omega^2)} \quad (51)$$

Three zones have to be considered (cf. Fig. 11)

* Zone 1: $\Omega^2 < \Omega_1^2$ so $\Theta_1 < -1$

$$\alpha_1 = \Theta_1 + \sqrt{\Theta_1^2 - 1} \quad (52)$$

According to Eq. (51)

$$U_0 = \frac{2F}{kD} \frac{1}{(-\Theta_1 + \sqrt{\Theta_1^2 - 1})} \quad (53)$$

Using Eq. (49)

$$U_0 = \frac{2F}{k(1 - \Omega^2 + \sqrt{(1 - \Omega^2)^2 - D^2})} \quad (54)$$

In the static case, $\Omega^2 = 0$, therefore Eq. (54) becomes

$$U_0 = \bar{U} \frac{2}{1 + \sqrt{1 - D^2}} \quad \text{with} \quad \bar{U} = \frac{F}{k} = U_{stat} \quad (55)$$

Using equation Eq. (7), it can be demonstrated that this result is exactly the same as the statistical mean of the static response of a random spring

<figure><img src="image_13.png" /><figcaption>Fig. 11. Characteristic curve with its propagation band.</figcaption></figure>
$$U_0 = \frac{2F}{k\pi} \int_{-1}^{1} \frac{\sqrt{1-u^2}}{1+Du} du \quad (56)$$

* Zone 2: $\Omega^2 > \Omega_2^2$ so $\Theta_1 > 1$

$$\alpha_1 = \Theta_1 - \sqrt{\Theta_1^2 - 1} \quad (57)$$

The same calculation as that used previously gives the following result:

$$U_0 = \frac{2F}{k(1 - \Omega^2 - \sqrt{(1 - \Omega^2)^2 - D^2})} \quad (58)$$

* Zone 3: $\Omega_1^2 < \Omega^2 < \Omega_2^2$ so $\Theta_1 \in [-1; 1]$

$$\alpha_1 = \Theta_1 + i\sqrt{1 - \Theta_1^2} \quad (59)$$

According to Eqs. (49) and (51)

$$U_0 = \frac{2F}{k(1 - \Omega^2 + i\sqrt{D^2 - (1 - \Omega^2)^2})} \quad (60)$$

It is now possible to evaluate the amplitude at the resonance frequency when $\Omega^2 = 1$ (so $\Theta_1 = 0$)

$$U_0 = \frac{2U_{stat}}{iD} \quad (61)$$

This last value can be compared to that viscous damped single degree of freedom system

$$U_0(\Omega^2 = 1) = \frac{U_{stat}}{i2\xi} \quad (62)$$

where $\xi$ is the specific damping factor. This analysis gives an equivalence between the structural damping $\eta = 2\xi$ and $D/2$ the dispersion damping factor. This equivalent damping associated with the uncertainty can be considered as a numerical radiation loss in the fictitious semi-infinite periodic mechanical model shown in Fig. 12. As can be seen in Fig. 13, the mean response modulus $|U_0|$ is constant in the third zone of the frequency domain. The real and imaginary complex parts of the mean response are given in Figs. 14 and 15. It can be seen that the imaginary parts are equal to zero in the first and second zones. In the case where we introduce a damping factor $c$ in the 1DOF system, the eigenvalue $\Theta_1$ is given by

$$A\Theta_1 = -k + \omega^2 m + i\omega c \quad (63)$$

<figure><img src="image_14.png" /><figcaption>Fig. 13. $|U_0|$ with respect to $\Omega$ for different values of $D$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Real part of $U_0$ with respect to $\Omega$ for different values of $D$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Imaginary part of $U_0$ with respect to $\Omega$ for different values of $D$.</figcaption></figure>

or in nondimensional form

$$D\Theta_1 = -1 + \Omega^2 + 2i\xi\Omega \quad (64)$$

Solution $\alpha_1$ is obtained with the complex second-order equation, Eq. (32). The mean response modulus $|U_0|$ has been plotted for different values of $\xi$, as shown in Fig. 16. The expression of the response variance, given in Eq. (45) becomes in this simple case

$$\mathcal{V}(U) = |U_0|^2 \frac{|\alpha_1|^2}{1 - |\alpha_1|^2} \quad (65)$$

<figure><img src="image_17.png" /><figcaption>Fig. 16. $|U_0|/U_{stat}$ with respect to $\Omega$ for different values of $\xi$ for $D = 0.01$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Normalized standard deviation σ with respect to Ω for different values of D for ξ = 0.02.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Two DOF system.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. First eigenvector.</figcaption></figure>

Thus, the normalized standard deviation σ is defined by

$$ \sigma = \frac{\sqrt{\mathcal{V}(U)}}{|U_0|} = \frac{|\alpha_1|}{\sqrt{1 - |\alpha_1|^2}} \quad (66) $$

The evolution of σ with respect to frequency Ω for many values of D is shown in Fig. 17. In the case of a 2DOF system, schema- tized in Fig. 18, the equilibrium equations can be written in the following matrix form:

$$ k \begin{bmatrix} -\Omega^2 + 2 + D\Theta & -(1 + D\Theta) \\ -(1 + D\Theta) & -\Omega^2 + 1 + D\Theta \end{bmatrix} \begin{bmatrix} U^1 \\ U^2 \end{bmatrix} = \begin{bmatrix} F \\ 0 \end{bmatrix} \quad (67) $$

with $\omega_0^2 = (k/m)$, $\Omega^2 = (\omega^2/\omega_0^2)$ and $D = (A/k)$.

The nullity of the determinant of this matrix gives the character- istic equation

$$ \Omega^4 + (-2D\Theta - 3)\Omega^2 + D\Theta + 1 = 0 \quad (68) $$

Eq. (68) gives only one value for Θ

$$ D\Theta_1 = \frac{\Omega^4 - 3\Omega^2 + 1}{2\Omega^2 - 1} = -\frac{1}{\lambda_1}D \quad (69) $$

Thus, there are two solutions for the eigenvalue problem given in Eq. (12), with only one different from zero

$$ \lambda_1 = \frac{-1}{\Theta_1}, \quad \lambda_2 = 0 \quad (70) $$

The first eigenvector is $Y_1$ corresponding to the unique eigenm- odes of perturbation

$$ Y_1 = \begin{pmatrix} \Omega^2 \\ 1 - \Omega^2 \end{pmatrix} r_1^{-1} \quad (71) $$

where $r_1$ is calculated with respect to the normalization Eq. (14): $r_1 = m^{1/2}(\Omega^4 + (1 - \Omega^2)^2)^{1/2}$.

The mode shape of the first eigenmodes of perturbation is shown in Fig. 19. The curve in Fig. 20 corresponding to D = 0.2 shows two propagation bands associated with Θ₁ and that Ω = $\sqrt{1/2}$ corresponds to Θ₁ → ∞. The αᵢ are calculated in each zone defined in Fig. 20. For example, in the two propagation bands

<figure><img src="image_21.png" /><figcaption>Fig. 20. Characteristic curve.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 21. (a) Absolute value of the first component of $\mathbf{U}_0$ called $|\mathbf{U}_0^1|$ with different values of $D$ for $\xi = 0$; and (b) absolute value of the second component of $\mathbf{U}_0$ called $|\mathbf{U}_0^2|$ with different values of $D$ for $\xi = 0$.</figcaption></figure>

$$ \alpha_1 = \frac{\Omega^4 - 3\Omega^2 + 1}{D(2\Omega^2 - 1)} + i\sqrt{1 - \frac{(\Omega^4 - 3\Omega^2 + 1)^2}{D^2(2\Omega^2 - 1)^2}} \quad (72) $$

The mean response is calculated by applying Eq. (50) at the two degrees of freedom system

$$ \mathbf{U}_0 = \left( \mathbf{I} + \frac{\alpha_1^2}{2\Omega^2 - 1} \begin{bmatrix} \Omega^2 & -\Omega^2 \\ 1 - \Omega^2 & -1 + \Omega^2 \end{bmatrix} \right) \bar{\mathbf{U}} \quad (73) $$

$$ \text{with } \bar{\mathbf{U}} = \bar{\mathbf{K}}^{-1} \begin{bmatrix} 0 \\ F \end{bmatrix} \quad (74) $$

The mean value modulus is given in Fig. 21 for several values of $D$. If viscous damping is introduced, the eigenvector $\mathbf{Y}_1$ and $\lambda_1$ are complex. In this case of the 2DOF system defined in Fig. 2, the mean value modulus is given in Fig. 22 for several values of $D$ and the damping coefficient $\xi$ equal to $2.10^{-2}$.

## Simplified Analysis in the Case of Complex Structures

### Use of the Modal Base

In the case of complex structures, the eigenmodes of perturbation can be approximated using the normal modes $\mathbf{X}_i$ defined by

$$ \mathbf{K}\mathbf{X}_i = \omega_i^2\mathbf{M}\mathbf{X}_i \quad \text{with} \quad \mathbf{X}_i^T\mathbf{M}\mathbf{X}_j = \delta_{ij} \quad (75) $$

The solution $\mathbf{U}$ is approximated as a function of the $n$ first modes associated with the frequency band of the analysis

<figure><img src="image_23.png" /><figcaption>Fig. 22. (a) Absolute value of the first component of $\mathbf{U}_0$ called $|\mathbf{U}_0^1|$ with different values of $D$ for $\xi = 0.02$; and (b) absolute value of the second component of $\mathbf{U}_0$ called $|\mathbf{U}_0^2|$ with different values of $D$ for $\xi = 0.02$.</figcaption></figure>

$$ \mathbf{U}(x) = \sum_{i=1}^{n} \mathbf{X}_i q_i = \mathbf{X} \mathbf{q} \quad (76) $$

The dynamic system equation becomes

$$ (\mathbf{\Lambda}^2 - \omega^2 \mathbf{I}) \mathbf{q} + \mathbf{B} \dot{\mathbf{q}} + x \mathbf{X}^T \mathbf{A} \mathbf{X} = \mathbf{f} \quad (77) $$

where $\mathbf{f} = \mathbf{T}\mathbf{X}\mathbf{F}$; and $\mathbf{\Lambda}^2 = \begin{bmatrix} \omega_1^2 & 0 \\ 0 & \omega_i^2 \end{bmatrix}$ is the spectral matrix.

In this case, $\mathbf{B}$ is the modal damping matrix associated with the damping matrix $\mathbf{C}$, which models the dissipation of this system

$$ \mathbf{B} = \mathbf{X}^T \mathbf{C} \mathbf{X} \quad (78) $$

The eigenmodes of perturbation are the solutions of the following truncated dynamic system:

$$ (\mathbf{\Lambda}^2 + i\omega\mathbf{B} - \omega^2\mathbf{I})\mathbf{Q}_i = -\Theta_i \bar{\mathbf{A}}\mathbf{Q}_i \quad \text{with} \quad \bar{\mathbf{A}} = \mathbf{X}^T\mathbf{A}\mathbf{X} \quad (79) $$

In this case, $\mathbf{Q}_i$ and $\Theta_i$ are frequency dependent. The approximated eigenmodes of perturbation are given by

$$ \mathbf{Y}_i = \mathbf{X}\mathbf{Q}_i \quad (80) $$

The computation time of the eigenmodes of perturbation with respect to $\omega$ is considerably reduced because of the modal truncation. If we assume that the matrixes $\mathbf{A}$ and $\mathbf{B}$ are diagonalizable in the modal base, $\mathbf{Y}_i$ and $\mathbf{X}_i$ are identical. As in the case of the Basile hypothesis for the damping matrix, there is no reason to make this hypothesis for the disturbance matrix $\mathbf{A}$. However, this simplification is justified if the matrix $\mathbf{A}$ is small compared to $\mathbf{K}$. Under this condition, we have
$$U_0 = \left( \mathbf{I} + \sum_{i=1}^{n} \alpha_i^2 \frac{1}{a_i} \mathbf{X}_i \mathbf{X}_i^T \bar{\mathbf{A}} \right) \bar{\mathbf{U}} \quad (81)$$

The $\Theta_i$ is easily calculated with the following relationships:

$$\omega_i^2 - \omega^2 + ib_i\omega = -\Theta_i \bar{A}_{ii} \quad (82)$$

where $\bar{A}_{ii} = \mathbf{X}_i^T \mathbf{A} \mathbf{X}_i$; $b_i = \mathbf{X}_i^T \mathbf{C} \mathbf{X}_i$; and $\alpha_i$ is the solution with its modulus <1 from the following second-order equation:

$$\alpha_i^2 - 2\Theta_i\alpha_i + 1 = 0 \quad (83)$$

Given that matrix **B** is considered diagonal in agreement with the Basile hypothesis, the solution $\bar{\mathbf{U}}$ can be easily approximated by

$$\bar{\mathbf{U}} = \left( \sum_{j=1}^{n} \frac{\mathbf{X}_j^T \mathbf{X}_j}{-\omega^2 + \omega_j^2 + i\omega b_j} \right) \mathbf{F} \quad (84)$$

By using this relationship in Eq. (81), we obtain

$$\mathbf{U}_0 = \sum_{j=1}^{n} \frac{\mathbf{X}_j \mathbf{X}_j^T \mathbf{F}}{-\omega^2 + \omega_j^2 + i\omega b_j} + \sum_{j=1}^{n} \sum_{i=1}^{n} \frac{\alpha_i^2}{a_i(-\omega^2 + \omega_j^2 + i\omega b_j)} \mathbf{X}_i \mathbf{X}_i^T \bar{\mathbf{A}} \mathbf{X}_j \mathbf{X}_j^T \mathbf{F} \quad (85)$$

Thus, the mean value $\mathbf{U}_0$ corresponds to the sum of the mean value of each mode of the damped system

$$\mathbf{U}_0 = \sum_{j=1}^{n} \frac{1 + \alpha_j^2}{-\omega^2 + \omega_j^2 + i\omega b_j} \mathbf{X}_j \mathbf{X}_j^T \mathbf{F} \quad (86)$$

Of course, this very simple expression of mean response is usable only when the stiffness disturbance and the damping of the dynamic system are weak.

## Illustration Example

In the case of the 2DOF system shown in Fig. 2, this simplification was tested with $D = 0.1$ and $\xi = 0.02$, and it can be seen in Fig. 23 that the mean modulus response approximated by Eq. (86) is close to the exact solution, although matrices **A** and **B** are not diagonalizable in the modal base of the 2DOF system.

## Conclusion

The concept of perturbation eigenmodes defined in this study made it possible to express the dynamic steady-state of a system as a function of a structural parameter. Thanks to this direct relationship, it is possible to apply the Monte Carlo method at low cost in the case of a dynamic system having a random parameter. Considering that the probability density of the parameter is associated with the Chebyshev polynomials of the second kind, the notion of eigenmodes of perturbation led to a closed form of the random steady-state solution. This result was obtained by linking the PCE method to a semi-infinite periodic system. In particular, this new method allowed interpreting the damping effect observed on the system mean response like a loss of energy associated with the wave propagation in a fictitious mechanical system subject to the Sommerfeld boundary condition. More precisely, it appears that each eigenmode of perturbation is associated with a propagation mode in a fictitious semiperiodic system. In order to simplify

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

<figure><img src="image_25.png" /><figcaption>Fig. 23. (a) Absolute value of the first component of $\mathbf{U}_0$ called $|U_0^1|$; and (b) absolute value of the second component of $\mathbf{U}_0$ called $|U_0^2|$.</figcaption></figure>

the calculation of the eigenmodes of perturbation in the case of complex structures, an approximate method based on the use of a truncated modal base has been proposed. The study was limited to the case of a single random parameter with the main objective of presenting an original vision of a class of random dynamic systems and its extension to the case of several parameters was linked to the response of a periodic medium with several dimensions, which requires the development of very specific methods in the field of acoustics.

## Appendix. Determination of un et $H_i$ Defined in the Closed Form Solution

To calculate coefficient $H_i$ and the vector $\bar{\mathbf{U}}_n$, we use Eq. (24) and replace the $\Theta_i$ by its expression with the $\alpha_i$ given in Eq. (32)

$$\mathbf{U}(x) = \bar{\mathbf{U}}P_0(x) + \sum_{i=1}^{m} \frac{2\alpha_i x}{1 + \alpha_i^2 - 2\alpha_i x} \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \quad (87)$$

The ordinary generating function of the Chebyshev polynomials of the second kind allows us to write (Mason 1993)

$$\frac{2\alpha_i x}{1 + \alpha_i^2 - 2\alpha_i x} = \sum_{n=0}^{\infty} 2\alpha_i^{n+1} x P_n(x) \quad (88)$$

The recurrence relation Eq. (4) applied to the previous equation gives

$$\sum_{n=1}^{\infty} 2\alpha_i^{n+1} x P_n(x) = \sum_{n=1}^{\infty} \alpha_i^{n+1} (P_{n+1}(x) + P_{n-1}(x)) \quad (89)$$
Eq. (87) becomes

$$ 
\begin{aligned}
\mathbf{U}(x) = \bar{\mathbf{U}}P_0(x) &+ \sum_{n=1}^{\infty} \sum_{i=1}^{m} \alpha_i^{n+1} (P_{n+1}(x) + P_{n-1}(x)) \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \\
& + \sum_{i=1}^{m} 2\alpha_i x P_0(x) \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}}
\end{aligned}
\quad (90) $$

with the identities

$$ P_0(x) = 1 \quad \text{and} \quad P_1(x) = 2x \quad (91) $$

By changing the subscript of $P_{n+1}(x)$ and $P_{n-1}(x)$, we obtain

$$ 
\begin{aligned}
\mathbf{U}(x) = \bar{\mathbf{U}} &+ \sum_{n=0}^{\infty} \sum_{i=1}^{m} \alpha_i^{n+2} (P_{n+2}(x) + P_n(x)) \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \\
& + \sum_{i=1}^{m} \alpha_i P_1(x) \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}}
\end{aligned}
\quad (92) $$

By calculating the coefficients of the polynomials $P_n(x)$ with respect to Eq. (3), we obtain

$$ \forall n > 0 \quad \mathbf{U}_n = \sum_{i=1}^{m} \alpha_i^n (1 + \alpha_i^2) \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \quad (93) $$

By identification with Eq. (36), the following relationship can be written:

$$ H_i = (1 + \alpha_i^2) \frac{\mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}}}{a_i} \quad (94) $$

Otherwise, the terms $\bar{\mathbf{U}}_n$ for $n \ge 1$ are equal to zero thanks to Eq. (92).

The first term $\mathbf{U}_0$ is given by

$$ n = 0 \quad \mathbf{U}_0 = \bar{\mathbf{U}} + \sum_{i=1}^{m} \alpha_i^2 \frac{1}{a_i} \mathbf{Y}_i \mathbf{Y}_i^T \mathbf{A} \bar{\mathbf{U}} \quad (95) $$

## Data Availability Statement

All data, models, and code generated or used during the study appear in the published article.

## References

Argoul, P., and L. Jezequel. 1989. "Improvement of a nonparametric identification procedure used in nonlinear dynamics." *J. Appl. Mech.* 56 (3): 697-703. https://doi.org/10.1115/1.3176149.

Baker, G. A., and P. Graves-Morris. 1996. *Padé Approximants*. 2nd ed. Cambridge, UK: Cambridge University Press.

Brezinski, C. 1996. "Extrapolation algorithms and Padé approximations: A historical survey." *Appl. Numer. Math.* 20 (3): 299-318.

Ghanem, R., and P. D. Spanos. 1991. *Stochastic finite elements: A spectral approach*. New York: Springer.

Huang, S. P., S. T. Quek, and K. K. Phoon. 2001. "Convergence study of the truncated Karhunen-Loève expansion for simulation of stochastic processes." *Int. J. Numer. Methods Eng.* 52 (9): 1029-1043. https://doi.org/10.1002/nme.255.

Jacquelin, E., S. Adhikari, J.-J. Sinou, and M. I. Friswell. 2015a. "Polynomial chaos expansion and steady-state response of a class of random dynamical systems." *J. Eng. Mech.* 141 (4): 04014145. https://doi.org/10.1061/(ASCE)EM.1943-7889.0000856.

Jacquelin, E., S. Adhikari, J.-J. Sinou, and M. I. Friswell. 2015b. "Polynomial chaos expansion in structural dynamics: Accelerating the convergence of the first two statistical moment sequences." *J. Sound Vib.* 356 (Nov): 144-154. https://doi.org/10.1016/j.jsv.2015.06.039.

Jacquelin, E., O. Dessombz, J.-J. Sinou, S. Adhikari, and M. I. Friswell. 2017. "Polynomial chaos-based extended Padé expansion in structural dynamics." *Int. J. Numer. Methods Eng.* 111 (12): 1170-1191. https://doi.org/10.1002/nme.5497.

Keshavarzzadeh, V., R. G. Ghanem, S. F. Masri, and O. J. Aldraihem. 2014. "Convergence acceleration of polynomial chaos solutions via sequence transformation." *Comput. Methods Appl. Mech. Eng.* 271 (Apr): 167-184. https://doi.org/10.1016/j.cma.2013.12.003.

Kleiber, M., and T. D. Hien. 1992. *The stochastic finite element method: Basic perturbation technique and computer implementation*. Hoboken, NJ: Wiley.

Li, X.-B., F.-Q. Gong, and J. Deng. 2006. "Probabilistic density function estimation of geotechnical shear strength parameters using the second Chebyshev orthogonal polynomial." *J. Cent. South Univ.* 13 (3): 275-280. https://doi.org/10.1007/s11771-006-0123-4.

Liu, J., X. Li, and L. Wu. 2016. "An operational matrix technique for solving variable order fractional differential-integral equation based on the second kind of Chebyshev polynomials." *Adv. Math. Phys.* 2016: 9. https://doi.org/10.1155/2016/6345978.

Liu, X., L. Hong, and L. Yang. 2014. "Hopf bifurcations of a stochastic fractional-order van der Pol system." In *Abstract and applied analysis* 2014, 10. London: Hindawi.

Ma, S.-J., W. Xu, W. Li, and T. Fang. 2006. "Analysis of stochastic bifurcation and chaos in stochastic Duffing-van der Pol system via Chebyshev polynomial approximation." *Chin. Phys.* 15 (6): 1231-1238. https://doi.org/10.1088/1009-1963/15/6/017.

Mason, J. C. 1993. "Chebyshev polynomials of the second, third and fourth kinds in approximation, indefinite integration, and integral transforms." *J. Comput. Appl. Math.* 49 (1-3): 169-178. https://doi.org/10.1016/0377-0427(93)90148-5.

Matos, A. C. 1996. "Some convergence results for the generalized Padé-type approximants." *Numer. Algorithms* 11 (1): 255-269. https://doi.org/10.1007/BF02142501.

Wiener, N. 1938. "The homogeneous chaos." *Am. J. Math.* 60 (4): 897-936. https://doi.org/10.2307/2371268.

Xiu, D., and G. E. M. Karniadakis. 2002. "The Wiener-Askey polynomial chaos for stochastic differential equations." *SIAM J. Sci. Comput.* 24 (2): 619-644. https://doi.org/10.1137/S1064827501387826.

Yamazaki, F., M. Shinozuka, and G. Dasgupta. 1988. "Neumann expansion for stochastic finite element analysis." *J. Eng. Mech.* 114 (8): 1335-1354. https://doi.org/10.1061/(ASCE)0733-9399(1988)114:8(1335).