Stochastic Finite Element Method in Mechanical Vibration

1. Introduction

Material properties, geometry parameters and applied loads of the structure are assumed to be stochastic. Although the finite element method analysis of complicated structures has become a generally widespread and accepted numerical method in the world, regarding the given factors as constants can not apparently correspond to the reality of a structure.

The direct Monte Carlo simulation of the stochastic finite element method(DSFEM) requires a large number of samples, which requires much calculation time and occupies much computer storage space [1]. Monte Carlo simulation by applying the Neumann expansion (NSFEM) enhances computational efficiency and saves storage in such a way that the NSFEM combined with Monte Carlo simulation enhances the finite element model advantageously [2]. The preconditioned Conjugate Gradient method (PCG) applied in the calculation of stochastic finite elements can also enhance computational accuracy and efficiency [3]. The TSFEM assumes that random variables are dealt with by Taylor expansion around mean values and is obtained by appropriate mathematical treatment [4, 14]. According to first-order or second-order perturbation methods, calculation formulas can be obtained [2, 5, 6,8, 9, 13, 15, 16]. The result is called the PSFEM and has been adopted by many scholars.

The PSFEM is often applied in dynamic analysis of structures and the second- order perturbation technique has been proved to be accurate and efficient. Dynamic reliability of a large frame is calculated by the SFEM and response sensitivity is formulated in the context of stiffness and mass matrix condensation [7]. Nonlinear structural dynamics are developed by the PSFEM. Nonlinearities due to material and geometrical effects have also been included [8]. By forming a new dynamic shape function matrix, dynamic analysis of the spatial frame structure is presented by the PSFEM [9].

It is significant to extend this research to the dynamic state. Considering the influence of random factors, the mechanical vibrations for a linear system are illustrated by using the TSFEM and the CG.

2. Random variable

Material properties, geometry parameters and applied loads of machines are assumed to be independent random variables, and are indicated as $$a_1,a_2,\cdots ,$$ $$a_i,\cdots ,a_{n_1}$$ . Their means are $${\mu }_1,{\mu }_2,\cdots ,{\mu }_i,\cdots ,{\mu }_{n_1}$$ , and their variances are.. $${\sigma }_i{}^2,\cdots ,{\sigma }_{n_1}{}^2$$ . When they are subject to normal distributions, the standard method used to simulate them is to take advantage of well-tested computer programs. When they are subject to unknown distributions, the sample of $$a_1,a_2,\cdots ,a_i,\cdots ,a_{n_1}$$ can be generated from the following method:

$$P\left\{\left|x-{\mu }^{\prime }\right|\ge \epsilon \right\}\le \frac{{\sigma }^2}{{\epsilon }^2}$$ E1

where, $$x$$ is a random variable, $${\mu }^{\prime }$$ is the mean, $$\sigma$$ is the standard deviation,and $$\epsilon$$ is an arbitrary positive number. Eq.1 is called the Chebyschev inequality.

The Chebyschev inequality can also be written as

$$P\left\{\left|x-{\mu }^{\prime }\right|\epsilon \right\}\ge 1-\frac{{\sigma }^2}{{\epsilon }^2}$$ E2

After substituting $$\epsilon =6{\sigma }_i$$ , $$x=a_i$$ , Eq.2 becomes

$$P\left\{\left|a_i-{\mu }_i\right|6{\sigma }_i\right\}\ge 0.9722$$ E3

where

$$a_{i}6{\sigma }_i+{\mu }_i$$ E4

or

$$a_i-6{\sigma }_i+{\mu }_i$$ E5

If it is assumed that $$z^{\prime }$$ is a random number within the open interval (0,1), then

$$a_i=6{\sigma }_i{z}^{\prime }+{\mu }_i$$ E6

or

$$a_i=-6{\sigma }_i{z}^{\prime }+{\mu }_i$$ E7

