Published: 30 September 2013

# A fast algorithm for calculation of random response of elastic thin plate under distributed random load

Jiang Jinhui1
Chen Guoping2
Zhang Fang3
1, 2, 3Institute of Vibration Engineering Research, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
1, 2, 3State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing 210016, China
Corresponding Author:
Jiang Jinhui
Views 33

#### Abstract

The calculation of response power spectrum of elastic thin plates, the common load-carrying structures in engineering, under distributed random load has been suffering from many problems, such as low efficiency and complicated computation which cause a lot of inconvenience for engineering applications. Therefore, a fast algorithm is introduced in this article to address above problems. Firstly, the function of the cross power spectrum function is projected onto orthogonal domain of multi-dimension Legendre polynomial, and the infinite distributing information is expressed by finite parameters. Secondly, the dynamic response computation of distributed random load is transformed into that of distributed harmonic load. Thirdly, the response power spectrum density for each frequency point can be obtained by vector multiplication. The fast algorithm is accurate; and its efficiency is much higher than a regular algorithm. Moreover, the fast algorithm conveniently illustrates the coherences of distributed random load. Finally, the fast algorithm is compared with a regular algorithm computation. The high efficiency and precision of the fast algorithm is demonstrated in this paper.

## 1. Introduction

Thin plate is a very common structure in engineering. Aircrafts and satellites, for instance, consist of many thin plates or stressed structures. In these circumstances, thin plates are mostly subjected to distributed dynamic load [1], especially distributed random load. Though it is relatively simple to calculate the direct problem under distributed determinate load, the calculation of structural dynamic response under distributed random load is quite complicated and inefficient [2, 3]. The calculation of structural dynamic response under distributed random load is usually simplified in engineering, which consequently magnifies the error [4-8]. In addition, the calculation of dynamic random response focuses on the structure due to multi-points excitation [9-13]. Few studies involve the calculation of structure random response due to the distributed random load. This paper proposes a new fast algorithm for the calculation of response power spectrum density of thin plates under distributed random load with the prerequisite of linear elasticity. The fast algorithm is an accurate and highly efficient method; and the above-mentioned difficulties in calculation of structural dynamic response under distributed random load can be satisfactorily resolved by means of the fast algorithm, so it may be useful for engineering problem involving response analysis of distributed random load.

There are mainly three methods to calculate structural response power spectrum density [14]: CQC (Complete Quadratic Combination) [15-18], SRSS (Square Root of the Sum of Squares) [19-21], and PEM (Pseudo-Excitation Method) [22, 23] proposed by Lin Jiahao. CQC is an accurate algorithm since it not only takes the square items of the main mode shapes, but also considers the coupling items. However, its efficiency is terrible. For elastic thin plates, the formula of CQC method is displayed as Eq. (1):

1
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{mn}\left(-\omega \right){H}_{ij}\left(\omega \right){Q}_{mnij}\left(\omega \right),$

