Structural stochastic responses determination via a sample-based stochastic finite element method
SStructural stochastic responses determination via asample-based stochastic finite element method
Zhibao Zheng a,b, ∗ a School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China b Key Lab of Structures Dynamic Behavior and Control, Harbin Institute of Technology,Ministry of Education, Harbin 150090, China
Abstract
This paper presents a new stochastic finite element method for computingstructural stochastic responses. The method provides a new expansion ofstochastic response and decouples the stochastic response into a combina-tion of a series of deterministic responses with random variable coefficients.A dedicated iterative algorithm is proposed to determine the deterministicresponses and corresponding random variable coefficients one by one. Thealgorithm computes the deterministic responses and corresponding randomvariable coefficients in their individual space and is insensitive to stochasticdimensions, thus it can be applied to high dimensional stochastic problemsreadily without extra difficulties. More importantly, the deterministic re-sponses can be computed efficiently by use of existing Finite Element Method(FEM) solvers, thus the proposed method can be easy to embed into exist-ing FEM structural analysis softwares. Three practical examples, includ-ing low-dimensional and high-dimensional stochastic problems, are given to ∗ Corresponding author
Email address: [email protected] (Zhibao Zheng)
Preprint submitted to Elsevier February 11, 2021 a r X i v : . [ m a t h . NA ] F e b emonstrate the accuracy and effectiveness of the proposed method. Keywords:
Stochastic finite element method; High dimensions; Largescale; Stochastic responses;
1. Introduction
Ordinary or partial differential equations (PDEs) are powerful tools todescribe many real life engineering and scientific processes. A wide body ofnumerical methods based on finite differences, finite elements, and boundaryelements are available to approximately solve the governing equations for theresponse quantities of interest. In particular, Finite element method (FEM)has become state-of-the-art since it offers a simple way to solve very highresolution models in various computational physics problems, ranging fromstructural mechanics, thermodynamics to nano-bio mechanics [1]. Neverthe-less, the FEM is deterministic in nature and is therefore limited to describethe general characteristics of a real life system. The considerable influence ofinherent uncertainties on system behavior has led the scientific communityto recognize the importance of a stochastic approach to realistic engineeringsystems [2]. More than ever, the goal becomes to represent and propagateuncertainties from the available data to the desired results through PDEswithin the framework of stochastic equations [3, 4].The modelling of uncertainties consists in defining a suitable probabilityspace (Θ , Σ , P ), where Θ denotes the space of elementary events, Σ is a σ -algebra defined on Θ and P is a probability measure. In this paper, weconsider the structural stochastic response u ( θ ) of the problem is a stochasticfunction, with value in a certain function space, which has to verify almost2urely the stochastic partial differential equations (SPDEs) discretized bystochastic finite element equations [5] as K ( θ ) u ( θ ) = F ( θ ) (1)where K ( θ ) is an operator representing properties of the physical modelunder investigation, which can be considered the stochastic stiffness matrixand F ( θ ) is a stochastic load vector. Randomness on the model can beformalized as a dependency of the operator and loads on the elementaryevent θ ∈ Θ and it’s a great challenge to solve Equation (1) in the highdimensional stochastic space Θ.As an extension of deterministic FEM, stochastic finite element method(SFEM) [6, 7] has become a common tool for the solution of Equation (1).Given the representation of uncertain system parameters and environmentalsource in terms of random fields, it becomes possible to integrate discretiza-tion methods for the response and random fields to arrive at a system ofrandom algebraic equations. Two prominent variants of the SFEM are thenon-intrusive methods and the Galerkin-type methods. Although variousnon-intrustive methods, e.g., Monte Carlo simulation [8], or regression andprojection methods [9], can be readily applied to compute the response statis-tics to an arbitrary degree of accuracy, this is the method of last resort sincethe attendant computational cost can be prohibitive for real life problems.The Galerkin-type spectral methods [5, 6, 10], which are developed forlinear SFEM, provide an explicit functional relationship between the ran-dom input and output, hence allow easy evaluation of the statistics of thestochastic system response. These methods transform Equation (1) aris-ing from spatial discretization of SPDEs into a deterministic finite element3quation by stochastic Galerkin projection, but the size of the determinis-tic finite element equation is significantly higher than that of the originalSPDEs. Although several iterative solvers have been developed to decreasethe substantial computational requisite [11, 12], the difficulty to build efficientpreconditioners and memory requirements still limit their use to small-scaleand low-dimensional stochastic problems. Also, the Curse of Dimensionalityin stochastic spaces makes these methods more inefficient. For this line ofapproach to be successful in practice, it is crucial to have general-purposeand highly efficient numerical schemes for the solution of stochastic finiteelement equation (1).In this article, we develop a highly efficient numerical method for the ex-plicit and high precision solution of Equation (1) with application to struc-tural responses that involve uncertainties. An universal construct of solutionto stochastic finite element (SFE) equations is firstly developed, which isindependent on the types of SFE problems. Based on the construct of thissolution, we further develop a numerical algorithm for solving SFE equa-tions. The representations of the stochastic solutions are applicable forhigh-dimensional stochastic problems, and more importantly, the stochas-tic analysis and deterministic analysis in the solution procedure can thus beimplemented in individual space. In this way, the proposed algorithm forthe solution of SFE equation integrate the advantages of the non-instrusivemethods and the Galerkin-type methods simultaneously, and thus have greatpotential for uncertainty quantifications in structural analysis.The paper is organized as follows: Section 2 briefly introduces the seriesexpansion methods of random fields simulation and the derivation of stochas-4ic finite element equations. A new method outlines for solving stochasticfinite element equations is described in Section 3. Following this, the al-gorithm implementation of the proposed method is elaborated in Section 4.Three practical problems are used to demonstrate the accuracy and effective-ness of the proposed method in Section 5. Some conclusions and prospectsare discussed in Sections 6.
2. Stochastic finite element method
In the framework of SFEM for structural analysis, uncertain physical pa-rameters usually consists of the Young modulus, Poisson’s ratio, yield stress,cross section geometry of physical systems, earthquake loading, wind loads,etc. In most cases, due to the lack of relevant experimental data, assump-tions are made regarding probabilistic characteristics of random fields, sunchas Gaussian or non-Gaussian, stationary or non-stationary, etc [13, 14]. Thefirst step in applying the FEM to problems involving one or more of therandom parameters is to model random fields based on the assumptions ofprobabilistic characteristics, thus random fields discretization is a key stepin the numerical solutions of stochastic finite element equations. In order toderive stochastic finite element equations effectively, explicit expressions ofrandom fields are also crucial. In general, we represent random fields by anenumerable set of random variables, and the series expansion of a second-order random field ω ( x, θ ), which is indexed on a bounded domain D , can beexpressed as ω ( x, θ ) = M (cid:88) i =0 ξ i ( θ ) ω i ( x ) (2)5here { ξ i ( θ ) } Mi =0 and { ω i ( x ) } Mi =0 are random variables and deterministic func-tions, respectively, and M is the number of retained items. Equation (2) canbe obtained by some methods for discretization of random fields. Variousdiscretization techniques are available in the literature for approximatingrandom fields including shape function methods, optimal linear estimation,weighted integral methods, orthogonal series expansion [15, 16].As a special case of the orthogonal series expansion, Karhunen-Lo´eveexpansion is the most commonly used method in SFEM, and it has a formas ω ( x, θ ) = ω ( x ) + M (cid:88) i =1 ξ i ( θ ) (cid:112) λ i ω i ( x ) (3)where ω ( x ) is the mean function of the random field ω ( x, θ ), { λ i } Mi =1 and { ω i ( x ) } Mi =1 are eigenvalues and eigenfunctions of the covariance function C ωω ( x , x ) of the random field ω ( x, θ ), and they are solutions of the ho-mogenous Fredholm integral equation of the second kind [6, 17], (cid:90) D C ωω ( x , x ) ω i ( x ) dx = λ i ω i ( x ) (4)Due to the symmetry and the positive definiteness of covariance kernel C ωω ( x , x ),the eigenfunctions { ω i ( x ) } Mi =1 form a complete orthogonal set satisfying theequation (cid:90) D ω i ( x ) ω j ( x ) dx = δ ij (5)where δ ij is the Kronecker delta function. An explicit expression for therandom variables { ξ i ( θ ) } Mi =1 in Equation (2) can be obtained by ξ i ( θ ) = 1 √ λ i (cid:90) D [ ω ( x, θ ) − ω ( x )] ω i ( x ) dx (6)6hich is a set of uncorrelated standardized random variables and satisfy E { ξ i ( θ ) } = 0 , E { ξ i ( θ ) ξ j ( θ ) } = δ ij (7)where E {·} is the expectation operator.The Karhunen-Lo`eve expansion (3) offers a unified and powerful tool [17,18] for representing stationary and nonstationary, Gaussian and non-Gaussianrandom fields with explicitly known covariance functions. Karhunen-Lo`eveexpansion is optimal among series expansion methods in the global meansquare error with respect to the number of random variables in the repre-sentation, which means that only a few terms M are required in order tocapture most of randomness, thus it has received much attentions in manydisciplines. In stochastic finite element analysis, it has been widely used todiscretize the random fields representing the randomness of structures and ex-citations. It is worth mentioning that the implementation of Karhunen-Lo`eveexpansion requires solutions of the integral equation (4) with the covariancefunction as the integral kernel. Although only a limited number of analyticaleigen-solutions are available [6], the solution of the integral equation can benumerically approximated for random fields with arbitrary covariance func-tions. For random fields that are defined on two- and three-dimensional do-mains, the finite element method becomes the only available method for thediscretization of the multi-dimensional integral eigenvalue problems [19, 20].In this paper, the generation of the finite element mesh for random fields issame to that for responses. 7 .2. Stochastic finite element equations We simply recall the deterministic finite element method of the relevantformulation before dealing with stochastic problems. The deterministic finiteelement method in linear elasticity defined on Ω eventually derive a N × N linear system Ku = F (8)where N is the number of degrees of freedom, K , u , F are global stiffnessmatrix, displacement vector and load vector, respectively. By assembling theelement stiffness matrices k e , the global stiffness matrix K can be obtainedas k e = (cid:90) Ω e B T DBd Ω e (9)where B and D stand for the strain matrix and the elasticity matrix, respec-tively.We suppose that the material Young’s modulus is a random field [15] andcan be written as the form in Equation (3), D ( x, θ ) = D (cid:34) ω ( x ) + M (cid:88) i =1 ξ i ( θ ) (cid:112) λ i ω i ( x ) (cid:35) (10)where D is a constant matrix, and random variables { ξ i ( θ ) } Mi =1 construct a M -dimensional stochastic space. By substituting Equation (10) into Equa-tion (9), the element stiffness matrix thus becomes as, k e ( θ ) = k e + M (cid:88) i =1 ξ i ( θ ) k ei (11)where k e is the mean element stiffness matrix given by k e = (cid:90) Ω e ω ( x ) B T D Bd Ω e (12)8nd k ei are deterministic matrices given by k ei = (cid:90) Ω e (cid:112) λ i ω i ( x ) B T D Bd Ω e (13)The stochastic global stiffness matrix K ( θ ) in the stochastic finite elementequation (1) is obtained by assembling the stochastic element stiffness ma-trices k e ( θ ), K ( θ ) = M (cid:88) i =0 ξ i ( θ ) K i (14)where ξ ( θ ) ≡ K i are obtained by assembling elementmatrices k ei in the way similar to the deterministic case. In a similar way, wecan get the stochastic global load vector as F ( θ ) = Q (cid:88) l =0 η l ( θ ) F l (15)After assembling the stochastic global stiffness matrix K ( θ ) and the stochas-tic global load vector F ( θ ), the stochastic finite element equation (1) can berewritten as (cid:32) M (cid:88) i =0 ξ i ( θ ) K i (cid:33) u ( θ ) = Q (cid:88) l =0 η l ( θ ) F l (16)The high precision solution of Equation (16) is one of the most impor-tant problems of the stochastic finite element method. Spectral stochasticfinite element method (SSFEM) [21, 6] is a popular method in the past fewdecades. SSFEM represents the stochastic response u ( θ ) through polyno-mial chaos expansion (PCE) and transform Equation (16) into a determin-istic finite element equation by stochastic Galerkin projection. The size ofthe deterministic finite element equation depends directly on the number ofterms retained in the PCE and the number of degrees of freedom N , and the9omputational cost for the solution of this system is much larger than thatof the original problem. Although several improved methods [11, 22] havebeen developed to decrease computational costs, the Curse of Dimensionalitystill limit SSFEM to low-dimensional stochastic problems, thus it is crucialto develop a new method for the solution of Equation (16).
3. A new method for solving stochastic finite element equations
In order to avoid the difficulties of SSFEM, in this section, we propose anew method for solving the sochastic finite element equation (16) defined inlow- and high-dimensional stochastic spaces. A natural idea is to representthe stochastic solution u ( θ ) of Equation (16) by use of random field expan-sions, however common methods are inactive since we almost know nothingabout u ( θ ) except the governing equation (16). Inspired by Karhunen-Lo`eveexpansion (3) and the general spectral decomposition [23], we construct the u ( θ ) as u ( θ ) = ∞ (cid:88) i =1 λ i ( θ ) d i (17)where { λ i ( θ ) } ∞ i =1 are random variables and { d i } ∞ i =1 are deterministic dis-cretized basis vectors. Similar to the orthogonal conditions Equation (5)and (7) of Karhunen-Lo`eve expansion, the following bi-orthogonal conditionis introduced d Ti d j = δ ij , E { λ i ( θ ) λ j ( θ ) } = κ i δ ij (18)where E {·} is the expectation operator and κ i = E { λ i ( θ ) } .It is shown in expansion (17) that the solution space of u ( θ ) is decoupledinto a stochastic space and a deterministic space and it allows to compute10 λ i ( θ ) } ∞ i =1 in the stochastic space and { d i } ∞ i =1 in the deterministic space, re-spectively. In this way, the difficulties in expanding the unknown solutionrandom field of Equation (1) can be overcome. One only requires to seeka set of deterministic orthogonal vectors { d i } ∞ i =1 and the corresponding un-correlated random variables { λ i ( θ ) } ∞ i =1 such that the expanded solution inEquation (17) satisfies the Equation (1). In practical, we truncate Equation(17) at the k -th term as, u k ( θ ) = k (cid:88) j =1 λ j ( θ ) d j (19)As mentioned above, neither { d i } ki =1 nor { λ i ( θ ) } ki =1 is known a priori, a nat-ural choice is to successively determine these unknown couples { λ i ( θ ) , d i } one after another via iterative methods. In order to compute the couple( λ k ( θ ) , d k ), we suppose that the approximate solution u k − ( θ ) has been ob-tained, then substituting Equation (19) into Equation (1) yields, K ( θ ) (cid:34) k − (cid:88) j =1 λ j ( θ ) d j + λ k ( θ ) d k (cid:35) = F ( θ ) (20)If random variable λ k ( θ ) has been determined (or given an initial value),the d k can be determined using stochastic Galerkin method and a dedicatediteration [23], this corresponds E (cid:40) λ k ( θ ) K ( θ ) (cid:34) k − (cid:88) j =1 λ j ( θ ) d j + λ k ( θ ) d k (cid:35)(cid:41) = E { λ k ( θ ) F ( θ ) } (21)Considering Equation (14) and (15), the Equation (21) about d k can besimplified as, (cid:32) M (cid:88) i =0 c ikk K i (cid:33) d k = Q (cid:88) l =0 b kl F l − M (cid:88) i =0 k − (cid:88) j =1 c ijk K i d j (22)11here c ijk = E { ξ i ( θ ) λ j ( θ ) λ k ( θ ) } , b kl = E { η l ( θ ) λ k ( θ ) } (23)Once d k has been determined in Equation (22), the random variable λ k ( θ )can be subsequently updated via the similar procedure. This requires tomultiply d k on both sides of Equation (20) to yield d Tk K ( θ ) (cid:34) k − (cid:88) j =1 λ j ( θ ) d j + λ k ( θ ) d k (cid:35) = d Tk F ( θ ) (24)Considering Equation (14) and (15), the Equation (24) about λ k ( θ ) can besimplified as, (cid:32) M (cid:88) i =0 g ikk ξ i ( θ ) (cid:33) λ k ( θ ) = Q (cid:88) l =0 h kl η l ( θ ) − M (cid:88) i =0 k − (cid:88) j =1 g ijk ξ i ( θ ) λ j ( θ ) (25)where g ijk = d Tk K i d j , h kl = d Tk F l (26)The classical SSFEM is to represent the stochastic solution of nodes { u i ( θ ) } Ni =1 in terms of a set of polynomial chaos and transforms the originalstochastic finite element equation into a deterministic finite element equa-tion with size N × ( M + p )! M ! p ! , where ( · )! represents the factorial operator, N , M and p are the number of system degrees of freedom, the number of randomvariables and the order of polynomial chaos expansion, respectively. The sizeof the deterministic finite element equation is significantly higher than thatof the original stochastic finite element equation. For instance, the size is1 × when N = 1000, M = 10 and p = 4, which leads to the Curse of Di-mensionality, and is prohibitive for problems with high stochastic dimensionsand large scales. 12he method in this paper decouples the original stochastic finite elementequation into a deterministic finite element equation (22) with size N andone-dimensional stochastic algebraic equation (25). The iteration process ofEquation (21) and (24) was used in the paper [23] to solve linear stochas-tic partial differential equations and subsequently applied to time-dependentand nonlinear problems [24, 25]. The key of this method is to transform theoriginal SPDE to a deterministic PDE and a stochastic algebraic equationlike Equation (25). The method for solving Equation (25)-like is to repre-sent the random variable λ k ( θ ) in terms of a set of polynomial chaos andtransforms Equation (25)-like into a deterministic equation with size ( M + p )! M ! p ! ,which greatly alleviates the Curse of Dimensionality, but is still prohibitivefor problems with high stochastic dimensions. In order to avoid this difficulty,we develop a simulation method to determine λ k ( θ ). For each realization of { θ ( r ) } Rr =1 , the λ k ( θ ( r ) ) can be obtained by solving (25) as, λ k (cid:0) θ ( r ) (cid:1) = Q (cid:80) l =0 h kl η l (cid:0) θ ( r ) (cid:1) − M (cid:80) i =0 k − (cid:80) j =1 g ijk ξ i (cid:0) θ ( r ) (cid:1) λ j (cid:0) θ ( r ) (cid:1) M (cid:80) i =0 g ikk ξ i ( θ ( r ) ) (27)It is important to note that Equation (27) has become a one-dimensionallinear algebraic equation about λ k ( θ ( r ) ). Compared to classic methods, wedo not need to choose the type and order of polynomial chaos. The totalcomputational cost for determining { λ k ( θ ( r ) ) } Rr =1 is very low even for highstochastic dimensions, which hopefully avoid the Curse of Dimensionality.Then statistical methods are readily introduced to obtain λ k ( θ ) from samples { λ k ( θ ( r ) ) } Rr =1 . Hence, this method will be particularly appropriate for a wideclass of high-dimensional stochastic problems in practice.13 . Algorithm implementationAlgorithm 1 Algorithm for solving linear stochastic finite element equations while ε global > ε do initialize λ (0) k ( θ ) while ε local > ε do compute d ( j ) k by solving Equation (22) orthogonalization d ( j ) k ⊥ d i , i = 1 , · · · , k − d ( j ) k = d ( j ) k (cid:13)(cid:13)(cid:13) d ( j ) k (cid:13)(cid:13)(cid:13) compute λ ( j ) k ( θ ) by Equation (27) orthogonalization λ ( j ) k ( θ ) ⊥ λ i ( θ ) , i = 1 , · · · , k − compute local error ε local , j = j + 1 update u ( θ ) as u k ( θ ) = k − (cid:80) i =1 λ i ( θ ) d i + λ k ( θ ) d k compute global error ε global , k = k + 1The above procedure for solving the stochastic finite element equation (1)is summarized in Algorithm 1, which consists of a outer loop procedure anda inner loop procedure. The inner loop, which is from step 3 to 8, is used todetermine the couple of ( λ k ( θ ) , d k ). With an initial random variable λ (0) k ( θ )given in step 2, d ( j ) k can be determined in step 4 and 5, where superscript j represents the j -th round of iteration. With the obtained d ( j ) k , the randomvariable λ ( j ) k ( θ ) is then updated in step 6 and 7. While the outer loop, whichis from step 1 to 10, corresponds to recursively building the set of couplesand then generates a set of couples such that the approximate solution instep 9 satisfies Equation (1). 14ote that both d ( j ) k and λ ( j ) k ( θ ) require orthogonalizations such that thebi-orthogonal condition in Equation (18) holds along the whole process, herewe use the Gram-Schmidt Orthogonalization method in step 5 and 7. It iswritten as, d ( j ) k = d ( j ) k − k − (cid:80) i =1 (cid:16) d ( j ) Tk d i (cid:17) d i λ ( j ) k ( θ ) = λ ( j ) k ( θ ) − k − (cid:80) i =1 E (cid:110) λ ( j ) k ( θ ) λ i ( θ ) (cid:111) E { λ i ( θ ) } λ i ( θ ) , k ≥ u ( θ ). The truncation criterion in step 1 is considered as aglobal error, which is defined as, ε global = E { ∆ u k ( θ ) } E { u k ( θ ) } = E { λ k ( θ ) } d Tk d kk (cid:80) i =1 k (cid:80) j =1 E { λ i ( θ ) λ j ( θ ) } d Ti d j = E { λ k ( θ ) } k (cid:80) i =1 E { λ i ( θ ) } (29)which measures the contribution of the k -th couple ( λ k ( θ ) , d k ) to the stochas-tic solution u ( θ ) and converges to the final solution when it achieves therequired precision. Further, the stop criterion for computing each couple( λ k ( θ ) , d k ) is considered as a local error and defined as, ε local = (cid:13)(cid:13)(cid:13) d ( j ) k − d ( j − k (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d ( j − k (cid:13)(cid:13)(cid:13) = (cid:13)(cid:13)(cid:13) d ( j ) k − d ( j − k (cid:13)(cid:13)(cid:13) (30)which measures the difference between d ( j ) k and d ( j − k and the calculation isstopped when d ( j ) k is almost the same as d ( j − k .15 . Applications The numerical implementation of the proposed method is illustrated withthe aid of three practical applications. The first application consists of anelectric pylon frame with stochastic material properties and a stochastic load.The second application is a roof truss under stochastic wind loads definedin low-dimensional and high-dimensional stochastic spaces, which is to il-lustrate the efficiency of applying the proposed method to high-dimensionalstochastic problems. The third example tests the ability of the proposedmethod for dealing with a large-scale engineering problem given by comput-ing the deformation of a tunnel under the action of self-weight. These threeexampls serve to verify the validity and accuracy of the proposed method anddemonstrate that there is the same solution construct for different problems.
In this problem, we consider a frame system as shown in Figure 1, whichis a electric pylon frame consisting of 91 elements with square cross-sections.A load P is applied vertically downward at the far right tip of the arm of thepylon. Clamped boundary conditions are applied at the base of the framemodel. Spatial nodes of the electric pylon frame model are defined in Table1. All elements of the electric pylon frame are constructed of 300M steel andhave identical cross-sectional areas. Deterministic material properties aregiven as, mass density ρ = 7 . , cross-sectional area A = 4cm , Young’smodulus E = 200GPa.The response of the electric pylon forced under a load P deeply depends onthese parameters. In order to better reflect the structural response influenced16
10 20 30 40 500102030405060 P Figure 1: 91-element electric pylon frameTable 1: Nodal definitions of the electric pylon frame node x y node x y node x y node x y EA = ( ξ ( θ ) + 0 . ξ ( θ )) EA, EI = ( ξ ( θ ) + 0 . ξ ( θ )) EI (31)and consider a stochastic load P ( θ ) as P ( θ ) = (1 + ξ ( θ ) + ξ ( θ )) P (32)where P = 1000N. Indepdent random variables { ξ i ( θ ) } i =1 in Equation (31)and (32) satisfylog { ξ i ( θ ) } i =1 ∼ N (0 , . , ξ ( θ ) , ξ ( θ ) ∼ N (0 , .
1) (33)Similar to the derivation of Equation (16), a stochastic finite elementequation for this problem can be obtained as, (cid:32) (cid:88) i =1 ξ i ( θ ) K i (cid:33) u ( θ ) = (cid:32) (cid:88) i =5 ξ i ( θ ) (cid:33) F (34)In order to solve Equation (34) by use of Algorithm 1, the convergencecriterias are set as ε = ε = 10 − , R = 1 × random samples, i.e. (cid:8) ξ i (cid:0) θ ( r ) (cid:1)(cid:9) × r =1 , i = 1 , · · · ,
2, are given in Equation (34) and { λ (0) k ( θ ( r ) ) } × r =1 are adopt in step 2 in Algorithm 1. In this example, only two retained termsin Equation (19), the displacement components { d i ( x, y ) } i =1 and probabil-ity density functions (PDFs) of corresponding random variables { λ i ( θ ) } i =1 shown in Figure 2, can achieve the required precision, which demonstratesthe high efficiency of the proposed method. It is seen from Figure 2 that, thesecond random variable λ ( θ ) is very small and makes almost no contribu-tions to the approximate solution u ( θ ). In practical, we determine whether18
20 400204060 0 20 4002040600 5 1000.10.20.30.4 -2 0 210 -5 Figure 2: Displacement components { d i } i =1 (top) and PDFs of corresponding randomvariables { λ i ( θ ) } i =1 (bottom) -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 00204060 Figure 3: Comparison of PDFs between the Monte Carlo simulation and the proposedmethod d and the random variable λ ( θ ) are necessary.Further we compare the proposed method with existing methods, includ-ing Monte Carlo simulation [8] and spectral stochastic finite element method(SSFEM) [21, 6]. Here Hermite Polynomial Chaos (PC) of 6 standard Gaus-sian random variables are adopted in the SSFEM, and the order of PC is setas p = 3 and p = 4. We test the computational efficiency of these methodsby use of a personal laptop (dual-core, Intel i7, 2.40GHz) and the computa-tional times of the proposed method, PC ( p = 3), PC ( p = 4) and 1 × standard Monte Carlo simulations are 3.4s, 71.9s, 474.9s and 1412.7s, re-spectively, which demonstrates the high efficiency of the proposed method.Based on above methods, the resulted approximate PDFs of the response ofthe far right tip of the arm of the electric pylon are seen from Figure 3. Theresult of the two-term approximation of the proposed method is in very goodaccordance with that from the Monte Carlo simulation, while the PC methodrequires fourth order ( p = 4) to achieve a similar accuracy. In addition, ourmethod is based on random samples, thus can avoid choosing the order p ofPC basis. We observe in practice that the number of random samples hasless influence on the computational cost. In the general case, the sample sizeis enough when it is sufficient to describe the statistical characteristics of therandom variables. In this example, we consider the stochastic response of a roof truss undera stochastic wind load acting vertically downward on the roof. The roof truss,as shown in Figure 4, includes 185 spatial nodes and 664 elements, where ma-20erial properties of all members are set as Young’s modulus E = 209GPa andcross-sectional areas A = 16cm . The stochastic wind load is a random field × × × . Figure 4: Model of the roof truss with the covariance function C F F ( x , y ; x , y ) = σ F e −| x − x | / l x −| y − y | / l y ,where the variance function σ F = 0 .
15, the correlation lengths l x = l y = 24,and it can be expanded by use of Karhunen-Lo`eve expansion Equation (3)with a M -term truncated as f ( x, y, θ ) = M (cid:88) i =0 ξ i ( θ ) f i ( x, y ) (35)where ξ ( θ ) ≡ f ( x, y ) = 10kN. { f i ( x, y ) } Mi =1 isobtained by solving Equation (3). Based on the expansion Equation (35) ofthe stochastic wind load, the following stochastic finite element equation isobtained, Ku ( θ ) = M (cid:88) i =0 ξ i ( θ ) F i (36)In this example, the initializations give the random samples (cid:8) ξ i (cid:0) θ ( r ) (cid:1)(cid:9) × r =1 , i =1 , · · · , M and the initial random variable samples { λ (0) k ( θ ( r ) ) } × r =1 , and set21 d d d d d d d ( a ). Displacement components { d i } i =1 λ -0.01 0 0.010100200300400 λ -4 -2 0 2 4 × -3 λ -2 0 2 × -3 λ -1 0 1 2 × -3 λ -10 -5 0 5 × -4 λ -4 -2 0 2 × -4 λ -2 0 2 × -4 λ ( b ). PDFs of { λ i ( θ ) } i =1 k -6 -4 -2 E rr o r ( c ). Iterative errors of k retained items Figure 5: Solutions of the couples { λ i ( θ ) , d i } i =1 and iterative errors of the solving process the convergence criterias as ε = ε = 10 − . We first consider a low-dimensional case by choosing M = 10. It is seen from Figure 5c that thedisplacement components { d i } and corresponding random variables { λ i ( θ ) } can be determined after 8 iterations, which demonstrates the fast conver-gence rate of the proposed method. Correspondingly, the number of couples( λ k ( θ ) , d k ) that constitute the stochastic response is adopted as k = 8. Asshown in Figure 5a and Figure 5b, with the increasing of the number of22ouples, the ranges of corresponding random variables are more closely ap-proaching to zero, indicating that the contribution of the higher order randomvariables to the approximate solution decays dramatically.For the maximum displacement of the whole roof truss, the resulted ap-proximate PDF compared with 1 × standard Monte Carlo simulations(MCS) is seen in Figure 6, which indicates that the result of eight-termapproximation is in very good accordance with that from the Monte Carlosimulation. According to our experience, further increasing the number ofcouples will not significantly improve the accuracy since the series in Equa-tion (19) has converged and thus the first few couples dominate the solutionof the problem. This example demonstrates the success of our proposed con-struct of the stochastic solution and Algorithm 1 for the solution of practicalproblems. -8 -7 -6 -5 -410 -3 Monte CarloProposed Method
Figure 6: Comparison of PDFs between theMCS and Algorithm 1
10 100 200 300 400 500 600 700 800 900 1000020406080100
Figure 7: Time costs of different stochasticdimensions M = 10 ∼ One of the main purposes the proposed method is to solve high-dimensionalstochastic problems. Here we introduce the high-dimensional stochastic prob-lems by choosing M = 100 ∼ This example is to compute the deformation of a tunnel under the actionof self-weight [26]. In order to reduce the size of the stochastic finite elementequation while ensuring the accuracy, triangle elements with gradients areused to generate a fine mesh for the tunnel structure and a coarse mesh forthe rock, totally including 2729 nodes and 5145 triangle elements, as shown inFigure 8. Material properties and thicknesses of all components are seen fromTable 2, here we consider the Young’s modulus of components as a randomfield with the mean value shown in Table 2 and the covariance function C EE ( x , y ; x , y ) = σ E e −| x − x | / l x −| y − y | / l y , where variance function σ E =0 .
1, correlation lengths l x = 10, l y = 20. Similar to Example 5.2, we modelthe Young’s modulus random field by use of Karhunen-Lo`eve expansion with10 terms, and derive a stochastic finite element equation.Given the random samples (cid:8) ξ i (cid:0) θ ( r ) (cid:1)(cid:9) × r =1 , i = 1 , · · · ,
10, the initial ran-dom variable samples { λ (0) k ( θ ( r ) ) } × r =1 and set the convergence criterias as ε = 10 − , ε = 10 − , displacements { d i } and corresponding random vari-ables { λ i ( θ ) } can be determined after 6 iterations, as shown in Figure 9 d ,which indicates the high efficiency of the proposed method. Figure 9( a –24 ock rock reinforcementconcrete spray (C25)backfillingconcrete (C10)concrete lining(C25)
15 59.40.953.5 10.5 r= Figure 8: Model of the tunnel (top) and the finite element mesh (bottom)Table 2: Descriptions of materials properties
Young’s modulus Poisson’s ratio mass density thickness(GPa) (kg/m ) (m)rock 2.0 0.25 2200rock reinforcement 2.6 0.20 2300 2.80concrete lining 28.5 0.20 2500 0.20backfilling concrete 18.5 0.20 2300 0.50concrete spray 28.5 0.20 2200 0.9525 .23e-021.56e-028.91e-032.24e-03-4.44e-03-1.11e-021.51e-023.00e-03-9.11e-03-2.12e-02-3.33e-02-4.54e-022.08e-021.60e-021.12e-026.44e-031.64e-03-3.16e-03 5.24e-03-1.90e-04-5.62e-03-1.11e-02-1.65e-02-2.19e-021.66e-028.59e-036.08e-04-7.37e-03-1.54e-02-2.33e-02 4.51e-022.73e-029.50e-03-8.31e-03-2.61e-02-4.39e-03 ( a ). Displacement { d i } i =1 in x direction -2.74e-20-8.27e-03-1.65e-02-2.48e-02-3.31e-02-4.13e-02 2.40e-021.33e-022.51e-03-8.25e-03-1.90e-02-2.98e-02 3.98e-022.18e-023.70e-03-1.44e-02-3.24e-02-5.05e-024.74e-022.95e-021.17e-02-6.07e-03-2.39e-02-4.17e-02 2.06e-02-1.25e-02-4.55e-02-7.86e-02-1.12e-01-1.45e-01 4.16e-021.87e-02-4.16e-03-2.70e-02-4.99e-02-7.28e-02 ( b ). Displacement { d i } i =1 in y direction λ -0.005 0 0.005 0.01 0.015050100150200 λ -0.006 -0.002 0 0.002 0.0060100200300400 λ -2 0 2 × -3 λ -5 0 5 × -4 λ -1 0 1 × -3 λ ( c ). PDFs of { λ i ( θ ) } i =1 -8 -6 -4 -2 ( d ). Iterative errors of k retained items Figure 9: Solutions of the couples { λ i ( θ ) , d i } i =1 and iterative errors of the solving process c ) shows the displacement components { d i } i =1 and PDFs of correspondingrandom variables { λ i ( θ ) } i =1 , where Figure 9 a and b are the displacementcomponents in the x direction (horizontal direction) and y direction (verticaldirection), respectively. Mean values and variances of the displacement areshown in Figure 10 a and b , and as a part of the whole displacement (shownin Figure 8 bottom), mean values and variances of the tunnel displacementare seen from Figure 10 a and b . Both tunnel displacements and rock dis-placements can be captured efficiently, which once again demonstrates theeffectiveness of the proposed method. Comparing with Figure 9( a – c ) and26igure 10 a , we observe that the first retained item, i.e. E { λ ( θ ) } d x and E { λ ( θ ) } d y , can roughly approximate the mean displacements in the x di-rection and y direction. For most cases, the mean displacement consideringuncertainties is very close to the displacement obtained from the determinis-tic case, but considering uncertainties can be better to reflect the variabilitiesof displacements, which is of great significance for structure design and eval-uation, such as reliability analysis and sensitivity analysis [27]. It is seenfrom Figure 10 b that, randomness in this example has more influence on thevariance of the tunnel displacement in the x direction and less influence onthat in the y direction, which provides a potential way for the design andevaluation of tunnel structures considering uncertainties.( a ). Means in the x and y direction ( b ). Variances in the x and y direction( a ). Means of the tunnel displacement ( b ). Variances of the tunnel displacement Figure 10: Means and variances of the displacement in the x and y direction
6. Conclusions
In this paper, we develop a method for solving stochastic finite elementequations and illustrate its accuracy and efficiency on three practical exam-ples. The proposed method sovles stochastic problems by use of a universal27olution construct and a dedicated iterative algorithm. It allows to solve high-dimensional stochastic problems with very low computational costs, whichhas been illustrated on numerical examples. Thus it appears as a powerfulway to avoid the Curse of Dimensionality. In addition, since the stochasticanalysis and deterministic analysis in the solving procedure are implementedin their individual spaces, the existing FEM and ODE codes can be readily in-corporated into the computatinoal procedure. In these senses, this method isparticulary appropriate for large-scale and high-dimensional stochastic prob-lems of practical interests and has great potential in uncertainty quantifi-cation of practical problems in science and engineering. In the follow-upresearch, it hopefully further applies the proposed method to a wider rangeof uncertainty quantification, such as reliability analysis, sensitivity analysis,etc.
Acknowledgments
This research was supported by the Research Foundation of Harbin In-stitute of Technology and the National Natural Science Foundation of China(Project 11972009). These supports are gratefully acknowledged.