Large numbers of samples of random variables $$a_1,a_2,\cdots ,a_i,\cdots ,a_{n_1}$$ are produced from Eqs.6 and 7 so as to resolve the stochastic finite element problem through Monte Carlo stimulation.

3. Dynamic analysis of finite element

For a linear system, the dynamic equilibrium equation is given by

$$\left[M\right]\left\{\ddot{\delta }\right\}+\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}$$ E8

where $$\left\{\ddot{\delta }\right\},\left\{\dot{\delta }\right\},\left\{\delta \right\}$$ are the acceleration, velocity and displacement vectors. $$\left[M\right],\left[K\right]$$ and $$\left[C\right]$$ are the global mass, stiffness and damping matrices obtained by assembling the element variables in global coordinate system.

In order to program easily, the comprehensive calculation steps of the Newmark method are as follows

1. The initial calculation

The matrices $$\left[K\right]$$ , $$\left[M\right]$$ and $$\left[C\right]$$ are formed.

The initial values $$\left\{{\delta }_t\right\},\left\{{\dot{\delta }}_t\right\},\left\{{\ddot{\delta }}_t\right\}$$ are given.

After selecting step $$\Delta t$$ and parameters $$\gamma ,\beta$$ , the following relevant parameters are calculated:

$$\gamma \ge 0.50$$

$$\beta \ge 0.25{\left(0.5+\gamma \right)}^2$$ E9

$$b_0=\frac{1}{\beta {\left(\Delta t\right)}^2}$$ $$b_1=\frac{\gamma }{\beta \Delta t}$$

$$b_2=\frac{1}{\beta \Delta t}$$ E10

$$b_3=\frac{1}{2\beta }-1$$ $$b_4=\frac{\gamma }{\beta }-1$$

$$b_5=\frac{\Delta t}{2}\left(\frac{\gamma }{\beta }-2\right)$$ E11

$$b_6=\Delta t\left(1-\gamma \right)$$

$$b_7=\gamma \Delta t$$ E12

The stiffness matrix is defined as

$$\left[\tilde{K}\right]=\left[K\right]+b_0\left[M\right]+b_1\left[C\right]$$ E13

The stiffness matrix inversion $${\left[\tilde{K}\right]}^{-1}$$ is solved.

2. Calculation of each step time

At time $$t+\Delta t$$ , the load vector is defined as

$$\left\{{\tilde{F}}_{t+\Delta t}\right\}=\left\{F_{t+\Delta t}\right\}+\left[M\right](b_0\left\{{\delta }_t\right\}+$$

$$b_2\left\{{\dot{\delta }}_t\right\}+b_3\left\{{\ddot{\delta }}_t\right\})$$ E14

$$+\left[C\right]\left(b_1\left\{{\delta }_t\right\}+b_4\left\{{\dot{\delta }}_t\right\}+b_5\left\{{\ddot{\delta }}_t\right\}\right)$$ E15

At time $$t+\Delta t$$ , the displacement vector is given by

$$\left\{{\delta }_{t+\Delta t}\right\}={\left[\tilde{K}\right]}^{-1}\left\{{\tilde{F}}_{t+\Delta t}\right\}$$ E16

At time $$t+\Delta t$$ , the velocity vector and acceleration vector are obtained as

$$\left\{{\ddot{\delta }}_{t+\Delta t}\right\}=b_0\left(\left\{{\delta }_{t+\Delta t}\right\}-\left\{{\delta }_t\right\}\right)-b_2\left\{{\dot{\delta }}_t\right\}-b_3\left\{{\ddot{\delta }}_t\right\}$$ E17

$$\left\{{\dot{\delta }}_{t+\Delta t}\right\}=\left\{{\dot{\delta }}_t\right\}+b_6\left\{{\ddot{\delta }}_t\right\}+b_7\left\{{\ddot{\delta }}_{t+\Delta t}\right\}$$ E18