where ${Q}_{mnij}\left(\omega \right)={\int }_{0}^{a}{\int }_{0}^{b}{\int }_{0}^{a}{\int }_{0}^{b}{S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right)d{x}_{1}d{y}_{1}d{x}_{2}d{y}_{2}\text{,}$${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ and ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ represent cross-spectrum density of response and cross-spectrum density of excitation between points $\left({x}_{1},{y}_{1}\right)$ and $\left({x}_{2},{y}_{2}\right)\text{,}$ respectively. ${\phi }_{mn}\left({x}_{1},{y}_{1}\right)$ is mode shape of mode $mn$ at point $\left({x}_{1},{y}_{1}\right)$, ${H}_{ij}\left(\omega \right)$ is modal response frequency function of mode $ij$, and $\left[0,a\right]$, $\left[0,b\right]$ represent the effect range of variants $x$ and $y$, respectively. All the following symbols are treated same.

In order to reduce the calculation of CQC method, a simplified method which ignores all the coupling items in Eq. (1) is applied in engineering. It is called SRSS method. For elastic thin plates, the formula of SRSS method is shown as Eq. (2):

2
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{mn}\left({x}_{2},{y}_{2}\right){\left|{H}_{mn}\left(\omega \right)\right|}^{2}$
$\bullet {\int }_{0}^{a}{\int }_{0}^{b}{\int }_{0}^{a}{\int }_{0}^{b}{S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{mn}\left({x}_{2},{y}_{2}\right)d{x}_{1}d{y}_{1}d{x}_{2}d{y}_{2}.$

However, SRSS, when applied to structures with close vibration modes or large damping, suffers from large errors due to its assumption of independence among responses of all vibration modes.

PEM has been developed for mulit-input-multi-output stationary/non-stationary random vibration response problem since 1980s, thus it is a CQC algorithm with high efficiency and simple procedure [9-13]. Essentially, it approximates the random input in basic dynamic equation with a virtual determinate harmonic input and conducts vector multiplication for every single frequency to obtain response power spectrum density. Maily, this method is discussed for multiple inputs. A fast algorithm of PEM for one-dimensional continuous beam was developed [8], but the calculation of response power spectrum density of thin plates under distributed random load has not yet been involved.

The fast algorithm is proposed for accurate calculation of response power spectrum density of structures under line-distributed (one-dimensional) or plane-distributed (two-dimensional) random load on the thin plate. First, one-dimensional Legendre orthogonal polynomials are introduced to establish multi-dimensional ones. Second, the function of cross-power spectrum density, which is complicated to distribute on thin plate structures, is decomposed in orthogonal domain of multi-dimensional Legendre polynomials. Third, based on CQC method, the calculation of response power spectrum density of distributed random load is replaced with that of distributed determinate load; and by conducting vector multiplication for every single frequency, the response power spectrum density is thus gained. Regular approximated and the fast algorithms are programmed and evaluated in efficiency and precision to validate the advantages of the fast algorithm. The result of Nastran (the software for finite element analysis) demonstrates the correctness of the fast algorithm. Besides, the fast algorithm overcomes Nastran’s shortage that it can only compute the load model whose distributed load is completely coherent. In fact, not all distributed load is completely coherent, for example, fluctuating winds; if the stationary distributed random excitation is partial coherent, Nastran does not work. So, it is necessary to develop the fast algorithm which can calculate the structure random response under arbitrarily coherent distributed random load.

## 2. Theory

### 2.1. Multi-dimensional Legendre orthogonal polynomials

Located within a real span $\left[-1,1\right]$ and with weighting function $\rho \left(x\right)=1\text{,}$ Legendre orthogonal polynomials [24, 25] can be described as follows:

3
${L}_{0}\left(x\right)=1,{L}_{n}\left(x\right)=\frac{1}{{2}^{n}n!}\frac{{d}^{n}}{d{x}^{n}}\left\{{\left({x}^{2}-1\right)}^{n}\right\},\left(n=1,2,...\right).$

In $\left[0,s\right]$, recurrence formula of Legendre polynomials can be derived as below:

$\left\{\begin{array}{l}{L}_{0}\left(x\right)=1,\\ {L}_{1}\left(x\right)=\frac{2x}{s}-1,\\ \dots \\ {L}_{n+1}\left(x\right)=\frac{2n+1}{n+1}\left(\frac{2x}{s}-1\right){L}_{n}\left(x\right)-\frac{n}{n+1}{L}_{n-1}\left(x\right),\end{array}\right\\mathrm{}\mathrm{}x\in \left[0\text{\hspace{0.17em}},\text{\hspace{0.17em}}s\right],n=1,2,\dots \mathrm{}.$

Based on the above one-dimensional Legendre orthogonal polynomials, two-dimensional orthogonal polynomials can be established as ${L}_{mn}\left(x,y\right)={L}_{m}\left(x\right){L}_{n}\left(y\right)$. According to the orthogonality of one-dimensional Legendre orthogonal polynomials, the orthogonality of two-dimensional ones is easily proved as follows:

4
${\int }_{0}^{{l}_{1}}{\int }_{0}^{{l}_{2}}{L}_{mn}\left(x,y\right){L}_{ij}\left(x,y\right)dxdy=\left\{\begin{array}{ll}\frac{{l}_{1}{l}_{2}}{{\left(2n+1\right)}^{2}},& m=i&n=j,\\ 0,& \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}.\end{array}\right\$

It can also be proved that the established two-dimensional Legendre orthogonal polynomials possess the characteristics of parity and recurrence. As a result, $\left\{{L}_{mn}\left(x,y\right)\right\}$ constructs a series of two-dimensional orthogonal basis functions. Furthermore, multi-dimensional orthogonal polynomials can be constructed similarly as below:

5
$\left\{{L}_{{j}_{1}{j}_{2}\cdots {j}_{n}}\left({x}_{1},{x}_{2},\cdots ,{x}_{n}\right)\right\}={L}_{{j}_{1}}\left({x}_{1}\right){L}_{{j}_{2}}\left({x}_{2}\right)\cdots {L}_{{j}_{n}}\left({x}_{n}\right).$

In Eq. (5), $n$ represents the number of dimensions of orthogonal basis functions. The orthogonality, parity and recurrence of multidimensional orthogonal basis functions are easy to prove. In this paper, two-dimensional and four-dimensional orthogonal polynomials are involved.

If a continuous function $f\left(x,y\right)$ is square integrable ($f\left(x,y\right)\in {L}_{2}&f\left(x,y\right)\in C$, where ${L}_{2}$, $C$ is square integrable function space and continuous function space, respectively), it can be expressed via the expansion with two-dimensional Legendre orthogonal polynomials in the following form:

6
$f\left(x,y\right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{L}_{mn}\left(x,y\right){a}_{mn},$

where ${a}_{mn}$ is a expansion coefficient and can be determined easily through the application of an inner (dot) product. The partial sums of expansion series can be expressed as the best approximation of the function, known as best weighted norm approximation:

7
$g\left(x,y\right)=\sum _{m=1}^{M}\sum _{n=1}^{N}{L}_{mn}\left(x,y\right){a}_{mn}.$

Let $\epsilon$ be an arbitrary positive number, $M$, $N$ can be found to satisfy the following inequality:

8
${‖f\left(x,y\right)-g\left(x,y\right)‖}^{2}<\epsilon .$

$M$, $N$ can be easily determined via a routine written by LabVIEW.

### 2.2. General model for cross power spectrum density of distributed random load of elastic thin plates

As shown in Figure 1, the cross power spectrum density of distributed random load, which is stationary Gaussian random process in the time domain, can be described as Eq. (9):

9
${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)\sqrt{{S}_{f}\left({x}_{1},{y}_{1},\omega \right){S}_{f}\left({x}_{2},{y}_{2},\omega \right)}.$

Here, ${S}_{f}\left(x,y,\omega \right)$ is stationary auto-spectral density of the distributed random load at point $\left(x,y\right)$, which is function of space and frequency variable. For example, Davenport spectrum for fluctuating winds velocity:

10
$\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\left|\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)\right|{e}^{i\theta \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)}.$

$\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ is a coherence function and satisfies Eq. (11) below:

11
$|\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)|\le 1.$

The argument $\left(\theta \right)$ shows as:

12
$\theta \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\text{arctan}\left[\frac{\mathrm{I}\mathrm{m}\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)}{\mathrm{R}\mathrm{e}\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)}\right].$

