Stochastic Finite Element for Structural Vibration

Abstract

This paper proposes a new method of calculating stochastic field. It is an improvement of the midpoint method of stochastic field. The vibration equation of a system is transformed to a static problem by using the Newmark method and the Taylor expansion is extended for the structural vibration analysis with uncertain factors. In order to develop computational efficiency and allow for efficient storage, the Conjugate Gradient method (CG) is also employed. An example is given, respectively, and calculated results are compared to validate the proposed methods.

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, regarding the given factors as constants cannot apparently correspond to the reality of a structure. In order to enhance computational accuracy, the influence of random factors must be considered.

Many physics parameters of material possess spatial variability, such as Young’s modulus and Poisson’s ratio, so we should regard them as stochastic fields. Stochastic field discretization is the problem that various stochastic finite element methods need to resolve, but discrete form of stochastic field plays the decisive influence on the calculation and computational accuracy of stochastic finite element. The simplest discretization is the midpoint method (MSF) [1]. The stochastic field is described by a single random variable replacing the value of the field at a central point of the mesh. The local average method of the stochastic field describes the stochastic field of an element in terms of the spatial average. The local average method of rectangle element is described by the mean, variance, and covariance [2]. It can be extended for 3D [3]. The stochastic field of nonrectangular element is described by the mean vector and covariance matrix using Gaussian quadrature [4]. The stochastic field can be described by the shape function and nodal values, and it is necessary to know the related function [5]. Making use of the Karhunen-Loeve expansion, stochastic field is represented by series expansion [6]. When stiffness matrices are deduced, a weighted integral method is adopted to consider stochastic field of material parameters [7, 8]. When stochastic field is expressed by series expansion, the optimal linear-estimation method is applied to make the error of variance minimum [9].

The direct Monte Carlo simulation of the stochastic finite element method (DSFEM) requires a large number of samples, which requires much calculation time [10]. 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 [11]. The preconditioned Conjugate Gradient method (PCG) applied in the calculation of stochastic finite elements can also enhance computational accuracy and efficiency [12]. According to first-order or second-order perturbation methods, calculation formulas can be obtained [5, 11, 14–16, 18–20, 22]. The result is called the PSFEM and has been adopted by many authors.

The PSFEM is often applied in dynamic analysis of structures and the second-order perturbation technique has been proved to be efficient [5]. Dynamic reliability of a frame is calculated by the SFEM and response sensitivity is formulated in the context of stiffness and mass matrix condensation [16]. When load actions are treated as stochastic processes, vibration of the structure is resolved by the PSFEM [17]. The PSFEM is an adequate tool for nonlinear structural dynamics. Nonlinearities due to material and geometrical effects have also been included [18]. By forming a new dynamic shape function matrix, dynamic analysis of the spatial frame structure is presented by the PSFEM [19]. The NSFEM is introduced in dynamic analysis within the framework of a Monte Carlo simulation [20]. The NSFEM is applied to the dynamic response of a random structure system and results are compared with those from the PSFEM and the DSFEM [21]. With the aid of the PSFEM, a stochastic formulation for nonlinear dynamic analysis of a structure is presented [22].

An improved midpoint method of stochastic field (IMSF) is presented. The IMSF is more accurate than the MSF. The Newmark method transforms differential equations into linear equations. The IMSF is used and the structural vibrations for a linear system are computed by the Taylor expansion method, the CG method, and the PCG method of stochastic finite element (TSFEM, CG, PCG). The TSFEM, the CG, and the PCG based on the MSF are called the MTSFEM, the MCG, and the MPCG. An example demonstrates the superiority of the proposed methods.

2. Improved Midpoint Method of Stochastic Field

When finite element method is used, structure is divided into small elements whose number is appropriate. In this paper, Young’s modulus is assumed to be a Gaussian process. When element is appropriately small, Young’s modulus of an element is described by a variable. The Young’s modulus of structure is described by a group of variables. Without loss of generality, it is supposed that the structure is divided into elements of 𝑚 nodes and 𝑛 nodes. The Young’s modulus of node within element 𝑒 of 𝑚 nodes is expressed by 𝑎_𝑒𝑚1,𝑎_𝑒𝑚2,…,𝑎_𝑒𝑚𝑚. The Young’s modulus of midpoint within element 𝑒 of 𝑚 nodes is expressed by 𝑎_𝑒𝑚𝑙. The Young’s modulus of node within element 𝑓 of 𝑛 nodes is expressed by 𝑎_𝑓𝑛1,𝑎_𝑓𝑛2,…,𝑎_𝑓𝑛𝑛. The Young’s modulus of midpoint within element 𝑓 of n nodes is expressed by 𝑎_𝑓𝑛𝑙.

The Young’s modulus within element 𝑒 is defined as𝑎_𝑒=𝑎_𝑒𝑚1+𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙𝑚+1. (2.1)