Vectors $$\left\{{\delta }_{t+i_1\Delta t}\right\},\left\{{\dot{\delta }}_{t+i_1\Delta t}\right\},\left\{{\ddot{\delta }}_{t+i_1\Delta t}\right\}$$ are solved at time $$t+i_1\Delta t\left(i_1=\mathrm{2,3,}\cdots ,n_3\right)$$ step-by -step.

4. Analysis of mechanical vibration based on CG

Eq.11 can be expressed as

$$\left[\tilde{K}\right]\left\{{\delta }_{t+\Delta t}\right\}=\left\{{\tilde{F}}_{t+\Delta t}\right\}$$ E19

$$N_1$$ samples of random variables $$a_1,a_2,\cdots ,a_i,\cdots ,a_{n_1}$$ are produced. $$N_1$$ matrices $$\left[\tilde{K}\right]$$ and $$N_1$$ Eq.14 are also generated. For a linear vibration, Eq.14 is the system of linear equations. The CG method is an effective method for solving the large system of linear equations according to the following steps:

1. First, select an approximate solution as the initial value

$${\delta }^{\left(0\right)}{}_{t+\Delta t}={\left({\delta }^{\left(0\right)}{}_{{\left(t+\Delta t\right)}_1},{\delta }^{\left(0\right)}{}_{{\left(t+\Delta t\right)}_2},\cdots ,{\delta }^{\left(0\right)}{}_{{\left(t+\Delta t\right)}_{N^{\prime }}}\right)}^T$$ E20

2. Calculate the first residual vector

$$r^{\left(0\right)}=\left\{{\tilde{F}}_{t+\Delta t}\right\}-\left[\tilde{K}\right]{\delta }^{\left(0\right)}{}_{t+\Delta t}$$ E21

and vector

$$p^{\left(0\right)}={\left[\tilde{K}\right]}^T{r}^{\left(0\right)}$$ E22

where, $${\left[\tilde{K}\right]}^T$$ is the transposed matrix of $$\left[\tilde{K}\right]$$

3. For $$\tilde{i}=\mathrm{0,1,2,}\cdots ,n_2-1$$ , iterate step-by-step as follows

$${\alpha }_{\tilde{i}}=\frac{\left(\left[\tilde{K}\right]p^{\left(\tilde{i}\right)},r^{\left(\tilde{i}\right)}\right)}{\left(\left[\tilde{K}\right]p^{\left(\tilde{i}\right)},\left[\tilde{K}\right]p^{\left(\tilde{i}\right)}\right)}=\frac{\left(p^{\left(\tilde{i}\right)},{\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}\right)}\right)}{\left(\left[\tilde{K}\right]p^{\left(\tilde{i}\right)},\left[\tilde{K}\right]p^{\left(\tilde{i}\right)}\right)}=\frac{\left({\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}\right)},{\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}\right)}\right)}{\left(\left[\tilde{K}\right]p^{\left(\tilde{i}\right)},\left[\tilde{K}\right]p^{\left(\tilde{i}\right)}\right)}$$ E23

$${\left\{{\delta }_{t+\Delta t}\right\}}^{\left(\tilde{i}+1\right)}={\left\{{\delta }_{t+\Delta t}\right\}}^{\left(\tilde{i}\right)}+{\alpha }_{\tilde{i}}{p}^{\left(\tilde{i}\right)}$$ E24

$$r^{\left(\tilde{i}+1\right)}=r^{\left(\tilde{i}\right)}-{\alpha }_{\tilde{i}}\left[\tilde{K}\right]p^{\left(\tilde{i}\right)}$$ E25

$${\beta }_{\tilde{i}+1}=\frac{\left({\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}+1\right)},{\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}+1\right)}\right)}{\left({\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}\right)},{\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}\right)}\right)}$$ E26

$$p^{\left(\tilde{i}+1\right)}={\left[\tilde{K}\right]}^T{r}^{\left(\tilde{i}+1\right)}+{\beta }_{\tilde{i}+1}{p}^{\left(\tilde{i}\right)}$$ E27