If the auto-power spectrum density and the coherence function of every point are given, the cross-power spectrum density function of distributed random load can be uniquely determined as follows.

Fig. 1Elastic thin plate mode under distributed random load f(x,y,t)

All loading points of distributed random load are completely correlative when$|\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)|=1\text{;}$ all loading points of distributed random load are completely independent when $|\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)|=0$ and$\mathrm{}{x}_{1}\ne {x}_{2}$, ${y}_{1}\ne {y}_{2}\text{;}$ part of loading points of distributed random load is correlative when $|\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)|<1\text{,}$ which is the common situation. Especially, when correlation of distributed random load is non-spatial, the coherence function is $\theta \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=0$ or $\pi$, and $\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=±1$, the auto-spectral density at every spatial point is$\mathrm{}{S}_{f}\left({x}_{1},{y}_{1},\omega \right)={p}^{2}\left({x}_{1},{y}_{1}\right){S}_{f}\left(\omega \right)$, ${S}_{f}\left({x}_{2},{y}_{2},\omega \right)={p}^{2}\left({x}_{2},{y}_{2}\right){S}_{f}\left(\omega \right)$ and the cross-spectral density is ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=p\left({x}_{1},{y}_{1}\right)p\left({x}_{2},{y}_{2}\right){S}_{f}\left(\omega \right)\text{.}$ For example, if coherence function of distributed random load is $\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={e}^{-{\left({x}_{1}-{x}_{2}\right)}^{2}-{\left({y}_{1}-{y}_{2}\right)}^{2}}{e}^{j\omega \left({x}_{1}-{x}_{2}\right)/{y}_{1}{y}_{2}}\text{,}$${S}_{f}\left({x}_{1},{y}_{1},\omega \right)$ is auto-power spectrum density of the load, and its cross-spectral density function is ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={e}^{-{\left({x}_{1}-{x}_{2}\right)}^{2}-{\left({y}_{1}-{y}_{2}\right)}^{2}}{e}^{j\omega \left[\left({x}_{1}-{x}_{2}\right){\left({y}_{1}-{y}_{2}\right)}^{2}\right]}\sqrt{{S}_{f}\left({x}_{1},{y}_{1},\omega \right){S}_{f}\left({x}_{2},{y}_{2},\omega \right)}\text{,}$ the points of distributed random load will be partially coherent, and the statistical features of the load will only be related to the displacement in the space and will satisfy that ${S}_{f}\left(x,y,\omega \right)>0$; and ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={S}_{f}^{*}\left({x}_{2},{y}_{2},{x}_{1},{y}_{1},\omega \right)$ can be obtained when ${x}_{1}={x}_{2}=x$, ${y}_{1}={y}_{2}=y$.

### 2.3. A fast algorithm based on multi-dimensional Legendre orthogonal polynomials

There are five variants in plane-distributed random load ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ applied on thin plates. To determinate $\omega$, hypothetically $x$ and $y$ are defined within a certain range, auto-power spectrum density ${S}_{f}\left(x,y,\omega \right)$ is a continuous function of $x$ and $y$, the real and the imaginary parts of both $\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ and ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ are also continuous functions. After projecting these continuous functions onto orthogonal basis functions, their optimal approximation is found as:

13
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{mn}\left(-\omega \right){H}_{ij}\left(\omega \right)\cdot {Q}_{mnij},$

where ${Q}_{mnji}\left(\omega \right)={\int }_{0}^{a}{\int }_{0}^{b}{\int }_{0}^{a}{\int }_{0}^{b}{S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right)d{x}_{1}d{y}_{1}d{x}_{2}d{y}_{2}$.

14
${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={S}_{f}^{r}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)+i{S}_{f}^{i}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right),$
15
${S}_{f}^{r}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }\sum _{u=1}^{\infty }\sum _{v=1}^{\infty }{L}_{stuv}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){a}_{stuv}\left(\omega \right),$
16
${S}_{f}^{i}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }\sum _{u=1}^{\infty }\sum _{v=1}^{\infty }{L}_{stuv}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){b}_{stuv}\left(\omega \right),$