Its mean is𝜇_𝑒𝑎=𝜇(^𝑒𝑚1+𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙)=𝜇𝑚+1_{𝑎𝑒𝑚1}+𝜇_{𝑎𝑒𝑚2}+⋯+𝜇_{𝑎𝑒𝑚𝑚}+𝜇_{𝑎𝑒𝑚𝑙},𝑚+1 (2.2) where the means of Young’s modulus at the first node, the second node,…, the 𝑚th node of element 𝑒 are expressed by 𝜇_{𝑎𝑒𝑚1},𝜇_{𝑎𝑒𝑚2},…,𝜇_{𝑎𝑒𝑚𝑚}.The mean of Young’s modulus at the midpoint of element 𝑒 is expressed by 𝜇_{𝑎𝑒𝑚𝑙}.

The Young’s modulus within element 𝑓 is defined as𝑎_𝑓=𝑎_𝑓𝑛1+𝑎_𝑓𝑛2+⋯+𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙𝑛+1. (2.3)

Its mean is 𝜇_𝑓𝑎=𝜇(^𝑓𝑛1+𝑎_𝑓𝑛2+⋯+𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙)=𝜇𝑛+1_{𝑎𝑓𝑛1}+𝜇_{𝑎𝑓𝑛2}+⋯+𝜇_{𝑎𝑓𝑛𝑛}+𝜇_{𝑎𝑓𝑛𝑙},𝑛+1 (2.4) where the means of Young’s modulus at the first node, the second node,…,the 𝑛th node of element 𝑓 are expressed by 𝜇_{𝑎𝑓𝑛1},𝜇_{𝑎𝑓𝑛2},…,𝜇_{𝑎𝑓𝑛𝑛}. The mean of Young’s modulus at the midpoint of element 𝑓 is expressed by 𝜇_{𝑎𝑓𝑛𝑙}.

The covariance of Young’s modulus between two elements that each has 𝑚 nodes is obtained byCov(𝑎_𝑒,𝑎_𝑒′𝑎)=Cov(^𝑒𝑚1+𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎𝑚+1^𝑒′𝑚1+𝑎_𝑒′𝑚2+⋯+𝑎_{𝑒′𝑚𝑚}+𝑎_{𝑒′𝑚𝑙})=1𝑚+1(𝑚+1)²[Cov(𝑎_𝑒𝑚1,𝑎_𝑒′𝑚1+𝑎_𝑒′𝑚2+⋯+𝑎_{𝑒′𝑚𝑚}+𝑎_{𝑒′𝑚𝑙})+Cov(𝑎_𝑒𝑚2+⋯𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎_𝑒′𝑚1+𝑎_𝑒′𝑚2+⋯+𝑎_{𝑒′𝑚𝑚}+𝑎_{𝑒′𝑚𝑙}=1)](𝑚+1)²[Cov(𝑎_𝑒𝑚1,𝑎_𝑒′𝑚1)+Cov(𝑎_𝑒𝑚1,𝑎_𝑒′𝑚2)+⋯+Cov(𝑎_𝑒𝑚1,𝑎_{𝑒′𝑚𝑚})+Cov(𝑎_𝑒𝑚1,𝑎_{𝑒′𝑚𝑙}+1)](𝑚+1)²Cov(𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎_𝑒′𝑚1+𝑎_𝑒′𝑚2+⋯+𝑎_{𝑒′𝑚𝑚}+𝑎_{𝑒′𝑚𝑙})=1(𝑚+1)²(^𝑚_𝑔1^𝑚_=1𝑔2=1Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑒′𝑚𝑔2}1))+(𝑚+1)²(^𝑚_𝑔1=1Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑒′𝑚𝑙}+1))(𝑚+1)²(^𝑚_𝑔2=1Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑒′𝑚𝑔2})+Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑒′𝑚𝑙})), (2.5) where Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑒′𝑚𝑔2}) = the covariance of Young’s modulus between node 𝑔₁(𝑔₁=1,2,…,𝑚) of element 𝑒 and node 𝑔₂(𝑔₂=1,2,…,𝑚) of element 𝑒^′, Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑒′𝑚𝑙}) = the covariance of Young’s modulus between node 𝑔₁ of element 𝑒 and the midpoint of element 𝑒^′, Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑒′𝑚𝑔2})  = the covariance of Young’s modulus between the midpoint of element 𝑒 and node 𝑔₂ of element 𝑒^′, and Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑒′𝑚𝑙}) = the covariance of Young’s modulus between the midpoint of element 𝑒 and the midpoint of element 𝑒^′.

The covariance of Young’s modulus between two elements is given in Appendix A.

Using covariance matrix, the correlation of Young’s modulus between any two elements is given by𝐶_𝑎𝑎=Cov(𝑎₁,𝑎₁)Cov(𝑎₁,𝑎₂)⋯Cov(𝑎₁,𝑎_𝑁)Cov(𝑎₂,𝑎₁)Cov(𝑎₂,𝑎₂)⋯Cov(𝑎₂,𝑎_𝑁)⋮⋮⋮Cov(𝑎_𝑁,𝑎₁)Cov(𝑎_𝑁,𝑎₂)⋯Cov(𝑎_𝑁,𝑎_𝑁). (2.6)