The process can be stopped only if $$r^{n_2}$$ is small enough.

Vectors $${\left\{{\delta }_{t+\Delta t}\right\}}_1,{\left\{{\delta }_{t+\Delta t}\right\}}_2,\cdots ,{\left\{{\delta }_{t+\Delta t}\right\}}_{N_1}$$ are solutions of $$N_1$$ Eqs.14.

The mean of $$\left\{{\delta }_{t+\Delta t}\right\}$$ is given by

$$\mu \left\{{\delta }_{t+\Delta t}\right\}=\frac{{\left\{{\delta }_{t+\Delta t}\right\}}_1+{\left\{{\delta }_{t+\Delta t}\right\}}_2+\cdots +{\left\{{\delta }_{t+\Delta t}\right\}}_{N_1}}{N_1}$$ E28

The variance of $$\left\{{\delta }_{t+\Delta t}\right\}$$ is given by

$$Var\left\{{\delta }_{t+\Delta t}\right\}=\frac{1}{N_1-1}{\displaystyle \sum _{i=1}^{N_1}{\left({\left\{{\delta }_{t+\Delta t}\right\}}_i-\mu \left\{{\delta }_{t+\Delta t}\right\}\right)}^2}$$ E29

Similarly, the mean and variance of the vector $$\left\{{\delta }_{t+i_1\Delta t}\right\}$$ can be solved for at time $$t+i_1\Delta t\left(i_1=\mathrm{2,3,}\cdots ,n_3\right)$$ step-by-step.

At time $$t^{\prime }=t+i_2\Delta t\left(i_2=\mathrm{1,2,}\cdots ,n_4\right)$$ ,the strain and stress vectors for element $$d$$ are

$$\left\{\epsilon \right\}=\left[B\right]\left\{{\delta }_{t^{\prime }}{}^d\right\}$$ E30

and

$$\left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}$$ E31

where, $$\left[D\right]$$ =the material response matrix of element $$d$$ , $$\left[B\right]$$ =the gradient matrix of element $$d$$ and $$\left\{{\delta }_{t^{\prime }}{}^d\right\}$$ =the element $$d$$ nodal displacement vector at time $$t^{\prime }$$ .

Substituting Eq.25 into Eq.26, the stress for element $$d$$ is given by

$$\left\{\sigma \right\}=\left[D\right]\left[B\right]\left\{{\delta }_{t^{\prime }}{}^d\right\}$$ E32

Substituting..samples of random variables $$a_1,a_2,\cdots ,a_i,\cdots ,a_{n_1}$$ into Eq.27, the vectors $${\left\{\sigma \right\}}_1,{\left\{\sigma \right\}}_2,\cdots ,{\left\{\sigma \right\}}_{N_1}$$ can be obtained.

The mean of $$\left\{\sigma \right\}$$ is given by

$$\mu \left\{\sigma \right\}=\frac{{\left\{\sigma \right\}}_1+{\left\{\sigma \right\}}_2+\cdots +{\left\{\sigma \right\}}_{N_1}}{N_1}$$ E33

The variance of $$\left\{\sigma \right\}$$ is given by

$$Var\left\{\sigma \right\}=\frac{1}{N_1-1}{{\displaystyle \sum _{i=1}^{N_1}\left({\left\{\sigma \right\}}_i-\mu \left\{\sigma \right\}\right)}}^2$$ E34

The CG method belongs to method of iteration with the advantage of quick convergence. For practical purpose, PCG is applied to accelerate the convergence.

5. Analysis of mechanical vibration based on TSFEM

Independent random variables of the system are regarded as $$a_1,a_2,\cdots ,a_i,\cdots ,a_{n_1}$$ .

The partial derivative of Eq.14 with respect to $$a_i$$ is given by

