Abstract
The calculation of response power spectrum of elastic thin plates, the common loadcarrying 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 multidimension 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 [48]. In addition, the calculation of dynamic random response focuses on the structure due to multipoints excitation [913]. 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 abovementioned 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) [1518], SRSS (Square Root of the Sum of Squares) [1921], and PEM (PseudoExcitation 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):
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}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ and ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ represent crossspectrum density of response and crossspectrum density of excitation between points $({x}_{1},{y}_{1})$ and $({x}_{2},{y}_{2})\text{,}$ respectively. ${\phi}_{mn}({x}_{1},{y}_{1})$ is mode shape of mode $mn$ at point $({x}_{1},{y}_{1})$, ${H}_{ij}\left(\omega \right)$ is modal response frequency function of mode $ij$, and $[0,a]$, $[0,b]$ 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):
$\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 mulitinputmultioutput stationary/nonstationary random vibration response problem since 1980s, thus it is a CQC algorithm with high efficiency and simple procedure [913]. 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 onedimensional 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 linedistributed (onedimensional) or planedistributed (twodimensional) random load on the thin plate. First, onedimensional Legendre orthogonal polynomials are introduced to establish multidimensional ones. Second, the function of crosspower spectrum density, which is complicated to distribute on thin plate structures, is decomposed in orthogonal domain of multidimensional 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. Multidimensional Legendre orthogonal polynomials
Located within a real span $[1,1]$ and with weighting function $\rho \left(x\right)=1\text{,}$ Legendre orthogonal polynomials [24, 25] can be described as follows:
In $[0,s]$, 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}_{n1}\left(x\right),\end{array}\right.\mathrm{}\mathrm{}x\in \left[0\text{\hspace{0.17em}},\text{\hspace{0.17em}}s\right],n=\mathrm{1,2},\dots \mathrm{}.$
Based on the above onedimensional Legendre orthogonal polynomials, twodimensional orthogonal polynomials can be established as ${L}_{mn}(x,y)={L}_{m}\left(x\right){L}_{n}\left(y\right)$. According to the orthogonality of onedimensional Legendre orthogonal polynomials, the orthogonality of twodimensional ones is easily proved as follows:
It can also be proved that the established twodimensional Legendre orthogonal polynomials possess the characteristics of parity and recurrence. As a result, $\left\{{L}_{mn}\right(x,y\left)\right\}$ constructs a series of twodimensional orthogonal basis functions. Furthermore, multidimensional orthogonal polynomials can be constructed similarly as below:
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, twodimensional and fourdimensional orthogonal polynomials are involved.
If a continuous function $f(x,y)$ is square integrable ($f\left(x,y\right)\in {L}_{2}\&f(x,y)\in C$, where ${L}_{2}$, $C$ is square integrable function space and continuous function space, respectively), it can be expressed via the expansion with twodimensional Legendre orthogonal polynomials in the following form:
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:
Let $\epsilon $ be an arbitrary positive number, $M$, $N$ can be found to satisfy the following inequality:
$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):
Here, ${S}_{f}(x,y,\omega )$ is stationary autospectral density of the distributed random load at point $(x,y)$, which is function of space and frequency variable. For example, Davenport spectrum for fluctuating winds velocity:
$\rho ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ is a coherence function and satisfies Eq. (11) below:
The argument $\left(\theta \right)$ shows as:
If the autopower spectrum density and the coherence function of every point are given, the crosspower 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$\left\rho \right({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \left)\right=1\text{;}$ all loading points of distributed random load are completely independent when $\left\rho \right({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \left)\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 $\left\rho \right({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega \left)\right<1\text{,}$ which is the common situation. Especially, when correlation of distributed random load is nonspatial, the coherence function is $\theta ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )=0$ or $\pi $, and $\rho ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )=\pm 1$, the autospectral density at every spatial point is$\mathrm{}{S}_{f}({x}_{1},{y}_{1},\omega )={p}^{2}({x}_{1},{y}_{1}){S}_{f}\left(\omega \right)$, ${S}_{f}({x}_{2},{y}_{2},\omega )={p}^{2}({x}_{2},{y}_{2}){S}_{f}\left(\omega \right)$ and the crossspectral density is ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )=p({x}_{1},{y}_{1})p({x}_{2},{y}_{2}){S}_{f}\left(\omega \right)\text{.}$ For example, if coherence function of distributed random load is $\rho ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )={e}^{{({x}_{1}{x}_{2})}^{2}{({y}_{1}{y}_{2})}^{2}}{e}^{j\omega ({x}_{1}{x}_{2})/{y}_{1}{y}_{2}}\text{,}$${S}_{f}({x}_{1},{y}_{1},\omega )$ is autopower spectrum density of the load, and its crossspectral density function is ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )={e}^{{({x}_{1}{x}_{2})}^{2}{({y}_{1}{y}_{2})}^{2}}{e}^{j\omega \left[\right({x}_{1}{x}_{2}\left){({y}_{1}{y}_{2})}^{2}\right]}\sqrt{{S}_{f}({x}_{1},{y}_{1},\omega ){S}_{f}({x}_{2},{y}_{2},\omega )}\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}(x,y,\omega )>0$; and ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )={S}_{f}^{*}({x}_{2},{y}_{2},{x}_{1},{y}_{1},\omega )$ can be obtained when ${x}_{1}={x}_{2}=x$, ${y}_{1}={y}_{2}=y$.
2.3. A fast algorithm based on multidimensional Legendre orthogonal polynomials
There are five variants in planedistributed random load ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ applied on thin plates. To determinate $\omega $, hypothetically $x$ and $y$ are defined within a certain range, autopower spectrum density ${S}_{f}(x,y,\omega )$ is a continuous function of $x$ and $y$, the real and the imaginary parts of both $\rho ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ and ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ are also continuous functions. After projecting these continuous functions onto orthogonal basis functions, their optimal approximation is found as:
where ${Q}_{mnji}\left(\omega \right)={\int}_{0}^{a}{\int}_{0}^{b}{\int}_{0}^{a}{\int}_{0}^{b}{S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega ){\phi}_{mn}({x}_{1},{y}_{1}){\phi}_{ij}({x}_{2},{y}_{2})d{x}_{1}d{y}_{1}d{x}_{2}d{y}_{2}$.
where ${L}_{mnij}({x}_{1},{y}_{1},{x}_{2},{y}_{2})={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 fourdimensional Legendre orthogonal polynomial, which is defined as a basis function. Therefore:
where ${c}_{stuv}\left(\omega \right)={a}_{stuv}\left(\omega \right)+i{b}_{stuv}\left(\omega \right)$.
The response crosspower spectrum density between points $({x}_{1},{y}_{1})$ and $({x}_{\text{2}},{y}_{\text{2}})$ is:
$\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}({x}_{1},{y}_{1}){H}_{ij}\left(\omega \right){\int}_{0}^{a}{\int}_{0}^{b}{L}_{uv}({x}_{2},{y}_{2}){\phi}_{ij}({x}_{2},{y}_{2})d{x}_{2}d{y}_{2}\right].$
Given that:
Hence:
Fitting ${S}_{y}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ with certain precision requires an order number of ${s}_{0}\times {t}_{0}\times {u}_{0}\times {v}_{0}$, which can be written in matrices as:
$=\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\}$
$\xb7{\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 autopower spectrum density of response.
Based on Eq. (19), ${\mathfrak{R}}_{st}({x}_{1},{y}_{1},\omega )$ represents the structural response in frequency domain under distributed harmonic load, ${L}_{st}({x}_{1},{y}_{1}){e}^{j\omega t}$. For every discrete ${\omega}_{i}$, the orthogonal polynomial’s coefficients for fitting ${S}_{y}({x}_{1},{y}_{1},{x}_{2},{y}_{2},{\omega}_{i})$ can be computed and structural response under harmonic load can be constructed with the coefficients and twodimensional 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 pseudoexcitation method for twodimensional distributed random loads.
If the distributed random load is linearly distributed at the position $y={y}_{0}$ of elastic thin plate, the crossspectral density function of any two points, $({x}_{1},{y}_{1})$ and $({x}_{2},{y}_{\text{2}})$, can be described as ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )={S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )\delta ({y}_{1}{y}_{0})\delta ({y}_{2}{y}_{0})$. The algorithm above can be rewritten as:
$\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}({x}_{1},{y}_{0},{x}_{2},{y}_{0},\omega ){\phi}_{mn}({x}_{1},{y}_{0}){\phi}_{ij}({x}_{2},{y}_{0})d{x}_{1}d{x}_{2}\end{array}\right).$
Then the variant number of ${S}_{f}({x}_{1},{y}_{0},{x}_{2},{y}_{0},\omega )$ is reduced to three. For a determinate $\omega $, within a certain range, it can be obtained as:
$=\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}({x}_{1},{x}_{2},\omega ){c}_{st}\left(\omega \right){\phi}_{mn}({x}_{1},{y}_{0}){\phi}_{ij}({x}_{2},{y}_{0})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]$
$\xb7\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:
where:
Based on the Eq. (26), ${\mathfrak{R}}_{st}({x}_{1},{y}_{1},\omega )$ 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}({x}_{1},{y}_{1},{x}_{2},{y}_{2},{\omega}_{i})$ can be computed and structural response under harmonic load can be constructed by onedimensional 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 planedistributed random load.
2.4. Comparison with regular algorithm in the computation speed
Computation of response in frequency domain under twodimensional distributed random load is conducted with CQC, SRSS and the fast algorithm proposed in this paper. It is assumed that ${\phi}_{mn}(x,y){\phi}_{mn}(x,y)$ and ${H}_{mn}\left(\omega \right)$ are both known. Moreover, the former $q\times 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 GaussLegendre 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 $[a,b]$ 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(x,y)dx=\sum _{i=1}^{{i}_{0}}\sum _{j=1}^{{j}_{0}}{A}_{i}{B}_{j}f({x}_{i},{y}_{j})$, where ${x}_{i}$ and ${y}_{j}$ are Gauss points in intervals $[a,b]$ and $[c,d]$ 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(x,y,u,v)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({x}_{m},{y}_{n},{u}_{i},{v}_{j})\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}(4{n}^{4}+4)$; it’s $qr(4{n}^{4}+4)$ with SRSS method; with fast algorithm it is $m\left[qr\right(2{n}^{2}+4\left)\right]$, where $m$ is the order of orthogonal polynomials required for fitting ${S}_{f}({x}_{1},{x}_{2},{y}_{1},{y}_{2},\omega )$.
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}({x}_{1},{y}_{1},{x}_{\text{2}},{y}_{2},\omega )$ in orthogonal domain and gain coefficients for fitting orthogonal polynomials.
3) Compute ${\mathfrak{R}}_{n}({x}_{1},{y}_{1},\omega )$ 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 multidimension 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\xd7}{\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(x,y,t)$ is a twodimensional distributed random load, and ${S}_{f}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )$ represents the crossspectral 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}({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )=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 flatspectrum.
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:
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 ({x}_{1},{y}_{1},{x}_{2},{y}_{2},\omega )={e}^{[{({x}_{1}{x}_{2})}^{2}+{({y}_{1}{y}_{2})}^{2}]}{e}^{j\left[\right({x}_{1}{x}_{2})/{y}_{1}{y}_{2}]}$, and ${S}_{f}({x}_{1},{y}_{1},\omega )=1+{x}_{1}^{2}+{y}_{1}^{2}$.
Any two points are partially coherent when $({x}_{1},{y}_{1})\ne ({x}_{2},{y}_{2})$. First eight orders of modes are taken into computation. The orders of orthogonal fitting are $4\times 4\times 4\times 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 between fast algorithm and SRSS  Relative error between fast algorithm and CQC 
Largest relative error  24.97 %  24.97 %  6.357E8 
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 twodimensional distributed random load and partially coherent twodimensional 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 twodimensional 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. 8185.

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. 366370.

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. 1821.

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. 375387.

K. Soyluk. Comparison on random vibration methods for multisupport seismic excitation analysis of longspan bridges. Engineering Structures, Vol. 26, Issue 11, 2004, p. 17761778.

Lin J. H., Williams F. W. Computation and analysis of multiexcitation random seismic responses. Engineering Computations, Vol. 9, Issue 5, 1992, p. 561574.

Xue S. D., Li M. H., Cao Z., Zhang Y. G. Random vibration analysis of lattice shell subjected to multidimensional earthquake inputs. Proc. Int. Conf. on Advances in Structural Dynamics, Elsevier, 2000, p. 777784.

Jiang Jinhui, Chen Guoping, Zhang Fang. Fast algorithm for structural frequency domain response subjected to onedimension distributed random load. Journal of Vibration and Shock, Vol. 29, Issue 8, 2010, p. 5559.

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. 5359.

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. 805830.

W. X. Zhong. Review of a highefficiency algorithm series for structural random responses. Prog. Nat. Sci., Vol. 6, Issue 3, 1996, p. 257268.

M. F. Dimentberg, B. Ryzhik, L. Sperling. Random vibrations of a damped rotating shaft. Journal of Sound and Vibration, Vol. 279, 2005, p. 275284.

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. 597641.

J. H. Lin, Y. H. Zhang, Y. Zhao. Pseudo excitation method and some recent developments. Proceedings of the Twelfth East AsiaPacific Conference on Structural Engineering and Construction, Hong Kong, January, 2011, p. 24532458.

A. D. Kiureghian, A. Neuenhofer. Response spectrum method for multisupport seismic excitations. Earthquake Engrg. and Struct. Dynamics, Vol. 21, Issue 8, 1992, p. 713740.

H. Z. Ernesto, E. H. Vanmarcke. Seismic random vibration analysis of multisupport structural systems. ASCE J. Eng. Mech., Vol. 120, 1994, p. 11071128.

Wei Guo, Zhiwu Yu. Pseudo excitation method based on multiple degree of freedom modal equation. Advanced Materials Research, 2012, p. 317320.

Yu R. F., Wang L., Liu H., Zhou X. Y. The CQC3 method for multi component seismic responses of nonclassically 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 nonstationary random waves in viscoelastic layered halfspace. Soil Dynamics and Earthquake Engineering, Vol. 28, Issue 4, 2008, p. 305320.

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. 35463558.

Y. L. Xu, N. Zhang, H. Xia. Vibration of coupled train and cablestayed bridge systems in cross winds. Eng. Struct., Vol. 26, Issue 10, 2004, p. 13891406.

Li J., Liao S. T. Response analysis of stochastic parameter structures under nonstationary random excitation. Comput. Mech., Vol. 27, Issue 1, 2001, p. 6168.

Lin J. H., Zhao Y., Zhao Y. H. Accurate and highly efficient algorithms for structural stationary/nonstationary random responses. Comput. Methods Appl. Mech. Eng., Vol. 191, 2001, p. 103111.

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.
About this article
The research described in this paper was supported by the Fundamental Research Funds for the Central Universities No. NS2012080, the Aeronautical Science Foundation of China under Grant No. 2012ZA52001, and Research Fund for the Doctoral Program of Higher Education of China No. 20123218120005.