A Gaussian vector ^⇀𝑎=[𝑎₁,𝑎₂,…,𝑎_𝑁]^𝑇 is generated^⇀𝑎=𝐿𝑍. (2.7)𝑍=[𝑍₁,𝑍₂,…,𝑍_𝑁]^𝑇 consists of 𝑁 Gaussian random variables with mean zero and unit standard deviation. The Cholesky matrix 𝐿 can be obtained through a decomposition of the covariance matrix; therefore,𝜇[𝑍𝑍^𝑇]=𝐼,𝐿𝐿^𝑇=𝐶_𝑎𝑎. (2.8)𝐼 is the identity matrix. The generation of vector ^⇀𝑎 must satisfy the covariance matrix𝜇[^⇀𝑎^⇀𝑎^𝑇]=𝜇[𝐿𝑍(𝐿𝑍)^𝑇]=𝐿𝜇[𝑍𝑍^𝑇]𝐿^𝑇=𝐶_𝑎𝑎. (2.9)

Once the decomposition has been completed, different samples of vector ^⇀𝑎 can be acquired easily by (2.7). Thus, it is possible that Monte Carlo simulation resolves problem of stochastic finite element.

3. Dynamic Analysis of Finite Element

For a linear system, the dynamic equilibrium equation is given by[𝑀]{𝜹̈}+[𝐶]{𝜹̇}+[𝐾]{𝜹}={𝐹}, (3.1) where {𝛿̈},{𝛿̇},{𝛿} are the acceleration, velocity, and displacement vectors. [𝑀],[𝐾], and [𝐶] are the global mass, stiffness, and damping matrices obtained by assembling the element variables in global coordinate system.

By using the Newmark method, (3.1) becomes{𝛿_𝑡+Δ𝑡}=[𝐾]⁻¹{𝐹_𝑡+Δ𝑡}, (3.2) where {𝛿_𝑡+Δ𝑡}，[𝐾] and {𝐹_𝑡+Δ𝑡} indicate the displacement vector, stiffness matrix and load vector at time 𝑡+Δ𝑡. [𝐾] and {𝐹_𝑡+Δ𝑡} are given in Appendix B.

4. Dynamic Analysis of Structure Based on CG

Equation (3.2) can be rewritten as[𝐾]{𝛿_𝑡+▵𝑡}={𝐹_𝑡+▵𝑡}. (4.1) Using (2.7), (2.8), and (2.9), 𝑁₁ samples of vector ^⇀𝑎 are produced. 𝑁₁ matrices [𝐾] and 𝑁₁ (4.1) are generated. For linear vibrations, (4.1) is a system of linear equations. The CG is an adequate method for solving large systems of linear equations. It can be accomplished as follows.

first, select appropriate solution as initial values {𝛿_𝑡+Δ𝑡}⁽⁰⁾=(𝛿⁽⁰⁾_{(𝑡+Δ𝑡)1},𝛿⁽⁰⁾_{(𝑡+Δ𝑡)2},…,𝛿⁽⁰⁾_{(𝑡+Δ𝑡)𝑁1})^𝑇; (4.2)

calculate the first residual vector 𝑟⁽⁰⁾={𝐹_𝑡+Δ𝑡}−[𝐾]{𝛿_𝑡+Δ𝑡}⁽⁰⁾ (4.3) and vector 𝑝⁽⁰⁾=[𝐾]^𝑇𝑟⁽⁰⁾ (4.4) where [𝐾]^𝑇is the transposed matrix;

for 𝑖̃=0,1,2,…,𝑛₂−1, iterate step by step as follows: 𝛼_̃𝑖=([𝐾]𝑝⁽𝑖₎̃,𝑟⁽𝑖₎̃)([𝐾]𝑝⁽𝑖₎̃,[𝐾]𝑝⁽𝑖₎̃)=(𝑝⁽𝑖₎̃,[𝐾]^𝑇𝑟⁽𝑖₎̃)([𝐾]𝑝⁽𝑖₎̃,[𝐾]𝑝⁽𝑖₎̃)=([𝐾]^𝑇𝑟⁽𝑖₎̃,[𝐾]^𝑇𝑟⁽𝑖₎̃)([𝐾]𝑝⁽𝑖₎̃,[𝐾]𝑝⁽𝑖₎̃),{𝛿_𝑡+▵𝑡}⁽𝑖̃₊₁₎={𝛿_𝑡+▵𝑡}⁽𝑖₎̃+𝛼_̃𝑖𝑝⁽𝑖₎̃,𝑟⁽𝑖̃₊₁₎=𝑟⁽𝑖₎̃−𝛼_̃𝑖[𝐾]𝑝⁽𝑖₎̃,𝛽_̃𝑖+1=([𝐾]^𝑇𝑟⁽𝑖̃₊₁₎,[𝐾]^𝑇𝑟⁽𝑖̃₊₁₎)([𝐾]^𝑇𝑟⁽𝑖₎̃,[𝐾]^𝑇𝑟⁽𝑖₎̃),𝑝⁽𝑖̃₊₁₎=[𝐾]^𝑇𝑟⁽𝑖̃₊₁₎+𝛽_̃𝑖+1𝑝⁽𝑖₎̃. (4.5)