where ${L}_{mnij}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2}\right)={L}_{m}\left({x}_{1}\right){L}_{n}\left({y}_{1}\right){L}_{i}\left({x}_{2}\right){L}_{j}\left({y}_{2}\right)$ is a four-dimensional Legendre orthogonal polynomial, which is defined as a basis function. Therefore:

17
${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }\sum _{u=1}^{\infty }\sum _{v=1}^{\infty }{L}_{stuv}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right){c}_{stuv}\left(\omega \right),$

where ${c}_{stuv}\left(\omega \right)={a}_{stuv}\left(\omega \right)+i{b}_{stuv}\left(\omega \right)$.

The response cross-power spectrum density between points $\left({x}_{1},{y}_{1}\right)$ and $\left({x}_{\text{2}},{y}_{\text{2}}\right)$ is:

18
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }\sum _{u=1}^{\infty }\sum _{v=1}^{\infty }{c}_{stuv}\left(\omega \right)$
$\bullet \left[\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){H}_{mn}\left(-\omega \right){\int }_{0}^{a}{\int }_{0}^{b}{L}_{st}\left({x}_{1},{y}_{1}\right){\phi }_{mn}\left({x}_{1},{y}_{1}\right)d{x}_{1}d{y}_{1}\right]$
$\bullet \left[\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\phi }_{ij}\left({x}_{1},{y}_{1}\right){H}_{ij}\left(\omega \right){\int }_{0}^{a}{\int }_{0}^{b}{L}_{uv}\left({x}_{2},{y}_{2}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right)d{x}_{2}d{y}_{2}\right].$

Given that:

19
${\mathfrak{R}}_{st}\left({x}_{1},{y}_{1},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){H}_{mn}\left(\omega \right){\int }_{0}^{a}{\int }_{0}^{b}{L}_{st}\left({x}_{1},{y}_{1}\right){\phi }_{mn}\left({x}_{1},{y}_{1}\right)d{x}_{1}d{y}_{1}.$

Hence:

20
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }\sum _{u=1}^{\infty }\sum _{v=1}^{\infty }{c}_{stuv}\left(\omega \right){\mathfrak{R}}_{st}^{*}\left({x}_{1},{y}_{1},\omega \right)·{\mathfrak{R}}_{uv}\left({x}_{2},{y}_{2},\omega \right).$

Fitting ${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ with certain precision requires an order number of ${s}_{0}×{t}_{0}×{u}_{0}×{v}_{0}$, which can be written in matrices as:

21
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$
$=\left\{\begin{array}{ccccc}{\mathfrak{R}}_{11}^{*}\left({x}_{1},{y}_{1}\right){\mathfrak{R}}_{11}\left({x}_{2},{y}_{2}\right)& \dots & {\mathfrak{R}}_{11}^{*}\left({x}_{1},{y}_{1}\right){\mathfrak{R}}_{1{v}_{0}}\left({x}_{2},{y}_{2}\right)& \dots & {\mathfrak{R}}_{{s}_{0}{t}_{0}}^{*}\left({x}_{1},{y}_{1}\right){\mathfrak{R}}_{{u}_{0}{v}_{0}}\left({x}_{2},{y}_{2}\right)\end{array}\right\}$
$·{\left\{\begin{array}{ccccc}{c}_{1111}& \dots & {c}_{111{v}_{0}}& \dots & {c}_{{s}_{0}{t}_{0}{u}_{0}{v}_{0}}\end{array}\right\}}^{T}.$

When ${x}_{1}={x}_{2}$ and ${y}_{1}={y}_{2}$, it yields the auto-power spectrum density of response.

Based on Eq. (19), ${\mathfrak{R}}_{st}\left({x}_{1},{y}_{1},\omega \right)$ represents the structural response in frequency domain under distributed harmonic load, ${L}_{st}\left({x}_{1},{y}_{1}\right){e}^{j\omega t}$. For every discrete ${\omega }_{i}$, the orthogonal polynomial’s coefficients for fitting ${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},{\omega }_{i}\right)$ can be computed and structural response under harmonic load can be constructed with the coefficients and two-dimensional basis functions. According to the superposition of linear structures, linearly superimposing the responses of the structure under distributed harmonic load can finally lead to response power spectrum density of the structure. In this way, structural responses under distributed random load can be represented as responses in the frequency domain of the structure under determinate distributed harmonic load. According to the deduction above, this method takes all coupling items of the vibration frequency into consideration, which therefore makes it an accurate algorithm. The method is named pseudo-excitation method for two-dimensional distributed random loads.

If the distributed random load is linearly distributed at the position $y={y}_{0}$ of elastic thin plate, the cross-spectral density function of any two points, $\left({x}_{1},{y}_{1}\right)$ and $\left({x}_{2},{y}_{\text{2}}\right)$, can be described as ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)\delta \left({y}_{1}-{y}_{0}\right)\delta \left({y}_{2}-{y}_{0}\right)$. The algorithm above can be rewritten as:

22
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{mn}\left(-\omega \right){H}_{ij}\left(\omega \right)$
$\bullet {\int }_{0}^{a}{\int }_{0}^{b}{\int }_{0}^{a}{\int }_{0}^{b}\left(\begin{array}{c}{S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)\delta \left({y}_{1}-{y}_{0}\right)\delta \left({y}_{2}-{y}_{0}\right)\\ \bullet {\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right)d{x}_{1}d{y}_{1}d{x}_{2}d{y}_{2}\end{array}\right)$
$=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\left(\begin{array}{c}{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{mn}\left(-\omega \right){H}_{ij}\left(\omega \right)\\ \bullet {\int }_{0}^{a}{\int }_{0}^{b}{S}_{f}\left({x}_{1},{y}_{0},{x}_{2},{y}_{0},\omega \right){\phi }_{mn}\left({x}_{1},{y}_{0}\right){\phi }_{ij}\left({x}_{2},{y}_{0}\right)d{x}_{1}d{x}_{2}\end{array}\right).$

Then the variant number of ${S}_{f}\left({x}_{1},{y}_{0},{x}_{2},{y}_{0},\omega \right)$ is reduced to three. For a determinate $\omega$, within a certain range, it can be obtained as:

23
${S}_{f}\left({x}_{1},{y}_{0},{x}_{2},{y}_{0},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }{L}_{st}\left({x}_{1},{x}_{2},\omega \right){c}_{st}\left(\omega \right),$
24
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$
$=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\left(\begin{array}{c}{\phi }_{mn}\left({x}_{1},{y}_{1}\right){\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{mn}\left(-\omega \right){H}_{ij}\left(\omega \right)\\ \bullet {\int }_{0}^{a}{\int }_{0}^{b}\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }{L}_{st}\left({x}_{1},{x}_{2},\omega \right){c}_{st}\left(\omega \right){\phi }_{mn}\left({x}_{1},{y}_{0}\right){\phi }_{ij}\left({x}_{2},{y}_{0}\right)d{x}_{1}d{x}_{2}\end{array}\right)$
$=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }{c}_{st}\left(\omega \right)\left[\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){H}_{mn}\left(-\omega \right){\int }_{0}^{a}{L}_{s}\left({x}_{1}\right){\phi }_{mn}\left({x}_{1},{y}_{0}\right)d{x}_{1}\right]$
$·\left[\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\phi }_{ij}\left({x}_{2},{y}_{2}\right){H}_{ij}\left(\omega \right){\int }_{0}^{b}{L}_{t}\left({x}_{2}\right){\phi }_{ij}\left({x}_{2},{y}_{0}\right)d{x}_{2}\right].$

Finally, it can be simplified as:

25
${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=\sum _{s=1}^{\infty }\sum _{t=1}^{\infty }{c}_{st}\left(\omega \right){\mathfrak{R}}_{st}^{*}\left({x}_{1},{y}_{1},\omega \right)·{\mathfrak{R}}_{uv}\left({x}_{2},{y}_{2},\omega \right),$

where:

26
${\mathfrak{R}}_{st}\left({x}_{1},{y}_{1},\omega \right)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\phi }_{mn}\left({x}_{1},{y}_{1}\right){H}_{mn}\left(\omega \right){\int }_{0}^{a}{L}_{s}\left({x}_{1}\right){\phi }_{mn}\left({x}_{1},{y}_{0}\right)d{x}_{1}.$

Based on the Eq. (26), ${\mathfrak{R}}_{st}\left({x}_{1},{y}_{1},\omega \right)$ represents structural response in frequency domain under distributed harmonic load, ${L}_{s}\left(x\right){e}^{j\omega t}\text{.}$ Similarly, for every discrete ${\omega }_{i}\text{,}$ the orthogonal polynomials’ coefficients for fitting ${S}_{y}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},{\omega }_{i}\right)$ can be computed and structural response under harmonic load can be constructed by one-dimensional basis functions. According to the superposition of linear structures, linearly superimposing responses of the structure under linearly distributed harmonic load can finally lead to response power spectrum density of the structure. Simplified computation under linearly distributed random load is an extreme case of that under plane-distributed random load.

### 2.4. Comparison with regular algorithm in the computation speed