$$\frac{\partial \left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i}={\left[\tilde{K}\right]}^{-1}\left(\frac{\partial \left\{{\tilde{F}}_{t+\Delta t}\right\}}{\partial a_i}-\frac{\partial \left[\tilde{K}\right]}{\partial a_i}\left\{{\delta }_{t+\Delta t}\right\}\right)$$ E35

where

$$\frac{\partial \left\{{\tilde{F}}_{t+\Delta t}\right\}}{\partial a_i}=\frac{\partial \left\{F_{t+\Delta t}\right\}}{\partial a_i}+\left[M\right]\left(b_0\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}+b_2\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}+b_3\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}\right)+$$ E36

$$\frac{\partial \left[M\right]}{\partial a_i}\left(b_0\left\{{\delta }_t\right\}+b_2\left\{{\dot{\delta }}_t\right\}+b_3\left\{{\ddot{\delta }}_t\right\}\right)$$ E37

$$\begin{array}{l}+\left[C\right]\left(b_1\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}+b_4\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}+b_5\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}\right)\\ +\frac{\partial \left[C\right]}{\partial a_i}\left(b_1\left\{{\delta }_t\right\}+b_4\left\{{\dot{\delta }}_t\right\}+b_5\left\{{\ddot{\delta }}_t\right\}\right)\end{array}$$ E38

After $$\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}|{}_{t=0}=q_0,\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}|{}_{t=0}={\dot{q}}_0,\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}|{}_{t=0}={\ddot{q}}_0$$ are given, Eq.31 can be calculated.

The partial derivative of Eq.30 with respect to $$a_i$$ is given by

$$\frac{{\partial }^2\left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i{}^2}={\left[\tilde{K}\right]}^{-1}$$

$$\left(\frac{{\partial }^2\left\{{\tilde{F}}_{t+\Delta t}\right\}}{\partial a_i{}^2}-2\frac{\partial \left[\tilde{K}\right]}{\partial a_i}\frac{\partial \left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i}-\frac{{\partial }^2\left[\tilde{K}\right]}{\partial a_i{}^2}\left\{{\delta }_{t+\Delta t}\right\}\right)$$ E39

where

$$\frac{{\partial }^2\left\{{\tilde{F}}_{t+\Delta t}\right\}}{\partial a_i{}^2}=\frac{{\partial }^2\left\{F_{t+\Delta t}\right\}}{\partial a_i{}^2}+\frac{{\partial }^2\left[M\right]}{\partial a_i{}^2}\left(b_0\left\{{\delta }_t\right\}+b_2\left\{{\dot{\delta }}_t\right\}+b_3\left\{{\ddot{\delta }}_t\right\}\right)+$$ E40

$$2\frac{\partial \left[M\right]}{\partial a_i}\left(b_0\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}+b_2\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}+b_3\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}\right)$$

$$+\left[M\right]\left(b_0\frac{{\partial }^2\left\{{\delta }_t\right\}}{\partial a_i{}^2}+b_2\frac{{\partial }^2\left\{{\dot{\delta }}_t\right\}}{\partial a_i{}^2}+b_3\frac{{\partial }^2\left\{{\ddot{\delta }}_t\right\}}{\partial a_i{}^2}\right)$$ E41

$$\begin{array}{l}+\frac{{\partial }^2\left[C\right]}{\partial a_i{}^2}\left(b_1\left\{{\delta }_t\right\}+b_4\left\{{\dot{\delta }}_t\right\}+b_5\left\{{\ddot{\delta }}_t\right\}\right)+2\frac{\partial \left[C\right]}{\partial a_i}\left(b_1\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}+b_4\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}+b_5\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}\right)\\ +\left[C\right]\left(b_1\frac{{\partial }^2\left\{{\delta }_t\right\}}{\partial a_i{}^2}+b_4\frac{{\partial }^2\left\{{\dot{\delta }}_t\right\}}{\partial a_i{}^2}+b_5\frac{{\partial }^2\left\{{\ddot{\delta }}_t\right\}}{\partial a_i{}^2}\right)\end{array}$$ E42