This process stops only if 𝑟^𝑛2 is small enough.

Vectors {𝛿_𝑡+Δ𝑡}₁,{𝛿_𝑡+Δ𝑡}₂,…,{𝛿_𝑡+Δ𝑡}_𝑁1are solutions of 𝑁₁ (4.1).

The mean of {𝛿_𝑡+Δ𝑡} is given by𝜇{𝛿_𝑡+Δ𝑡}={𝛿_𝑡+Δ𝑡}₁+{𝛿_𝑡+Δ𝑡}₂+⋯+{𝛿_𝑡+Δ𝑡}_𝑁1𝑁₁. (4.6) The variance of {𝛿_𝑡+Δ𝑡} is given byVar{𝛿_𝑡+Δ𝑡1}=𝑁₁−1^𝑁1_𝑖=1({𝛿_𝑡+Δ𝑡}_𝑖−𝜇{𝛿_𝑡+Δ𝑡})². (4.7) Similarly, the mean and variance of vector {𝛿_{𝑡+𝑖1Δ𝑡}} can be solved at time 𝑡+𝑖₁Δ𝑡 step by step.

At time 𝑡^′=𝑡+𝑖₂Δ𝑡(𝑖₂=1,2,…,𝑛₁), the stress for element 𝑑 is given by{𝜎}=[𝐷][𝐵]{𝛿^𝑑_𝑡′}, (4.8) where [𝐷] = the material response matrix of element 𝑑, [𝐵] = the gradient matrix of element 𝑑 and {𝛿_𝑡′^𝑑} = the element 𝑑 nodal displacement vector at time 𝑡^′.

Substituting 𝑁₁ samples of vector ^⇀𝑎 into (4.8), vectors {𝜎}₁,{𝜎}₂,…,{𝜎}_𝑁1can be obtained.

The mean of {𝜎} is given by𝜇{𝜎}={𝜎}₁+{𝜎}₂+⋯+{𝜎}_𝑁1𝑁₁. (4.9) The variance of {𝜎} is given by1Var{𝜎}=𝑁₁−1^𝑁1_𝑖=1({𝜎}_𝑖−𝜇{𝜎})². (4.10) The CG belongs to methods of iteration that converge quickly. For practical purposes, PCG is applied to accelerate convergence.

5. Dynamic Analysis of Structure Based on the TSFEM

Young’s modulus of the structure is given as function of 𝑁 random variables 𝑎₁,𝑎₂,…,𝑎_𝑁. The partial derivative of (4.1) with respect to 𝑎_𝑖 is given by𝜕{𝛿_𝑡+▵𝑡}𝜕𝑎_𝑖=[𝐾]⁻¹(𝜕{𝐹_𝑡+▵𝑡}𝜕𝑎_𝑖−𝜕[𝐾]𝜕𝑎_𝑖{𝛿_𝑡+▵𝑡}), (5.1) where𝜕{𝐹_𝑡+Δ𝑡}𝜕𝑎_𝑖=𝜕{𝐹_𝑡+Δ𝑡}𝜕𝑎_𝑖+𝜕[𝑀]𝜕𝑎_𝑖(𝑏₀{𝛿_𝑡}+𝑏₂{𝛿̇_𝑡}+𝑏₃{𝛿̈_𝑡})+[𝑀](𝑏₀𝜕{𝛿_𝑡}𝜕𝑎_𝑖+𝑏₂𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖+𝑏₃𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖)+𝜕[𝐶]𝜕𝑎_𝑖(𝑏₁{𝛿_𝑡}+𝑏₄{𝛿̇_𝑡}+𝑏₅{𝛿̈_𝑡})+[𝐶](𝑏₁𝜕{𝛿_𝑡}𝜕𝑎_𝑖+𝑏₄𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖+𝑏₅𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖). (5.2) After 𝜕{𝛿_𝑡}/𝜕𝑎_𝑖=𝑞₀,𝜕{𝛿̇_𝑡}/𝜕𝑎_𝑖=𝑞̇₀, and 𝜕{𝛿̈_𝑡}/𝜕𝑎_𝑖=𝑞̈₀ are given, (5.2) can be calculated.