Computation of response in frequency domain under two-dimensional distributed random load is conducted with CQC, SRSS and the fast algorithm proposed in this paper. It is assumed that ${\phi }_{mn}\left(x,y\right){\phi }_{mn}\left(x,y\right)$ and ${H}_{mn}\left(\omega \right)$ are both known. Moreover, the former $q×r$ orders of modes are involved in the computation, and the three methods have the same number of discrete points in frequency domain. The algorithm employed to integrate is Gauss-Legendre integral formula as ${\int }_{a}^{b}f\left(x\right)dx=\sum _{i=1}^{n}{A}_{i}f\left({t}_{i}\right)$ [24], where ${t}_{i}$ represents Gauss points in the interval $\left[a,b\right]$ and ${A}_{i}$ are Gauss coefficients corresponding to ${t}_{i}$. Based on this formula, Gauss double integral formula is established as ${\int }_{a}^{b}{\int }_{c}^{d}f\left(x,y\right)dx=\sum _{i=1}^{{i}_{0}}\sum _{j=1}^{{j}_{0}}{A}_{i}{B}_{j}f\left({x}_{i},{y}_{j}\right)$, where ${x}_{i}$ and ${y}_{j}$ are Gauss points in intervals $\left[a,b\right]$ and $\left[c,d\right]$ respectively; and ${A}_{i}$ and ${B}_{j}$ are Gauss coefficients corresponding to ${x}_{i}$ and${y}_{j}\text{,}$ respectively. Furthermore, Gauss quadruple integral formula is ${\int }_{a}^{b}{\int }_{c}^{d}{\int }_{e}^{f}{\int }_{g}^{h}f\left(x,y,u,v\right)dxdydudv=\sum _{m=1}^{{m}_{0}}\sum _{n=1}^{{n}_{0}}\sum _{i=1}^{{i}_{0}}\sum _{j=1}^{{j}_{0}}{A}_{m}{B}_{n}{C}_{i}{D}_{j}f\left({x}_{m},{y}_{n},{u}_{i},{v}_{j}\right)\text{,}$ given that the number of Gauss points in every interval is the same $n$.

Low computation efficiency mainly results from large amount of multiplication. The total number of multiplication for every discrete point in frequency domain with CQC method is ${q}^{2}{r}^{2}\left(4{n}^{4}+4\right)$; it’s $qr\left(4{n}^{4}+4\right)$ with SRSS method; with fast algorithm it is $m\left[qr\left(2{n}^{2}+4\right)\right]$, where $m$ is the order of orthogonal polynomials required for fitting ${S}_{f}\left({x}_{1},{x}_{2},{y}_{1},{y}_{2},\omega \right)$.

Table 1Comparison of the calculated amount of three methods (unit: number of multiplication)

 Parameters CQC SRSS Fast algorithm $q=$5, $r=$5, $n=$ 5, $m=$ 16 1565000 62600 21600 $q=$ 5, $r=$ 5, $n=$ 10, $m=$ 16 2500250 1000100 81600 $q=$ 10, $r=$ 10, $n=$ 10, $m=$ 81 4.0004E+8 4.0004E+6 1.6524E+6 $q=$ 10, $r=$ 10, $n=$ 20, $m=$ 81 6.40004E+9 6.40004E+7 6.5124E+6

As revealed in Table 1, the fast algorithm requires far less multiplication than CQC; and the smaller $m$ is, the more efficient the fast algorithm becomes. Generally speaking, the number of Gauss points that computation requires grows larger as $m$ gets bigger. Consequently, the fast algorithm is always more efficient than CQC and SRSS.

### 2.5. Computation instruction

Computation steps are as follows:

1) Calculate ${\phi }_{mn}\left(x\right)$ and ${H}_{mn}\left(\omega \right)$. Their analytical solutions can be employed for typical structures. As for complex structures, ${H}_{mn}\left(\omega \right)$ can be obtained by finite element method.

2) Decompose ${S}_{f}\left({x}_{1},{y}_{1},{x}_{\text{2}},{y}_{2},\omega \right)$ in orthogonal domain and gain coefficients for fitting orthogonal polynomials.

3) Compute ${\mathfrak{R}}_{n}\left({x}_{1},{y}_{1},\omega \right)$ for every discrete points in frequency domain.

4) Conduct linear superposition based on Eq. (21) or (25), and obtain final results.

In summary, the fast algorithm for solving the problem of random response of elastic thin plate under distributed random load has two key steps. First, it projects the function of the cross power spectrum function onto orthogonal basis of multi-dimension Legendre polynomial and calculates the projection coefficient. Second, it calculates harmonic response under harmonic distributed load formed by the orthogonal basis function. In fact, it converts random response problem into a harmonic response problem, and harmonic response calculation is much easier than random response calculation. The fast algorithm applies to regular thin elastic plate and simple elastic structure, and it will be suitable for complicated structure after modification.

## 3. Simulation example

An elastic rectangular thin plate of homogeneous material with four simply supported edges is taken as the structure model. It has 1 m length, 0.6 m width, and 0.02 m thickness. The elastic modulus $E=\text{2.1×}{\text{10}}^{\text{11}}\text{N/}{\text{m}}^{\text{2}}$, density $\rho =\text{7800}\text{}\text{kg/}{\text{m}}^{\text{3}}$, Poisson's ratio is 0.3, modal damping ratio is 0.02, and number of finite elements is 240. The natural frequencies of the model are listed in Table 2. $f\left(x,y,t\right)$ is a two-dimensional distributed random load, and ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)$ represents the cross-spectral density of any two points.

Table 2First eight orders of the plate model

 Orders 1 2 3 4 5 6 7 8 Natural frequencies (Hz) 183.88 326.85 565.75 588.29 720.27 898.67 941.63 1252.02

Example 1: ${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)=1.$

In this example, computation results of the fast algorithm and Nastran are compared, and efficiency is assessed. Due to incapability of defining coherence of distributed random load, Nastran can only compute distributed random load of complete coherence. Therefore, distributed random load in this example is set as complete coherence, and its power spectrum density is unit flat-spectrum.

