# Stochastic Finite Element for Structural Vibration

Authors: Mo Wenhui
Journal: Mathematical Problems in Engineering
 (2010)

---

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


---

## Body

## 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 ofm nodes and n nodes. The Young’s modulus of node within element e of m nodes is expressed by aem1,aem2,…,aemm. The Young’s modulus of midpoint within element e of m nodes is expressed by aeml. The Young’s modulus of node within element f of n nodes is expressed by afn1,afn2,…,afnn. The Young’s modulus of midpoint within element f of n nodes is expressed by afnl.The Young’s modulus within elemente is defined as(2.1)ae=aem1+aem2+⋯+aemm+aemlm+1.Its mean is(2.2)μe=μ(aem1+aem2+⋯+aemm+aemlm+1)=μaem1+μaem2+⋯+μaemm+μaemlm+1,
where the means of Young’s modulus at the first node, the second node,…, the mth node of element e are expressed by μaem1,μaem2,…,μaemm.The mean of Young’s modulus at the midpoint of element e is expressed by μaeml.The Young’s modulus within elementf is defined as(2.3)af=afn1+afn2+⋯+afnn+afnln+1.Its mean is(2.4)μf=μ(afn1+afn2+⋯+afnn+afnln+1)=μafn1+μafn2+⋯+μafnn+μafnln+1,
where the means of Young’s modulus at the first node, the second node,…,the nth node of element f are expressed by μafn1,μafn2,…,μafnn. The mean of Young’s modulus at the midpoint of element f is expressed by μafnl.The covariance of Young’s modulus between two elements that each hasm nodes is obtained by(2.5)Cov(ae,ae′)=Cov(aem1+aem2+⋯+aemm+aemlm+1,ae′m1+ae′m2+⋯+ae′mm+ae′mlm+1)=1(m+1)2[Cov(aem1,ae'm1+ae'm2+⋯+ae'mm+ae'ml)+Cov(aem2+⋯aemm+aeml,ae′m1+ae′m2+⋯+ae′mm+ae′ml)]=1(m+1)2[Cov(aem1,ae′m1)+Cov(aem1,ae′m2)+⋯+Cov(aem1,ae′mm)+Cov(aem1,ae′ml)]+1(m+1)2Cov(aem2+⋯+aemm+aeml,ae′m1+ae′m2+⋯+ae′mm+ae′ml)=1(m+1)2(∑g1=1m∑g2=1mCov(aemg1,ae′mg2))+1(m+1)2(∑g1=1mCov(aemg1,ae′ml))+1(m+1)2(∑g2=1mCov(aeml,ae′mg2)+Cov(aeml,ae′ml)),
where Cov(aemg1,ae'mg2) = the covariance of Young’s modulus between node g1(g1=1,2,…,m) of element e and node g2(g2=1,2,…,m) of element e′, Cov(aemg1,ae'ml) = the covariance of Young’s modulus between node g1 of element e and the midpoint of element e′, Cov(aeml,ae′mg2)  = the covariance of Young’s modulus between the midpoint of element e and node g2 of element e′, and Cov(aeml,ae′ml) = the covariance of Young’s modulus between the midpoint of element e and the midpoint of element e′.The covariance of Young’s modulus between two elements is given in AppendixA.Using covariance matrix, the correlation of Young’s modulus between any two elements is given by(2.6)Caa=(Cov(a1,a1)Cov(a1,a2)⋯Cov(a1,aN)Cov(a2,a1)Cov(a2,a2)⋯Cov(a2,aN)⋮⋮⋮Cov(aN,a1)Cov(aN,a2)⋯Cov(aN,aN)).A Gaussian vectora⇀=[a1,a2,…,aN]T is generated(2.7)a⇀=LZ.Z=[Z1,Z2,…,ZN]T consists of N Gaussian random variables with mean zero and unit standard deviation. The Cholesky matrix L can be obtained through a decomposition of the covariance matrix; therefore,(2.8)μ[ZZT]=I,LLT=Caa.I is the identity matrix. The generation of vector a⇀ must satisfy the covariance matrix(2.9)μ[a⇀a⇀T]=μ[LZ(LZ)T]=Lμ[ZZT]LT=Caa.Once the decomposition has been completed, different samples of vectora⇀ 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)[M]{δ̈}+[C]{δ̇}+[K]{δ}={F},
where {δ̈},{δ̇},{δ} are the acceleration, velocity, and displacement vectors. [M],[K], and [C] 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){δt+Δt}=[K̃]-1{F̃t+Δt},
where {δt+Δt}，[K̃] and {F̃t+Δt} indicate the displacement vector, stiffness matrix and load vector at time t+Δt. [K̃] and {F̃t+Δt} are given in Appendix B.