The partial derivative of (5.1) with respect to 𝑎_𝑗 is given by𝜕²{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗=[𝐾]⁻¹(𝜕²{𝐹_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗−𝜕[𝐾]𝜕𝑎_𝑖𝜕{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑗−𝜕[𝐾]𝜕𝑎_𝑗𝜕{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖−𝜕²[𝐾]𝜕𝑎_𝑖𝜕𝑎_𝑗{𝛿_𝑡+Δ𝑡}), (5.3) where𝜕²{𝐹_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗=𝜕²{𝐹_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝜕²[𝑀]𝜕𝑎_𝑖𝜕𝑎_𝑗(𝑏₀{𝛿_𝑡}+𝑏₂{𝛿̇_𝑡}+𝑏₃{𝛿̈_𝑡+})𝜕[𝑀]𝜕𝑎_𝑖(𝑏₀𝜕{𝛿_𝑡}𝜕𝑎_𝑗+𝑏₂𝜕{𝛿̇_𝑡}𝜕𝑎_𝑗+𝑏₃𝜕{𝛿̈_𝑡}𝜕𝑎_𝑗)+𝜕[𝑀]𝜕𝑎_𝑗(𝑏₀𝜕{𝛿_𝑡}𝜕𝑎_𝑖+𝑏₂𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖+𝑏₃𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖)+[𝑀](𝑏₀𝜕²{𝛿_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₂𝜕²{𝛿̇_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₃𝜕²{𝛿̈_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗)+𝜕²[𝐶]𝜕𝑎_𝑖𝜕𝑎_𝑗(𝑏₁{𝛿_𝑡}+𝑏₄{𝛿̇_𝑡}+𝑏₅{𝛿̈_𝑡+})𝜕[𝐶]𝜕𝑎_𝑖(𝑏₁𝜕{𝛿_𝑡}𝜕𝑎_𝑗+𝑏₄𝜕{𝛿̇_𝑡}𝜕𝑎_𝑗+𝑏₅𝜕{𝛿̈_𝑡}𝜕𝑎_𝑗)+𝜕[𝐶]𝜕𝑎_𝑗(𝑏₁𝜕{𝛿_𝑡}𝜕𝑎_𝑖+𝑏₄𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖+𝑏₅𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖)+[𝐶](𝑏₁𝜕²{𝛿_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₄𝜕²{𝛿̇_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₅𝜕²{𝛿̈_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗). (5.4) After 𝜕{𝛿_𝑡}/𝜕𝑎_𝑗=𝑞₁,𝜕{𝛿̇_𝑡}/𝜕𝑎_𝑗=𝑞̇₁, 𝜕{𝛿̈_𝑡}/𝜕𝑎_𝑗=𝑞̈₁, 𝜕²{𝛿_𝑡}/𝜕𝑎_𝑖𝜕𝑎_𝑗=𝑟₀, and 𝜕²{𝛿̇_𝑡}/𝜕𝑎_𝑖𝜕𝑎_𝑗=𝑟̇₀, 𝜕²{𝛿̈_𝑡}/𝜕𝑎_𝑖𝜕𝑎_𝑗=𝑟̈₀ are given, (5.4) can be calculated.

The displacement is expanded at the mean value point 𝑎=(𝑎₁,𝑎₂,…,𝑎_𝑖,…,𝑎_𝑛1)^𝑇 by means of a Taylor series. the mean of 𝛿_𝑡+▵𝑡 is obtained as𝜇{𝛿_𝑡+Δ𝑡}≈{𝛿_𝑡+Δ𝑡}|_𝑎=𝑎+12^𝑁𝑁_𝑖=1𝑗=1𝜕²{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗|_𝑎=𝑎Cov(𝑎_𝑖,𝑎_𝑗), (5.5) where 𝜇{𝛿_𝑡+Δ𝑡} expresses mean value 𝛿_𝑡+Δ𝑡 and Cov(𝑎_𝑖,𝑎_𝑗) is the covariance between 𝑎_𝑖 and 𝑎_𝑗.

The variance of 𝛿_𝑡+▵𝑡 is given byVar{𝛿_𝑡+Δ𝑡}≈^𝑁𝑁_𝑖=1𝑗=1𝜕{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖|_𝑎=𝑎⋅𝜕{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑗|_𝑎=𝑎⋅Cov(𝑎_𝑖,𝑎_𝑗). (5.6) The velocity vector {𝛿̇_𝑡+▵𝑡} and the acceleration vector {𝛿̈_𝑡+▵𝑡} are given in Appendix B. The partial derivative of 𝛿̈_𝑡+▵𝑡 with respect to 𝑎_𝑖 is given by𝜕{𝛿̈_𝑡+Δ𝑡}𝜕𝑎_𝑖=𝑏₀(𝜕{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖−𝜕{𝛿_𝑡}𝜕𝑎_𝑖)−𝑏₂𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖−𝑏₃𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖. (5.7) The partial derivative of 𝛿̇_𝑡+▵𝑡 with respect to 𝑎_𝑖 is given by𝜕{𝛿̇_𝑡+▵𝑡}𝜕𝑎_𝑖=𝜕{𝛿̇_𝑡}𝜕𝑎_𝑖+𝑏₆𝜕{𝛿̈_𝑡}𝜕𝑎_𝑖+𝑏₇𝜕{𝛿̈_𝑡+▵𝑡}𝜕𝑎_𝑖. (5.8) The partial derivative of (5.7) with respect to 𝑎_𝑗 is given by𝜕²{𝛿̈_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗=𝑏₀(𝜕²{𝛿_𝑡+Δ𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗−𝜕²{𝛿_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗)−𝑏₂𝜕²{𝛿̇_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗−𝑏₃𝜕²{𝛿̈_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗. (5.9) The partial derivative of (5.8) with respect to 𝑎_𝑗 is given by𝜕²{𝛿̇_𝑡+▵𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗=𝜕²{𝛿̇_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₆𝜕²{𝛿̈_𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗+𝑏₇𝜕²{𝛿̈_𝑡+▵𝑡}𝜕𝑎_𝑖𝜕𝑎_𝑗. (5.10) Equations (5.7), (5.8), (5.9) and (5.10), must be calculated for the following iteration.

Then, the mean and variance of displacement are obtained at time 𝑡+𝑖₁Δ𝑡(𝑖₁=2,3,…,𝑛₁) step by step.

The partial derivative of (4.8) with respect to 𝑎_𝑖 is given by𝜕{𝜎}𝜕𝑎_𝑖=𝜕[𝐷]𝜕𝑎_𝑖[𝐵]{𝛿^𝑑_𝑡′}+[𝐷]𝜕[𝐵]𝜕𝑎_𝑖{𝛿^𝑑_𝑡′}+[𝐷][𝐵]𝜕{𝛿^𝑑_𝑡′}𝜕𝑎_𝑖. (5.11) The partial derivative of (5.11) with respect to 𝑎_𝑗 is given by𝜕²{𝜎}𝜕𝑎_𝑖𝜕𝑎_𝑗=𝜕²[𝐷]𝜕𝑎_𝑖𝜕𝑎_𝑗[𝐵]{𝛿^𝑑_𝑡′}+𝜕[𝐷]𝜕𝑎_𝑖𝜕[𝐵]𝜕𝑎_𝑗{𝛿^𝑑_𝑡′}+𝜕[𝐷]𝜕𝑎_𝑖[𝐵]𝜕{𝛿^𝑑_𝑡′}𝜕𝑎_𝑗+𝜕[𝐷]𝜕𝑎_𝑗𝜕[𝐵]𝜕𝑎_𝑖{𝛿^𝑑_𝑡′𝜕}+[𝐷]²[𝐵]𝜕𝑎_𝑖𝜕𝑎_𝑗{𝛿^𝑑_𝑡′}+[𝐷]𝜕[𝐵]𝜕𝑎_𝑖𝜕{𝛿^𝑑_𝑡′}𝜕𝑎_𝑗+𝜕[𝐷]𝜕𝑎_𝑗[𝐵]𝜕{𝛿^𝑑_𝑡′}𝜕𝑎_𝑖+[𝐷]𝜕[𝐵]𝜕𝑎_𝑗𝜕{𝛿^𝑑_𝑡′}𝜕𝑎_𝑖𝜕+[𝐷][𝐵]²{𝛿^𝑑_𝑡′}𝜕𝑎_𝑖𝜕𝑎_𝑗. (5.12) The stress is expanded at mean value point 𝑎=(𝑎₁,𝑎₂,…,𝑎_𝑖,…,𝑎_𝑛1)^𝑇 by means of a Taylor series. By taking the expectation operator for two sides of the above (4.8), the mean of stress is obtained as𝜇{𝜎}≈{𝜎}|_𝑎=𝑎+12^𝑁𝑁_𝑖=1𝑗=1𝜕²{𝜎}𝜕𝑎_𝑖𝜕𝑎_𝑗|_𝑎=𝑎Cov(𝑎_𝑖,𝑎_𝑗), (5.13) where 𝜇{𝜎} expresses the mean of 𝜎 and Cov(𝑎_𝑖,𝑎_𝑗) is the covariance between 𝑎_𝑖 and 𝑎_𝑗.

The variance of 𝜎 is given byVar{𝜎}≈^𝑁𝑁_𝑖=1𝑗=1𝜕{𝜎}𝜕𝑎_𝑖|_𝑎=𝑎⋅𝜕{𝜎}𝜕𝑎_𝑗|_𝑎=𝑎⋅Cov(𝑎_𝑖,𝑎_𝑗). (5.14)

6. Numerical Example

Figure 1 shows a cantilever beam. The length is 1 m, the width is 0.2 m, and the height is 0.05 m. The load subjected to the cantilever beam is 100sin(100t)N. Its material is the concrete. It is divided into 400 rectangle elements that have 505 nodes and 400 midpoints. Young’s modulus is regarded as a stochastic process. For numerical calculation, the means of Young’s modulus at each node and the midpoint within an element are 𝑐₁(1.0+𝜃₁𝑥^′_𝑖/𝐿). Horizontal coordinates of each node and the midpoint within an element are 𝑥^′_𝑖. The covariance of Young’s modulus between any two nodes, between two midpoints and between each node and each midpoint are 𝑐₂(1.0+𝜃₂𝑥_𝑖/𝑙). 𝑐₁,𝑐₂,𝜃₁,𝜃₂,𝑙,𝐿 are constants. The distances between any two nodes, between two midpoints, and between each node and each midpoint are 𝑥_𝑖. Figure 2 shows the mean of vertical displacement at node 505. the DSFEM simulates 100 samples. It is common knowledge that The DSFEM approaches the accurate solution gradually with the increase of the number of simulations. The DSFEM uses the Cholesky decomposition to solve linear equations and provides the reference solution. Figure 3 shows the variance of vertical displacement at node 505. Figure 4 shows the mean of horizontal stress at node 5. Figure 5 shows the variance of horizontal stress at node 5. Table 1 shows results obtained from the PCG, the CG, the TSFEM, the MPCG, the MCG and the MTSFEM compare with those of the DSFEM within six seconds. The PCG adopts the preconditioned Conjugate Gradient method to solve linear equations. The errors of the PCG, the CG and the TSFEM are smaller than those of the MPCG, the MCG, and the MTSFEM. The maximum error is obtained by the MTSFEM.The minimum error is produced by the PCG.Table 2 compares the CPU times of the above-mentioned methods when the cantilever beam has vibrated for six seconds.The PCG requires the least amount of CPU time. The MTSFEM requires the greatest amount of CPU time.

Table 1 

Comparison of error.

Table 2 

Comparison of CPU time.

Figure 1 

A cantilever beam.

Figure 2 

The mean of vertical displacement at node 505.

Figure 3 

The variance of vertical displacement at node 505.

Figure 4 

The mean of horizontal stress at node 5.

Figure 5 

The variance of horizontal stress at node 5.

In order to test accuracy and computational efficiency of the above-mentioned methods, larger covariances of Young’s modulus are selected. Figure 6 shows the mean of vertical displacement at node 505. Figure 7 shows the variance of vertical displacement at node 505. Figure 8 shows the mean of horizontal stress at node 5. Figure 9 shows the variance of horizontal stress at node 5.Table 3 indicates the errors of the above-mentioned methods compare to results of the DSFEM. The results produced by the PCG and the CG are close to those produced by the DSFEM. The TSFEM and the MTSFEM cannot achieve satisfactory results.

Table 3 

Comparison of error for larger covariances.

Figure 6 

The mean of vertical displacement at node 505 for larger covariances.

Figure 7 

The variance of vertical displacement at node 505 for larger covariances.

Figure 8 

The mean of horizontal stress at node 5 for larger covariances.

Figure 9 

The variance of horizontal stress at node 5 for larger covariances.

7. Conclusion

In this paper, improved midpoint method has the advantage of high accuracy. It can be conveniently applied to DSFEM, PSFEM, NSFEM, and CG. The mechanical vibration in a linear system is investigated by using the Taylor expansion. When Young’s modulus is assumed to be a stochastic process, different samples of random variables are simulated. The combination of the CG method and Monte Carlo method makes this an effective method for analyzing a large vibration problem with the characteristics of high accuracy and quick convergence.

Appendices

A. Covariance of Young’s Modulus between Two Elements

The covariance of Young’s modulus between one element containing 𝑚 nodes and another element containing 𝑛 nodes is obtained byCov(𝑎_𝑒,𝑎_𝑓𝑎)=Cov(^𝑒𝑚1+𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎𝑚+1^𝑓𝑛1+𝑎_𝑓𝑛2+⋯+𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙)=1𝑚+1(𝑚+1)(𝑛+1)[Cov(𝑎_𝑒𝑚1,𝑎_𝑓𝑛1+𝑎_𝑓𝑛2+⋯𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙)+Cov(𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎_𝑓𝑛1+𝑎_𝑓𝑛2+⋯𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙=1)](𝑚+1)(𝑛+1)[Cov(𝑎_𝑒𝑚1,𝑎_𝑓𝑛1)+Cov(𝑎_𝑒𝑚1,𝑎_𝑓𝑛2)+⋯+Cov(𝑎_𝑒𝑚1,𝑎_𝑓𝑛𝑛)+Cov(𝑎_𝑒𝑚1,𝑎_𝑓𝑛𝑙+1)](𝑚+1)(𝑛+1)Cov(𝑎_𝑒𝑚2+⋯+𝑎_𝑒𝑚𝑚+𝑎_𝑒𝑚𝑙,𝑎_𝑓𝑛1+𝑎_𝑓𝑛2+⋯+𝑎_𝑓𝑛𝑛+𝑎_𝑓𝑛𝑙)=1((𝑚+1)(𝑛+1)^𝑚_𝑔1^𝑛_=1𝑔3=1Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑓𝑛𝑔3}+1))((𝑚+1)(𝑛+1)^𝑚_𝑔1=1Cov(𝑎_{𝑒𝑚𝑔1},𝑎_𝑓𝑛𝑙+1))((𝑚+1)(𝑛+1)^𝑛_𝑔3=1Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑓𝑛𝑔3}+1))(𝑚+1)(𝑛+1)Cov(𝑎_𝑒𝑚𝑙,𝑎_𝑓𝑛𝑙), (A.1) where Cov(𝑎_{𝑒𝑚𝑔1},𝑎_{𝑓𝑛𝑔3}) = the covariance of Young’s modulus between node 𝑔₁(𝑔₁=1,2,…,𝑚) of element 𝑒 and node 𝑔₃(𝑔₃=1,2,…,𝑛) of element 𝑓, Cov(𝑎_{𝑒𝑚𝑔1},𝑎_𝑓𝑛𝑙) = the covariance of Young’s modulus between node 𝑔₁ of element 𝑒 and the midpoint of element 𝑓, Cov(𝑎_𝑒𝑚𝑙,𝑎_{𝑓𝑛𝑔3}) = the covariance of Young’s modulus between the midpoint of element 𝑒 and node 𝑔₃ of element 𝑓, and Cov(𝑎_𝑒𝑚𝑙,𝑎_𝑓𝑛𝑙) = the covariance of Young’s modulus between the midpoint of element 𝑒 and the midpoint of element 𝑓.

The covariance of Young’s modulus between two elements that each has 𝑛 nodes is given byCov(𝑎_𝑓,𝑎_𝑓′1)=(𝑛+1)²(^𝑛_𝑔3^𝑛_=1𝑔4=1Cov(𝑎_{𝑓𝑛𝑔3},𝑎_{𝑓′𝑛𝑔4}+1))(𝑛+1)²(^𝑛_𝑔3=1Cov(𝑎_{𝑓𝑛𝑔3},𝑎_{𝑓′𝑛𝑙}+1))(𝑛+1)²(^𝑛_𝑔4=1Cov(𝑎_𝑓𝑛𝑙,𝑎_{𝑓′𝑛𝑔4})+𝐶𝑜𝑣(𝑎_𝑓𝑛𝑙,𝑎_{𝑓′𝑛𝑙})), (A.2) where Cov(𝑎_{𝑓𝑛𝑔3},𝑎_{𝑓′𝑛𝑔4}) = the covariance of Young’s modulus between node 𝑔₃(𝑔₃=1,2,…,𝑛) of element 𝑓 and node 𝑔₄(𝑔₄=1,2,…,𝑛) of element 𝑓^′, Cov(𝑎_{𝑓𝑛𝑔3},𝑎_{𝑓′𝑛𝑙}) = the covariance of Young’s modulus between node 𝑔₃ of element 𝑓 and the midpoint of element 𝑓^′, Cov(𝑎_𝑓𝑛𝑙,𝑎_{𝑓′𝑛𝑔4}) = the covariance of Young’s modulus between the midpoint of element 𝑓 and node 𝑔₄ of element 𝑓^′, and Cov(𝑎_𝑓𝑛𝑙,𝑎_{𝑓′𝑛𝑙}) = the covariance of Young’s modulus between the midpoint of element 𝑓 and the midpoint of element 𝑓^′.

B. Newmark Method

For ease of programming, the comprehensive calculation steps of the Newmark method are as follows.

In the initial calculation the matrices [𝐾], [𝑀], and [𝐶] are formed. The initial values {𝛿_𝑡}, {𝛿̇_𝑡}, {𝛿̈_𝑡} are given. After selecting step Δ𝑡 and parameters 𝛾,𝛽, the following relevant parameters are calculated: 𝛾≥0.50,𝛽≥0.25(0.5+𝛾)²,𝑏₀=1𝛽(Δ𝑡)²,𝑏₁=𝛾𝛽Δ𝑡,𝑏₂=1,𝑏𝛽Δ𝑡₃=12𝛽−1,𝑏₄=𝛾𝛽−1,𝑏₅=Δ𝑡2(𝛾𝛽𝑏−2),₆=Δ𝑡(1−𝛾),𝑏₇=𝛾Δ𝑡. (B.1) The stiffness matrix is defined as [𝐾]=[𝐾]+𝑏₀[𝑀]+𝑏₁[𝐶]. (B.2) The stiffness matrix inversion [𝐾]⁻¹ is solved.

Calculation of each step time At time 𝑡+▵𝑡, the load vector is defined as {𝐹_𝑡+Δ𝑡}={𝐹_𝑡+Δ𝑡}+[𝑀](𝑏₀{𝛿_𝑡}+𝑏₂{𝛿̇_𝑡}+𝑏₃{𝛿̈_𝑡})+[𝐶](𝑏₁{𝛿_𝑡}+𝑏₄{𝛿̇_𝑡}+𝑏₅{𝛿̈_𝑡}). (B.3) At time 𝑡+Δ𝑡, the displacement vector is given by {𝛿_𝑡+Δ𝑡}=[𝐾]⁻¹{𝐹_𝑡+Δ𝑡}. (B.4) At time 𝑡+Δ𝑡, the velocity vector and acceleration vector are obtained as {𝛿̈_𝑡+▵𝑡}=𝑏₀({𝛿_𝑡+Δ𝑡}−{𝛿_𝑡})−𝑏₂{𝛿̇_𝑡}−𝑏₃{𝛿̈_𝑡{𝛿̇},_𝑡+▵𝑡}={𝛿̇_𝑡}+𝑏₆{𝛿̈_𝑡}+𝑏₇{𝛿̈_𝑡+▵𝑡}. (B.5) Vectors {𝛿_{𝑡+𝑖1Δ𝑡}}, {𝛿̇_{𝑡+𝑖1▵𝑡}}, and {𝛿̈_{𝑡+𝑖1▵𝑡}} are solved at time 𝑡+𝑖₁Δ𝑡(𝑖₁=2,3,…,𝑛₁) step by step.