The first eight orders are taken into computation. The order of orthogonal fitting is 1, and the number of Gauss points is 10. Frequency span is 1000 Hz, and frequency resolution is 1 Hz. Time consumed for computation is listed in Table 3. The power spectrum density line of acceleration at Node 125 of finite element model is displayed in Figures 2 and 3.

Table 3Comparison of time consumed

 Algorithm CQC SRSS Fast algorithm Nastran Time consumed (s) 19505.4 3733.3 485.047 478.6

Response at Node 125 of finite element model is computed with CQC, SRSS, the fast algorithm method and Nastran, under the same other conditions. As shown in Table 3, the fast algorithm is far more efficient than CQC and SRSS, and equivalent to that of Nastran (time consumed for computation includes time for computing coefficients of orthogonal polynomials).

Fig. 2The acceleration power spectrum density response result at Node125 with the finite element method

Fig. 3Comparison of the acceleration power spectrum density result at Node 125 with four methods

Figures 2 and 3 show that result of the fast algorithm is almost the same as that of Nastran, except the slight difference ranging from 900 Hz to 1000 Hz. It is due to different modal truncation. They also manifest that result of SRSS, displayed as pink dotted line, has some errors which are quite large at some points. These errors get even larger when damping ratio gets higher or modes become closer. The errors between the acceleration power spectrum density at natural frequencies produced by CQC, SRSS, the fast algorithm and Nastran are listed in Table 4. According to Table 4, the fast algorithm has the same result as CQC, and it is fairly close to the result of Nastran, while SRSS result shows great error.

Table 4Comparison between CQC, SRSS, the fast algorithm, and Nastran

 Error types Relative error between CQC and Nastran Relative error between SRSS and Nastran Relative error between fast algorithm and Nastran ${f}_{1}=184Hz$ 0.827 % 1.853 % 0.827 % ${f}_{2}=588Hz$ 1.112 % 13.922 % 1.112 %

Example 2:

${S}_{f}\left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={e}^{-\left[{\left({x}_{1}-{x}_{2}\right)}^{2}+{\left({y}_{1}-{y}_{2}\right)}^{2}\right]}{e}^{j\left[\frac{{x}_{1}-{x}_{2}}{{y}_{1}{y}_{2}}\right]}\sqrt{\left(1+{x}_{1}^{2}+{y}_{1}^{2}\right)\left(1+{x}_{2}^{2}+{y}_{2}^{2}\right)}.$

