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

2.1. Multi-dimensional 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}_{n-1}\left(x\right),\end{array}\right.\mathrm{}\mathrm{}x\in \left[0,s\right],n=\mathrm{1,2},\dots \mathrm{}.$

Based on the above one-dimensional Legendre orthogonal polynomials, two-dimensional 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 one-dimensional Legendre orthogonal polynomials, the orthogonality of two-dimensional ones is easily proved as follows:

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}\right(x,y\left)\right\}$ constructs a series of two-dimensional orthogonal basis functions. Furthermore, multi-dimensional 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, two-dimensional and four-dimensional 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 two-dimensional 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 auto-spectral 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 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$\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 non-spatial, 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 auto-spectral 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 cross-spectral 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 auto-power spectrum density of the load, and its cross-spectral 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 multi-dimensional Legendre orthogonal polynomials

There are five variants in plane-distributed 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, auto-power 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 four-dimensional 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 cross-power spectrum density between points $(x_1,y_1)$ and $(x_{\text{2}},y_{\text{2}})$ is:

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:

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}(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 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, $(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:

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:

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 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}(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 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 $[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(4n^4+4)$; it’s $qr(4n^4+4)$ with SRSS method; with fast algorithm it is $m\left[qr\right(2n^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 )$.

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