After $$\frac{{\partial }^2\left\{{\delta }_t\right\}}{\partial a_i{}^2}|{}_{t=0}=r_0,\frac{{\partial }^2\left\{{\dot{\delta }}_t\right\}}{\partial a_i{}^2}|{}_{t=0}={\dot{r}}_0,\frac{{\partial }^2\left\{{\ddot{\delta }}_t\right\}}{\partial a_i{}^2}|{}_{t=0}={\ddot{r}}_0$$ are given, Eq.33 can be calculated.

The displacement is expanded at the mean value point $$\overline{a}={\left({\overline{a}}_1,{\overline{a}}_2,\cdots ,{\overline{a}}_i,\cdots ,{\overline{a}}_{n_1}\right)}^T$$ by means of a Taylor series. By taking the expectation operator for two sides of above Eq.11, the mean of the displacement is obtained as

$$\mu \left\{{\delta }_{t+\Delta t}\right\}\approx {\left\{{\delta }_{t+\Delta t}\right\}|}_{a=\overline{a}}+\frac{1}{2}{{\displaystyle \sum _{i=1}^{n_1}\frac{{\partial }^2\left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i{}^2}|}}_{a=\overline{a}}\cdot {\sigma }_i{}^2$$ E43

where, $$\mu \left\{{\delta }_{t+\Delta t}\right\}$$ expresses the mean value of $${\delta }_{t+\Delta t}$$ .

The variance of $${\delta }_{t+\Delta t}$$ is given by

$$Var\left\{{\delta }_{t+\Delta t}\right\}\approx {\displaystyle \sum _{i=1}^{n_1}{\left({\frac{\partial \left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i}|}_{a=\overline{a}}\right)}^2\cdot {\sigma }_i{}^2}$$ E44

The partial derivative of $${\ddot{\delta }}_{t+\Delta t}$$ with respect to $$a_i$$ is given by

$$\frac{\partial \left\{{\ddot{\delta }}_{t+\Delta t}\right\}}{\partial a_i}=b_0\left(\frac{\partial \left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i}-\frac{\partial \left\{{\delta }_t\right\}}{\partial a_i}\right)-b_2\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}-b_3\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}$$ E45

The partial derivative of $${\dot{\delta }}_{t+\Delta t}$$ with respect to $$a_i$$ is given by

$$\frac{\partial \left\{{\dot{\delta }}_{t+\Delta t}\right\}}{\partial a_i}=\frac{\partial \left\{{\dot{\delta }}_t\right\}}{\partial a_i}+b_6\frac{\partial \left\{{\ddot{\delta }}_t\right\}}{\partial a_i}+b_7\frac{\partial \left\{{\ddot{\delta }}_{t+\Delta t}\right\}}{\partial a_i}$$ E46

The partial derivative of Eq.36 with respect to $$a_i$$ is given by

$$\frac{{\partial }^2\left\{{\ddot{\delta }}_{t+\Delta t}\right\}}{\partial a_i{}^2}=b_0\left(\frac{{\partial }^2\left\{{\delta }_{t+\Delta t}\right\}}{\partial a_i{}^2}-\frac{{\partial }^2\left\{{\delta }_t\right\}}{\partial a_i{}^2}\right)$$

$$-b_2\frac{{\partial }^2\left\{{\dot{\delta }}_t\right\}}{\partial a_i{}^2}-b_3\frac{{\partial }^2\left\{{\ddot{\delta }}_t\right\}}{\partial a_i{}^2}$$ E47

The partial derivative of Eq.37 with respect to $$a_i$$ is given by

$$\frac{{\partial }^2\left\{{\dot{\delta }}_{t+\Delta t}\right\}}{\partial a_i{}^2}=\frac{{\partial }^2\left\{{\dot{\delta }}_t\right\}}{\partial a_i{}^2}+$$