In this example, the distributed random excitation is partially coherent, and the cross power spectrum density of distributed random load is defined as above. Derived from formula of Example 2, it is obtained that: $\rho \left({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \right)={e}^{-\left[{\left({x}_{1}-{x}_{2}\right)}^{2}+{\left({y}_{1}-{y}_{2}\right)}^{2}\right]}{e}^{j\left[\left({x}_{1}-{x}_{2}\right)/{y}_{1}{y}_{2}\right]}$, and ${S}_{f}\left({x}_{1},{y}_{1},\omega \right)=1+{x}_{1}^{2}+{y}_{1}^{2}$.

Any two points are partially coherent when $\left({x}_{1},{y}_{1}\right)\ne \left({x}_{2},{y}_{2}\right)$. First eight orders of modes are taken into computation. The orders of orthogonal fitting are $4×4×4×4$, and the number of Gauss points is 10. The start frequency is 100 Hz, and stop frequency 400 Hz. Under the same other conditions, computation is conducted with CQC, SRSS, the new method and Nastran. The efficiency of computation is demonstrated in Table 5.

Table 5Comparison of three methods in computation

 Algorithm CQC SRSS Fast algorithm Nastran Time consumed (s) 523847 40104.8 11231.5 Not applicable

Table 5 further proves excellent efficiency of the fast algorithm. Fig. 4 illustrates the power spectrum density of response at Node 125 computed with these three methods, and errors are listed in Table 6. The largest error between SRSS and CQC is above 20 %, because of leaving cross items of mode shapes out. The fast algorithm possesses the same high precision as CQC, as shown in Figure 4.

Fig. 4Comparison of the acceleration power spectrum density result at Node 125 with the three methods

Table 6Computational errors of different methods

 Error type Relative error between CQC and SRS Relative error betweenfast algorithm and SRSS Relative error between fast algorithm and CQC Largest relative error 24.97 % 24.97 % 6.357E-8

## 4. Conclusions

This paper proposes a fast algorithm for the calculation of structural response power spectrum density of thin plates under either one or two dimensional distributed random load. According to the simulation results of completely coherent two-dimensional distributed random load and partially coherent two-dimensional distributed random load, it is concluded that the proposed algorithm is fast as well as precise. Its efficiency is far higher than that of both CQC and SRSS. It is programmed for computation of structural responses in the frequency domain under two-dimensional distributed random load that is completely coherent, partially coherent or altogether incoherent. The efficiency is highly improved, and the precision is also guaranteed. It facilitates the computation of responses under distributed random load in engineering. In addition, it establishes a basis for the subsequent identification of distributed random load.

#### References

• Zhang Fang, Qin Yuantian, Deng Jihong. Research of identification technology of dynamic load distributed on the structure. Journal of Vibration Engineering, Vol. 19, Issue 1, 2006, p. 81-85.
• Cameron W., Coates Thamburaj P. Inverse method using finite strain measurement to determine flight load distribution functions. Journal of Aircraft, Vol. 45, Issue 2, 2008, p. 366-370.
• Coates C. W., Thamburaj P., Kim C. An inverse method for selection of Fourier coefficients for flight load identification. Proceedings of the 46th Annual AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Mechanics Conference, AIAA, Reston, 2005, p. 18-21.
• Lin J. H., Guo X. L., Zhi H., Howson W. P. Computer simulation of structural random loading identification. Computers and Structures, Vol. 79, Issue 4, 2001, p. 375-387.
• K. Soyluk. Comparison on random vibration methods for multi-support seismic excitation analysis of long-span bridges. Engineering Structures, Vol. 26, Issue 11, 2004, p. 1776-1778.
• Lin J. H., Williams F. W. Computation and analysis of multi-excitation random seismic responses. Engineering Computations, Vol. 9, Issue 5, 1992, p. 561-574.
• Xue S. D., Li M. H., Cao Z., Zhang Y. G. Random vibration analysis of lattice shell subjected to multi-dimensional earthquake inputs. Proc. Int. Conf. on Advances in Structural Dynamics, Elsevier, 2000, p. 777-784.
• Jiang Jinhui, Chen Guoping, Zhang Fang. Fast algorithm for structural frequency domain response subjected to one-dimension distributed random load. Journal of Vibration and Shock, Vol. 29, Issue 8, 2010, p. 55-59.
• P. Sniady, S. Biernat, R. Sieniawska, S. Zukowsky. Vibrations of the beam due to a load moving with stochastic velocity. Probabilist. Eng. Mech., Vol. 16, Issue 1, 2001, p. 53-59.
• S. S. Law, X. Q. Zhu. Bridge dynamic responses due to road surface roughness and braking of vehicle. J. Sound Vibration, Vol. 282, 2005, p. 805-830.
• W. X. Zhong. Review of a high-efficiency algorithm series for structural random responses. Prog. Nat. Sci., Vol. 6, Issue 3, 1996, p. 257-268.
• M. F. Dimentberg, B. Ryzhik, L. Sperling. Random vibrations of a damped rotating shaft. Journal of Sound and Vibration, Vol. 279, 2005, p. 275-284.
• Zhou X. Y., Yu R. F. Mode superposition method of non stationary seismic responses for non classically damped linear systems. Journal of Earthquake Engineering, Vol. 12, Issue 3, 2008, p. 597-641.
• J. H. Lin, Y. H. Zhang, Y. Zhao. Pseudo excitation method and some recent developments. Proceedings of the Twelfth East Asia-Pacific Conference on Structural Engineering and Construction, Hong Kong, January, 2011, p. 2453-2458.
• A. D. Kiureghian, A. Neuenhofer. Response spectrum method for multi-support seismic excitations. Earthquake Engrg. and Struct. Dynamics, Vol. 21, Issue 8, 1992, p. 713-740.
• H. Z. Ernesto, E. H. Vanmarcke. Seismic random vibration analysis of multi-support structural systems. ASCE J. Eng. Mech., Vol. 120, 1994, p. 1107-1128.
• Wei Guo, Zhiwu Yu. Pseudo excitation method based on multiple degree of freedom modal equation. Advanced Materials Research, 2012, p. 317-320.
• Yu R. F., Wang L., Liu H., Zhou X. Y. The CQC3 method for multi component seismic responses of non-classically damped irregular structures. American Structures Congress, New York, 2005.
• Q. Gao, J. H. Lin, W. X. Zhong, W. P. Howson, F. W. Williams. Propagation of partially coherent non-stationary random waves in viscoelastic layered half-space. Soil Dynamics and Earthquake Engineering, Vol. 28, Issue 4, 2008, p. 305-320.
• Zhichao Zhang, Yahui Zhang, Jiahao Lin, Yan Zhao, W. P. Howson, F. W. Williams. Random vibration of a train traversing a bridge subjected to traveling seismic waves. Engineering Structures, Vol. 33, Issue 12, 2011, p. 3546-3558.
• Y. L. Xu, N. Zhang, H. Xia. Vibration of coupled train and cable-stayed bridge systems in cross winds. Eng. Struct., Vol. 26, Issue 10, 2004, p. 1389-1406.
• Li J., Liao S. T. Response analysis of stochastic parameter structures under non-stationary random excitation. Comput. Mech., Vol. 27, Issue 1, 2001, p. 61-68.
• Lin J. H., Zhao Y., Zhao Y. H. Accurate and highly efficient algorithms for structural stationary/non-stationary random responses. Comput. Methods Appl. Mech. Eng., Vol. 191, 2001, p. 103-111.
• Sansone G. Orthogonal Functions. English Ed., New York, Dover, 1991.
• Charles C. Dyer, Peter S. S. Ip. An Elementary Introduction to Scientific Computing. Division of Physical Sciences, University of Toronto at Scarborough, January 2000.

16 May 2013
Accepted
04 September 2013
Published
30 September 2013
Keywords
elastic thin plate