## 4. Dynamic Analysis of Structure Based on CG

Equation (3.2) can be rewritten as(4.1)[K̃]{δt+∆t}={F̃t+∆t}.
Using (2.7), (2.8), and (2.9), N1 samples of vector a⇀ are produced. N1 matrices [K̃] and N1 (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(4.2){δt+Δt}(0)=(δ(0)(t+Δt)1,δ(0)(t+Δt)2,…,δ(0)(t+Δt)N1)T;calculate the first residual vector(4.3)r(0)={F̃t+Δt}-[K̃]{δt+Δt}(0)
and vector
(4.4)p(0)=[K̃]Tr(0)
where [K̃]Tis the transposed matrix;forĩ=0,1,2,…,n2-1, iterate step by step as follows:
(4.5)αĩ=([K̃]p(ĩ),r(ĩ))([K̃]p(ĩ),[K̃]p(ĩ))=(p(ĩ),[K̃]Tr(ĩ))([K̃]p(ĩ),[K̃]p(ĩ))=([K̃]Tr(ĩ),[K̃]Tr(ĩ))([K̃]p(ĩ),[K̃]p(ĩ)),{δt+∆t}(ĩ+1)={δt+∆t}(ĩ)+αĩp(ĩ),r(ĩ+1)=r(ĩ)-αĩ[K̃]p(ĩ),βĩ+1=([K̃]Tr(ĩ+1),[K̃]Tr(ĩ+1))([K̃]Tr(ĩ),[K̃]Tr(ĩ)),p(ĩ+1)=[K̃]Tr(ĩ+1)+βĩ+1p(ĩ).This process stops only ifrn2 is small enough.Vectors{δt+Δt}1,{δt+Δt}2,…,{δt+Δt}N1are solutions of N1 (4.1).The mean of{δt+Δt} is given by(4.6)μ{δt+Δt}={δt+Δt}1+{δt+Δt}2+⋯+{δt+Δt}N1N1.
The variance of {δt+Δt} is given by(4.7)Var{δt+Δt}=1N1-1∑i=1N1({δt+Δt}i-μ{δt+Δt})2.
Similarly, the mean and variance of vector {δt+i1Δt} can be solved at time t+i1Δt step by step.At timet′=t+i2Δt(i2=1,2,…,n1), the stress for element d is given by(4.8){σ}=[D][B]{δt′d},
where [D] = the material response matrix of element d, [B] = the gradient matrix of element d and {δt'd} = the element d nodal displacement vector at time t′.SubstitutingN1 samples of vector a⇀ into (4.8), vectors {σ}1,{σ}2,…,{σ}N1can be obtained.The mean of{σ} is given by(4.9)μ{σ}={σ}1+{σ}2+⋯+{σ}N1N1.
The variance of {σ} is given by(4.10)Var{σ}=1N1-1∑i=1N1({σ}i-μ{σ})2.
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 ofN random variables a1,a2,…,aN. The partial derivative of (4.1) with respect to ai is given by(5.1)∂{δt+∆t}∂ai=[K̃]-1(∂{F̃t+∆t}∂ai-∂[K̃]∂ai{δt+∆t}),
where(5.2)∂{F̃t+Δt}∂ai=∂{Ft+Δt}∂ai+∂[M]∂ai(b0{δt}+b2{δ̇t}+b3{δ̈t})+[M](b0∂{δt}∂ai+b2∂{δ̇t}∂ai+b3∂{δ̈t}∂ai)+∂[C]∂ai(b1{δt}+b4{δ̇t}+b5{δ̈t})+[C](b1∂{δt}∂ai+b4∂{δ̇t}∂ai+b5∂{δ̈t}∂ai).
After ∂{δt}/∂ai=q0,∂{δ̇t}/∂ai=q̇0, and ∂{δ̈t}/∂ai=q̈0 are given, (5.2) can be calculated.The partial derivative of (5.1) with respect to aj is given by(5.3)∂2{δt+Δt}∂ai∂aj=[K̃]-1(∂2{F̃t+Δt}∂ai∂aj-∂[K̃]∂ai∂{δt+Δt}∂aj-∂[K̃]∂aj∂{δt+Δt}∂ai-∂2[K̃]∂ai∂aj{δt+Δt}),
where(5.4)∂2{F̃t+Δt}∂ai∂aj=∂2{Ft+Δt}∂ai∂aj+∂2[M]∂ai∂aj(b0{δt}+b2{δ̇t}+b3{δ̈t})+∂[M]∂ai(b0∂{δt}∂aj+b2∂{δ̇t}∂aj+b3∂{δ̈t}∂aj)+∂[M]∂aj(b0∂{δt}∂ai+b2∂{δ̇t}∂ai+b3∂{δ̈t}∂ai)+[M](b0∂2{δt}∂ai∂aj+b2∂2{δ̇t}∂ai∂aj+b3∂2{δ̈t}∂ai∂aj)+∂2[C]∂ai∂aj(b1{δt}+b4{δ̇t}+b5{δ̈t})+∂[C]∂ai(b1∂{δt}∂aj+b4∂{δ̇t}∂aj+b5∂{δ̈t}∂aj)+∂[C]∂aj(b1∂{δt}∂ai+b4∂{δ̇t}∂ai+b5∂{δ̈t}∂ai)+[C](b1∂2{δt}∂ai∂aj+b4∂2{δ̇t}∂ai∂aj+b5∂2{δ̈t}∂ai∂aj).
After ∂{δt}/∂aj=q1,∂{δ̇t}/∂aj=q̇1, ∂{δ̈t}/∂aj=q̈1, ∂2{δt}/∂ai∂aj=r0, and ∂2{δ̇t}/∂ai∂aj=ṙ0, ∂2{δ̈t}/∂ai∂aj=r̈0 are given, (5.4) can be calculated.The displacement is expanded at the mean value pointa̅=(a̅1,a̅2,…,a̅i,…,a̅n1)T by means of a Taylor series. the mean of δt+∆t is obtained as(5.5)μ{δt+Δt}≈{δt+Δt}|a=a̅+12∑i=1N∑j=1N∂2{δt+Δt}∂ai∂aj|a=a̅Cov(ai,aj),
where μ{δt+Δt} expresses mean value δt+Δt and Cov(ai,aj) is the covariance between ai and aj.The variance ofδt+∆t is given by(5.6)Var{δt+Δt}≈∑i=1N∑j=1N∂{δt+Δt}∂ai|a=a̅·∂{δt+Δt}∂aj|a=a̅·Cov(ai,aj).
The velocity vector {δ̇t+∆t} and the acceleration vector {δ̈t+∆t} are given in Appendix B. The partial derivative of δ̈t+∆t with respect to ai is given by(5.7)∂{δ̈t+Δt}∂ai=b0(∂{δt+Δt}∂ai-∂{δt}∂ai)-b2∂{δ̇t}∂ai-b3∂{δ̈t}∂ai.
The partial derivative of δ̇t+∆t with respect to ai is given by(5.8)∂{δ̇t+∆t}∂ai=∂{δ̇t}∂ai+b6∂{δ̈t}∂ai+b7∂{δ̈t+∆t}∂ai.
The partial derivative of (5.7) with respect to aj is given by(5.9)∂2{δ̈t+Δt}∂ai∂aj=b0(∂2{δt+Δt}∂ai∂aj-∂2{δt}∂ai∂aj)-b2∂2{δ̇t}∂ai∂aj-b3∂2{δ̈t}∂ai∂aj.
The partial derivative of (5.8) with respect to aj is given by(5.10)∂2{δ̇t+∆t}∂ai∂aj=∂2{δ̇t}∂ai∂aj+b6∂2{δ̈t}∂ai∂aj+b7∂2{δ̈t+∆t}∂ai∂aj.
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 timet+i1Δt(i1=2,3,…,n1) step by step.The partial derivative of (4.8) with respect to ai is given by(5.11)∂{σ}∂ai=∂[D]∂ai[B]{δt′d}+[D]∂[B]∂ai{δt′d}+[D][B]∂{δt′d}∂ai.
The partial derivative of (5.11) with respect to aj is given by(5.12)∂2{σ}∂ai∂aj=∂2[D]∂ai∂aj[B]{δt′d}+∂[D]∂ai∂[B]∂aj{δt′d}+∂[D]∂ai[B]∂{δt′d}∂aj+∂[D]∂aj∂[B]∂ai{δt′d}+[D]∂2[B]∂ai∂aj{δt′d}+[D]∂[B]∂ai∂{δt′d}∂aj+∂[D]∂aj[B]∂{δt′d}∂ai+[D]∂[B]∂aj∂{δt′d}∂ai+[D][B]∂2{δt′d}∂ai∂aj.
The stress is expanded at mean value point a̅=(a̅1,a̅2,…,a̅i,…,a̅n1)T 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(5.13)μ{σ}≈{σ}|a=a̅+12∑i=1N∑j=1N∂2{σ}∂ai∂aj|a=a̅Cov(ai,aj),
where μ{σ} expresses the mean of σ and Cov(ai,aj) is the covariance between ai and aj.The variance ofσ is given by(5.14)Var{σ}≈∑i=1N∑j=1N∂{σ}∂ai|a=a̅·∂{σ}∂aj|a=a̅·Cov(ai,aj).


## 6. Numerical Example

Figure1 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 c1(1.0+θ1xi′/L). Horizontal coordinates of each node and the midpoint within an element are xi′. The covariance of Young’s modulus between any two nodes, between two midpoints and between each node and each midpoint are c2(1.0+θ2xi/l). c1,c2,θ1,θ2,l,L are constants. The distances between any two nodes, between two midpoints, and between each node and each midpoint are xi. 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.
The mean of vertical displacement at node 505The variance of vertical displacement at node 505The mean of horizontal stress at node 5The variance of horizontal stress at node 5PCG2.95%3.42%4.31%5.07%CG3.11%3.84%4.57%5.24%TSFEM5.12%6.01%6.28%7.17%MPCG7.14%7.02%8.67%9.94%MCG7.23%7.63%8.82%10.17%MTSFEM12.76%13.27%14.47%16.53%Table 2
Comparison of CPU time.
DSFEMPCGCGTSFEMMPCGMCGMTSFEMCPU time4 h 32 m 47 s1 h 15 m 43 s1 h 47 m 36 s10 h 45 m 34 s1 h 7 m 23 s1 h 36 m 31 s4 h 12 m 17 sh: hour; m: minute; s: second.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. Figure6 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.
The mean of vertical displacement at node 505The variance of vertical displacement at node 505The mean of horizontal stress at node 5The variance of horizontal stress at node 5PCG9.27%12.15%11.43%14.16%CG9.53%12.37%11.67%14.41%TSFEM43.19%64.28%56.42%76.39%MPCG31.42%45.13%47.16%63.17%MCG31.75%45.41%47.43%63.42%MTSFEM94.27%125.74%125.41%162.43%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.



---

*Source: 584863-2010-08-05.xml*