$$b_6\frac{{\partial }^2\left\{{\ddot{\delta }}_t\right\}}{\partial a_i{}^2}+b_7\frac{{\partial }^2\left\{{\ddot{\delta }}_{t+\Delta t}\right\}}{\partial a_i{}^2}$$ E48

The mean value and variance of the displacement are obtained at time $$t+i_1\Delta t\left(i_1=\mathrm{2,3,}\cdots ,n_3\right)$$ step-by-step.

The partial derivative of Eq.27 with respect to $$a_i$$ is given by

$$\frac{\partial \left\{\sigma \right\}}{\partial a_i}=\frac{\partial \left[D\right]}{\partial a_i}\left[B\right]\left\{{\delta }_{t^{\prime }}{}^d\right\}+\left[D\right]\frac{\partial \left[B\right]}{\partial a_i}\left\{{\delta }_{t^{\prime }}{}^d\right\}+\left[D\right]\left[B\right]\frac{\partial \left\{{\delta }_{t^{\prime }}{}^d\right\}}{\partial a_i}$$ E49

The partial derivative of Eq.40 with respect to $$a_i$$ is given by

$$\frac{{\partial }^2\left\{\sigma \right\}}{\partial a_i{}^2}=\frac{{\partial }^2\left[D\right]}{\partial a_i{}^2}\left[B\right]\left\{{\delta }_{t^{\prime }}{}^d\right\}+$$

$$2\frac{\partial \left[D\right]}{\partial a_i}\frac{\partial \left[B\right]}{\partial a_i}\left\{{\delta }_{t^{\prime }}{}^d\right\}+2\frac{\partial \left[D\right]}{\partial a_i}\left[B\right]\frac{\partial \left\{{\delta }_{t^{\prime }}{}^d\right\}}{\partial a_i}$$ E50

$$+\left[D\right]\frac{{\partial }^2\left[B\right]}{\partial a_i{}^2}\left\{{\delta }_{t^{\prime }}{}^d\right\}+2\left[D\right]\frac{\partial \left[B\right]}{\partial a_i}\frac{\partial \left\{{\delta }_{t^{\prime }}{}^d\right\}}{\partial a_i}$$ E51

$$+\left[D\right]\left[B\right]\frac{{\partial }^2\left\{{\delta }_{t^{\prime }}{}^d\right\}}{\partial a_i{}^2}$$ E52

The stress is expanded at mean value point $$\overline{a}={\left({\overline{a}}_1,{\overline{a}}_2,\cdots ,{\overline{a}}_i,\cdots ,{\overline{a}}_{n_1}\right)}^T$$ by means of a Taylor series. By taking the expectation operator for two sides of the above Eq.27, the mean of stress is obtained as

$$\mu \left\{\sigma \right\}\approx {\left\{\sigma \right\}|}_{a=\overline{a}}+\frac{1}{2}{{\displaystyle \sum _{i=1}^{n_1}\frac{{\partial }^2\left\{\sigma \right\}}{\partial a_i{}^2}|}}_{a=\overline{a}}\cdot {\sigma }_i{}^2$$ E53

where, $$\mu \left\{\sigma \right\}$$ expresses the mean value of $$\sigma$$ .

The variance of $$\sigma$$ is given by

$$Var\left\{\sigma \right\}\approx {\displaystyle \sum _{i=1}^{n_1}{\left({\frac{\partial \left\{\sigma \right\}}{\partial a_i}|}_{a=\overline{a}}\right)}^2\cdot {\sigma }_i{}^2}$$ E54

6. Numerical example

Figure 1 shows a four-bar linkage, or a crank and rocker mechanism. The establishment of differential equation system can be found in literature 10，11，12.The length of bar 1 is 0.075m, the length of bar 2 is 0.176m, the length of bar 3 is 0.29m,and the length of the bar 4 is0. 286m, the diameters of three bars are 0.02m. The torque T is 4Nm, the load F1 is 20sint N. The three bars are made of steel and they are regarded as three elements. Considering the boundary condition, there are 13 unit coordinates. Young’s modulus is regarded as a random variable. For numerical calculation, the means of the Young’s modulus within the three bars are. $$2\times {10}^{11}$$ . $$N/m^2$$ and the variances of the Young’s modulus are $${10}^{11}$$ $$N^2/m^4$$ .

Figure 2 shows the mean of the displacement at unit coordinate 11. Unit coordinate 11 is the deformation of the upper end of bar 3 in the vertical direction. The DSFEM simulates 1000 samples. The TSFEM produces an error of less than 0.5%. The CG produces an error of less than 0.1%. Figure 3 shows the variance of the displacement at unit coordinate 11. TSFEM produces an error of less than 1.0%. CG produces an error of less than 0.4%.Figure 4 shows the mean of stress at the top of bar 3. The TSFEM produces an error of less than 0.85%.The CG produces an error of less than 0.13%.Figure 5 shows the variance of stress at the top of bar 3. The TSFEM produces an error of less than 1%. The CG produces an error of less than 0.3%.The results obtained by the CG method and the TSFEM are very close to that obtained by the DSFEM. Table 1 indicates the comparison of CPU time when the mechanism has operated for six seconds.

Figure 6 shows a cantilever beam. The length, the width, the height, the Poisson’s ratio,the Young’s modulus and the load F are assumed to be random variables. Their means are 1m, 0.1m, 0.05m, 0.2, $$2\times {10}^{11}$$ $$N/m^2$$ ,100N.Their standard deviation are 0.2, 0.1, 0.1, 0.01, $${10}^9$$ , 0.1. Load subjected to the cantilever beam is Fsin(100t)N. It is divided into 400 rectangle elements that have 505 nodes. Figure 7 shows the mean of vertical displacement at node 505. DSFEM simulates 100 samples. The result obtained by the TSFEM produces an error of less than 2%. CG produces an error of less than 0.5%. Figure 8 shows the variance of vertical displacement at node 505.The TSFEM produces an error of less than 3.0%. CG produces an error of less than 0.8%.Figure 9 shows the mean of horizontal stress at node 5. The TSFEM produces an error of less than 2.4%. CG produces an error of less than 0.9%. Figure 10 shows the variance of horizontal stress at node 5. The TSFEM produces an error of less than 3.2%. CG produces an error of less than 1.3%. Table 2 indicates the comparison of CPU time when the cantilever beam has operated for six seconds.

Figure 1.

A four-bar linkage

Figure 2.

The mean of displacement at unit coordinate 11 for $${\sigma }_E{}^2={10}^{11}$$

Figure 3.

The variance of displacement at unit coordinate 11 for $${\sigma }_E{}^2={10}^{11}$$

Figure 4.

The mean of stress at the top of bar 3 for $${\sigma }_E{}^2={10}^{11}$$

Figure 5.

The variance of stress at the top of bar 3 for $${\sigma }_E{}^2={10}^{11}$$

DSFEM	TSFEM	CG
19 seconds	4 seconds	14 seconds

Table 1.

Comparison of CPU time for $${\sigma }_E{}^2={10}^{11}$$

Figure 6.

A cantilever beam

Figure 7.

The mean of vertical displacement at node 505

Figure 8.

The variance of vertical displacement at node 505

Figure 9.

The mean of horizontal stress at node 5

Figure 10.

The variance of horizontal stress at node 5

DSFEM	TSFEM	CG
3 hours 8 minutes 17 seconds	1 hour 45 minutes 10 seconds	40 minutes 24 seconds

Table 2.

Comparison of CPU time

7. Conclusions

Considering the influence of random factors, the mechanical vibration in a linear system is presented by using the TSFEM. Different samples of random variables are simulated. The combination of CG method and Monte Carlo method makes it become an effective method for analyzing the vibration problem with the characteristics of high accuracy and quick convergence.