Generalized stochastic cell-based smoothed finite element method (GS_CS-FEM) for solid mechanics

Finite Elements in Analysis and Design 63 (2013) 51–61

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Generalized stochastic cell-based smoothed finite element method (GS_CS-FEM) for solid mechanics G.R. Liu a, W. Zeng a,n, H. Nguyen-Xuan b,c a

School of Aerospace Systems, University of Cincinnati, 2851 Woodside Dr, Cincinnati, OH 45221, USA Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 700000, Viet Nam c Department of Mechanics, Faculty of Mathematics & Computer Science, University of Science, VNU-HCMC, Nguyen Van Cu Street, District 5, Ho Chi Minh City 700000, Viet Nam b

a r t i c l e i n f o

abstract

Article history: Received 30 December 2011 Received in revised form 24 August 2012 Accepted 27 August 2012 Available online 29 September 2012

The smoothed finite element method (S-FEM) was recently proposed to inject softening effects into and improve the accuracy of the standard FEM. In the S-FEM, the system stiffness matrix is obtained using strain smoothing technique over the smoothing domains associated with cells, nodes, edges or faces to establish models of desired properties. Stochastic FEM is regarded as an extension of the classical deterministic FEM to deal with the randomness of properties of input parameters in solid mechanics problems. In this paper, the cell-based S-FEM (or CS-FEM) is extended for stochastic analysis based on the generalized stochastic perturbation technique. Numerical examples are presented and obtained results are compared with the solution of Monte Carlo simulation. It is found that the present GS_CS-FEM method can improve the solution accuracy significantly for stochastic problems, in terms of the estimated means, variances and other probabilistic characteristics. & 2012 Elsevier B.V. All rights reserved.

Keywords: Smoothed finite element method (S-FEM) Stochastic FEM General perturbation technique Monte Carlo simulation

1. Introduction For several decades, the finite element method (FEM) has become one of the most popular numerical tools in solving practical problems in aeronautical, mechanical and civil engineering. Due to the complexity nature of the problems, lower order FEM is widely preferred. However, FEM using lower order elements has some inherent limitations leading to deficiency in many applications. For example, it exhibits an overestimation of stiffness matrix and underestimation of the internal strain energy. In addition, the element shape cannot be distorted too much due to the use of mapping in FEM [1]. To overcome these issues, the smoothed finite element method (S-FEM) was proposed recently by combining the strain smoothing technique often used in mesh-free methods [2] with the standard finite element method (FEM) [3,4]. The S-FEM may be formulated using the weakened weak (W2) formulation based on the G space theory [5]. It has been proven that a W2 model behaves ‘‘softer’’ than the corresponding standard Galerkin weak form model. Furthermore, upper bound solutions (to the exact solution for force-driving problems) can be also derived using W2 models with sufficient softness, such as when the NS-FEM is used. Therefore, we can now bound the solution from both sides using a

n

Corresponding author. Tel.: þ41 513 477 5889. E-mail address: [email protected] (W. Zeng).

0168-874X/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2012.08.007

standard weak form together with a proper W2 form [6], if so desired. In a more general form of meshfree setting, the smoothed point interpolation methods (S-PIMs) have also developed as typical W2 models. The S-PIM can be node-based (known as NS-PIM or LC-PIM) [7], edge-based (ES-PIM) [8], and cell-based (CS-PIM) [9]. The S-FEM can be, on the other hand, regarded as the linear form of S-PIM, and it is much simpler than S-PIM but reserving its major properties [10]. The S-FEM models can use the same FEM mesh, and work particularly well for lower order elements. Four different types of smoothing domains created leading to cell-based S-FEM (CS-FEM), node-based S-FEM (NSFEM), edge-based S-FEM (ES-FEM), and face-based S-FEM (FSFEM). Each of the four S-FEM models will have different features and properties [11], but they all have the following general properties: (1) There is no need to compute the derivatives of shape functions to form the stiffness matrix, because of the use of smoothing operation and field gradients are computed using shape functions directly; (2) No coordinate transformation or mapping is involved in S-FEM, and hence no restriction on the elements shape and hence mesh can be heavily distorted; (3) Many existing algorithms of FEM can be easily modified and applied to S-FEM with very little change; (4) The approximated accuracy and convergence rate can be improved with little increase in computational cost leading much higher computational efficiency, compared to the standard FEM [3]. Randomness of parameters is a natural characteristic in many engineering systems, and should be properly dealt with in

52

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

computational modeling and simulation. Such randomness creates uncertainty in the outputs of the numerical models. Deterministic FEM fails to consider the uncertainty, which limits its application [12]. From early 1980s, the concept of the stochastic FEM has been established, by combining the reliability evaluation approaches into FEM technology. During the last several decades, the stochastic FEM has become a powerful tool in computational stochastic mechanics, in dealing with largescale realistic engineering problems. Concerning the state-of-theart review of the past and recent developments in stochastic FEM area, we can refer to Ref. [13]. A nth order generalized stochastic finite element (GS-FEM) based on perturbation technique was proposed by Kamin´ski [14]. In large variations of random input parameters, the second order perturbation expansion was found ineffective. In Ref. [14], a 1D linear elastostatics problem with a single random variable showed that the accuracy of the expected values and variances can be improved using the nth perturbation order technique. The nth order approach also means it is possible that the probabilistic moments of the solution can be computed with a priori given accuracy [15]. In this paper, we establish a generalized stochastic CS-FEM (GS_CS-FEM), which combines the generalized nth order stochastic perturbation technique with the cell-based smoothed finite element approach for 2D solid mechanics problems. Two numerical examples are presented to demonstrate the effectiveness of the present method, in comparison with the Monte Carlo simulation approach. The paper is organized as follows. In Sections 2 and 3, the probability theory and the CS-FEM are briefly presented, respectively. The formulation of GS_CS-FEM using the generalized nth order perturbation stochastic approach is detailed in Section 4. Section 5 shows two numerical examples. Finally, some concluding remarks are made in Section 6.

2. Briefing of probability theory For a given set of random fields b(x) and its probability density function (PDF)pi(br), r ¼ 1,2,. . .,R, i ¼1,2, the first two probabilistic moments for the random fields br(xk) are defined as [16,17] Z þ1 0 br p1 ðbr Þdbr ð1Þ E½br  ¼ br ¼ 1

Covðbr ,bs Þ ¼ Srs ¼

Z

þ1 1

Z

þ1 1

0

0

ðbr br Þðbs bs Þp2 ðbr ,bs Þdbr dbs

ð2Þ

0

where br represents the first probabilistic moment of the random variable, Cov(br,bs) or Srs represents the covariance, p1(br) and p2(br,bs) denote the PDF and the joint PDF, respectively. For a real, single-valued continuous function of a random variable b, the expectation becomes the expected value, or mean of the variable b. Z þ1 0 E½b ¼ b ¼ bpðbÞdb: ð3Þ 1

The variance of b, expressed as Var(b), or s2b is defined as the mean square value of b about the mean Z þ1 Z þ1 0 0 ðbb Þ2 pðbÞdb: ð4Þ VarðbÞ ¼ E½ðbb Þ2  ¼ 1

1

Then the standard deviation (denoted by sb) and the coefficient of variation (denoted by g or COV) of a random variable b can be defined as: " #1=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi VarðbÞ sb ¼ VarðbÞ and g ¼ : ð5Þ ðE½bÞ2

Based on stochastic perturbation technique, if a small perturbation parameter e related to spatial expectations is adopted, the nth order truncated Taylor series expansion of the limit state function of a structural reliability analysis [14] can be expressed as: 0

f ðbÞ ¼ f ðbÞ þ

1 X 1 n ðnÞ e f ðbÞðDbÞn ffif 0 ðbÞ þ ef ,b ðbÞDb n! n¼1

þ  þ

1 n ðnÞ e f ðbÞðDbÞn n!

ð6Þ

where f(n)(b)¼(@nf(b)/@bn), eDb¼ e(b b0)and e2(Db)2 ¼ e2(b b0)2 are the first and second variation of b concerning the corresponding expected value b0, respectively. The nth order variation can be similarly written as:

en ðDbÞn ¼ en ðbb0 Þn

ð7Þ

0

where (U) denotes the mean value of the function value (U) taken at the expected value b0, and (U),b, (U),bb represents the first and second order partial derivative with respect to b at the point b0. Replacing br in Eq. (1) by f(b) in Eq. (6), we obtain the expected value of any limit state function f(b) with a specified small perturbation parameter e by expansion of Taylor series [14]: Z þ1   E f ðbÞ; b ¼ f ðbÞpðbÞdb 1 ! Z þ1 1 X 1 n ðnÞ 0 n ¼ f ðbÞ þ e f ðbÞðDbÞ pðbÞdb n! 1 n¼1 8 ! M X > R þ1 1 2n ð2nÞ > 0 2n > e f ðbÞðDbÞ pðbÞdb > 1 f ðbÞ þ > > ð2nÞ! > n¼1 > > > < for symmetric distribution functions ! ð8Þ ffi N X > 1 n ðnÞ > R þ1 0 n > > ð Þdb f ðbÞ þ e f ðbÞð D bÞ p b > 1 > ðnÞ! > > n¼1 > > : for asymmetric distribution functions Note that the approximation of expected values or variances satisfies a given priori precision with an admissible error via a proper selection of number of terms, M or N. For a small variation of the random variable with the symmetric PDF around its mean value, the expected value for the input random variable with symmetric probability density function in the second order perturbation approach can be expressed as [15]: Z þ1   E f ðbÞ; b ¼ f ðbÞpðbÞdb  Z1 þ1  1 ,bb 0 ,b f ðbÞ þ f ðbÞeDb þ f ðbÞe2 ðDbÞ2 pðbÞdb ffi 2 1 1 0 ,b ,bb ¼ f ðbÞ þ 0  f ðbÞe þ f ðbÞe2 Sbb ð9Þ 2 where Sbb stands for the second-order central moments, which is unique for a given random variable b. When the input random variable with the symmetric PDF was in a large variation, the following extension with a perturbation parameter e can be preferably adopted [15–19]:   1 ,bb 1 ,bbbb 0 ðbÞe2 m4 ðbÞ E f ðbÞ; b ¼ f ðbÞ þ f ðbÞe2 m2 ðbÞ þ f 2 4! 1 ,bbbbbb ðbÞe6 m6 ðbÞ þ    þ f 6! 1 @2 f 2 1 @4 f 2 0 e m2 ðbÞ þ e m4 ðbÞ ffif ðbÞ þ 2 @b2 4! @b4 þ

1 @6 f 6 e m6 ðbÞ þ    6! @b6

ð10Þ

where mn(b) denotes the nth order central probabilistic moment of b and the odd order terms vanish for a symmetry PDF (such as

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

Gaussian random distribution), and where higher than sixthorder terms are neglected. The sixth-order truncated expansion for a variance can be analogously employed [14]: Z þ1    1 ,bb 1 ,bbb 0 ,b Var f ðbÞ ¼ f þ f eDbþ f ðeDbÞ2 þ f ðeDbÞ3 2 3! 1 2 1 ,bbbb 1 ,bbbbb þ f ðeDbÞ4 þ f ðeDbÞ3 Eðf Þ pðbÞdb 4! 5! Z þ 1  2 ,b 2 f ðeDbÞ pðf ðbÞÞdb ffi 1

þ

1 4

Z

þ1

 2 ,bb f ðeDbÞ4 pðbÞdb

1

ð11Þ when the state function does not possess a symmetric feature, the odd orders of probabilistic moments should be nonzero. Even so, the procedure would be implemented in the same way.

3. Briefing of the cell-based smoothed finite element method 3.1. Governing equation and weak form Consider now a 2D static elasticity problem, the governing equation in the problem domain O bounded by G can be expressed in terms of stresses as: ð12Þ

which subjects to the boundary conditions ui ¼ u^ i on Gu and

sij nj ¼ t^i on Gt, where sij denotes the component of stress tensor and bi denotes the component of body force; nj is the unit vector normal to the boundary G. The Galerkin weak form of the FEM can be employed as Z Z drs ðuÞij Dijkl rs ðuÞkl dO t^i dui dG ¼ 0 ð13Þ O

Gt

where Dijkl denotes the elasticity tensor components. 3.2. Brief on the FEM In the FEM, if the problem field is divided into n elements through spatial discretization procedure, the displacement and the compatible strains e ¼ rSu within any element can be expressed as [20]: n X

h

u ðxÞ ¼

NI ðxÞdI

eh ðxC Þ ¼

BI ðxC ÞdI

2

N I0 x

3

0

6 BI ¼ 4 0 N I0 y

2

bIx

0

6 N I0 y 7 5¼40 N I0 x

3

bIy 7 5

bIy

ð16Þ

bIx

Then the element stiffness matrix for elasticity can be expressed as: Z BT DBdO ð17Þ Ke ¼ where Oe denotes the element domain. If we assemble each of element stiffness matrix Ke into a global stiffness matrix K, the global equilibrium equations of FEM yields Kd ¼ f

ð18Þ

where d stands for the global nodal displacement vector and f is the nodal force vector given from Z Z f I ¼ NTI bdO þ NTI tdG: ð19Þ O

G

3.3. Brief on the cell-based smoothed finite element method (CS-FEM) Discretization of the problem domain into elements for the CS-FEM is similar to the conventional FEM. However, the shape of elements could be much more versatile. The element can be general n-sided polygon [21] or other arbitrary (even concave) shapes, including triangular elements (T3) or quadrilateral elements (Q4) for 2D problems and tetrahedron elements (T4) for 3D problems [11,21]. When quadrilateral elements (Q4) (see, Fig. 1) are used, the meshing may be performed in exactly the same way as in the FEM. In the CS-FEM, the basic elements would be further subdivided into several number (nSC) of smoothing cells, where nSC depends on the stability condition [21] and accuracy requirements. The strain smoothing operation is performed over each smoothing cell, and the contribution to the stiffness matrix is evaluated. The element stiffness matrix is assembled using the stiffness matrix of each cell. A smoothing operation is performed using R

e uh ðx Þ ¼ r C

OC ru

R

ð15Þ

h

ðxÞdO

OC d O

¼

1 AC

Z

ruh ðxÞdO

ð20Þ

OC

R where OC denotes the smoothing cell, and AC ¼ OC dO. Note that O ¼ O1 [ O2 [    OSC and O1 \ O2 \    OSC ¼ |. Similar to the FEM, the approximation of displacement can be expressed as in Eq. (14). Using integration by parts and using the divergence theory, the smoothed gradient of displacement field is obtained as 1 AC

e uh ðx Þ ¼ r C

ð14Þ

I¼1 n X

matrix. For a 2D elastic problem, it would be given as:

Oe

Z 1 þ 1 ,b ,bbb 2 f f ðeDbÞ4 pðf ðbÞÞdb 3! 1  2 Z þ 1  2 1 ,bbb f ðeDbÞ6 pðbÞdb þ 3! 1 Z 1 þ 1 ,bb ,bbbb f f ðeDbÞ6 pðf ðbÞÞdb þ 4! 1 Z 1 þ 1 ,b ,bbbbb þ2  f f ðeDbÞ6 pðbÞdb 5! 1    2 1 2 ,b ¼ f e2 m2 ðbÞ þ ðf ,bb Þ2 þ f ,b f ,bbb e4 m4 ðbÞ 4 3! !  2 1 1 ,bb ,bbbb 2 ,b ,bbbbb 6 ,bbb 2 ðf Þ þ f f þ f f e m6 ðbÞ þ 3! 4! 5!

sij,j þ bi ¼ 0 in O

53

Z

uh ðxÞnðxÞdG ¼ GC

n 1 X uI AC I

Z

N I ðxÞnðxÞdG

ð21Þ

GC

where GC represents the boundary of the smoothing cell OC, n(x) is the vector matrix normal to the boundary GC. Then the smoothed strains are rewritten as:

I

where dI ¼ ½uI vI T is the nodal displacement vector. N is the shape function vector and B is the standard displacement gradient

eeh ðxC Þ ¼

n X I

e I ðxC ÞdI B

ð22Þ

54

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

Fig. 1. Division of a quadrilateral element into smoothing cells (SCs) in CS-FEM by connecting the mid-segment-points of opposite segments of smoothing domains [22]: (a) SC¼ 1; (b) SC ¼2; (c) SC ¼3; (d) 4 SC ¼4; (e) 8 SC¼ 8; and (f) 16 SC¼ 16.

where dI ¼ ½uI vI T is the nodal displacement vector. The smoothed e I ðxC Þ is calculated by strain matrix B 2

e ðx Þ b I1 C

6 e I ðxC Þ ¼ 6 0 B 4 e ðx Þ b I2 C

3 0 7 e ðx Þ 7: b I2 C 5 e ðx Þ b I1

ð23Þ

C

Obviously, the relationship between the smoothed strain e I ðxC Þ and standard strain matrix in S-FEM can be matrix B obtained as follows: Z Z Z 1 1 1 eI ¼ NI ðxÞnðxÞdG ¼ rs NI ðxÞdO ¼ BI ðxÞdO ð24Þ B A C OC GC AC OC AC e I is the averaged value of the standard Eq. (24) shows that B strain matrix BI(x) over the smoothing cell OC. Because a linear compatible displacement field is used along boundary GC, one Gauss point is sufficient to evaluate exactly the line integration on each segment of boundary GCi of OC. Therefore, Eq. (24) can be transformed to its algebraic form e ðx Þ ¼ b Ik C

M X



C C NI xGP nik li i

ð25Þ

where k¼1,2; M is the total number of the edge number of the domain OC, and xGP is the midpoint, i.e., Gaussian point, of i boundary segment of GCi , whose outward unit normal and length C are represented as nCi and li . The ‘‘smoothed’’ element stiffness matrix is assembled from e e , that is calculated using stiffness of each smoothing cell, K ee ¼ K

e T DB e C AC : B C

n X eh ðxC Þ @n e e I ðxC Þ @dI ¼ B n n @b @b I

ð28Þ

In the discrete global CS-FEM equilibrium Eq. (27), if some random quantities are integrated into the smoothed stiffness e and the force vector f, the following hierarchical matrix K GS_CS-FEM equations for elastostatics are obtained by employing the standard procedure:

e 0 d0 ¼ f 0 K

ð29Þ

– First-order (e1 terms) e 0 d,b þ K e ,b d0 ¼ f ,b K

ð30Þ

– Second-order (e2 terms) e 0 d,bb þ 2K e ,b d,b þ K e ,bb d0 ¼ f ,bb K

ð31Þ

ð26Þ

C¼1

e assembled from each The smoothed global stiffness matrix Kis e e . Then the discrete global CS-FEM of element stiffness matrix K equilibrium equations in displacement format can be written as: e ¼f Kd

The general perturbation approach is now applied based on smoothed FEM settings. For simplicity, we call the present methods as generalized stochastic cell-based smoothed finite element method (GS_CS-FEM), an analogous of the GS-FEM [23]. e I ðxC Þ that is stochastiConsider the smoothed strain matrix B cally independent from random variables. In general, the partial eh ðxC Þ with respect to a random derivative of the smoothed strain e variable b can then be derived from Eq. (22):

– Zeroth-order (e0 terms)

i¼1

nSC X

4. The formulation of the generalized stochastic CS-FEM

ð27Þ

where f is the general nodal force vector expressed as the same form in Eq. (19).

– Third-order (e3 terms) e 0 d,bbb þ3K e ,b d,bb þ 3K e ,bb d,b þ K e ,bbb d0 ¼ f ,bbb K

ð32Þ

– Nth-order (eN terms, notation using Pascal’s rule)  N  X N e ðkÞ ðNkÞ ðNÞ ¼f K d k k¼0

ð33Þ

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

55

zfflfflffl}|fflfflffl{k where the symbol (U)(k) denotes ðUÞ, bb. . .b , what means kth order partial derivative with respect to b evaluated at (U)0. e is Recall now that the smoothed global stiffness matrix K e e from Eq. (26), i.e., assembled from element stiffness matrix K e¼ K

nSC ne X X

e C AC e T DB B C

ð34Þ

e¼1 C

where ne is the total number of elements in domain O. Regarding Young’s modulus E as the random variable to be considered, the kth order derivatives of smoothed global stiffness matrix with respect to E are then written as: nSC ne X e X @K e C AC e T @D B ¼ B C @E @E e¼1C ¼1

Fig. 2. A cantilever subjected to a parabolic traction at the free end.

ð35Þ

k e when kZ2, the result @k K=@E ¼ 0 can be derived for elastic problems transparently. Noted that @k f=@E k ¼ 0 for any kZ1 as Young modulus E has no effect on the force vector f. Then the corresponding GS_CS-FEM equilibrium equations can be simplified into

Fig. 3. Domain discretization using 4-node quadrilateral elements of the cantilever (mesh 32  8).

0

– Zeroth-order (e terms) e 0 d0 ¼ f 0 K

ð36Þ

– First-order (e1 terms) e ,E d0 e 0 d,E ¼ K K

ð37Þ

– Second-order (e2 terms) e ,E d,E e 0 d,EE ¼ 2K K

ð38Þ

– Third-order (e3 terms) e 0 d,EEE ¼ 3K e ,E d,EE K

ð39Þ

– Nth-order (eN terms) N N1 d e ,E @ e 0 @ d ¼ N K K @E N @E N1

ð40Þ

A recursive procedure can be implemented to acquire the N thorder solution from the above equation series.

5. Numerical illustrations In order to analysis the properties of the GS_CS-FEM, two numerical examples will be presented to study the probabilistic output moments. The first example is a simple rectangular cantilever subjected to a parabolic traction at the free end and a plane stress condition is assumed. For the second problem, the infinite plate with a circular hole subjected to unidirectional tension is under a plane strain condition. 5.1. Cantilever beam subjected to a tip load A rectangular cantilever beam with length L and height H, which is subjected to a parabolic traction at the free end, is shown in Fig. 2. The analytical solution of displacements is [24]: " !# P H2 ux ¼ ð6L3xÞxy þ ð2 þ nÞy y2  6EI 4

Fig. 4. Comparison of the relative error in displacement v between SFEM and FEM.

" # P H2 x 2 2 uy ¼  þ ð3LxÞx þ3ny ðLxÞ ð4þ 5nÞ 6EI 4

ð41Þ

where the moment of inertia I of the beam is given by I¼H3/12. The Young’s modulus E ¼ 3:0  107 kPa(we adopt the script letter ‘‘E’’ to distinguish it from the notation of the expectation ‘‘E’’), Poisson ratio n ¼0.3, L¼12 m, H¼4 m, and P¼1000 N. The perturbation parameter e is chosen as an interval eA[0.8,1.2], and the input coefficient of variation of the randomized modulus is set as gðEÞ A ½0:0,0:3. In Fig. 3, the domain is discretized by 4-node quadrilateral elements and these elements can be further divided into different SCs, SC ¼1, 2, 3, 4, 8, 16, as shown in Fig. 1. The relative errors of deflection along centerline between FEM/CS-FEM (mesh 32  8) and analytical solution are demonstrated in Fig. 4. It is seen the deflection v computed by CS-FEM can be more accurate than the result of FEM using 4 Gauss points for full integration. And when the number of SC of an element increases from 4 to 16, the error of deflection will reduce gradually and obviously. To study the

56

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

convergence rate of the present method, the two norms called displacement norm and energy norm can be provided as follows: P h ndof 9ui ui 9 Errord ¼ P 9u 9 i ndof Z

1=2 1 ðeh eÞT Dðeh eÞ ð42Þ Errore ¼ 2LD O In Fig. 5, the depicted convergence rate shows CS-FEM can give almost comparable convergence rates compared to FEM both in displacement and energy norms, i.e., CS-FEM preserves the full superconvergence feature similar to FEM. Furthermore, the error in energy norm for the CS-FEM is always smaller than that of FEM, as plotted in Fig. 5(b). If the input quantity E consists of Gaussian random distribution, then all the central probabilistic moments can be obtained from the formulas [14]

m2m þ 1 ðE Þ ¼ 0,

m2m ðE Þ ¼ 1U3U5    ð2m1Þs2m ðEÞ

¼ ð2m1Þ! ðVarðEÞÞm

ð43Þ

Obviously, the maximum vertical deflection v in the centerline is at the center point of the free end. Applying Eqs. (10) and (11), the expectation and variance for this maximum vertical deflection

−2.5

−3

FEM−Q4−4GP(r=1.98) CSFEM −4SC(r=1.99) CSFEM −8SC(r=1.99) CSFEM−16SC(r=1.99)

vm should be 

 1 @2 vm 2 1 @4 vm 4 E vm e, gðE Þ ffi v0m ðE Þ þ e ðgEÞ2 þ e ð3gEÞ4 2 @E 2 4! @E 4

1 @6 vm 6 þ e 3  5gE 6 þ    6! @E 6

ð44Þ

    

 @vm 2 2 1 @vm 2 e ðgEÞ2 þ Var vm e, gðE Þ ¼ 4 @E @E   3 !! 2 @vm @ vm e4 ð3gEÞ4 þ 3! @E @E 3 0 ! !  2 3 !2 1 @ vm 1 @2 vm @4 vm þ@ þ 3! 4! @E 2 @E 3 @E 4 !!   5

2 @vm @ vm e6 3  5gE 6 þ    þ 5! @E @E 5

ð45Þ

The expected values of the maximum vertical deflection in centerline for 2nd, 4th, 6th, 8th, 10th order stochastic smoothed finite element approximation are provided in Fig. 6. The corresponding standard deviations are collected in Fig. 7, and the output coefficients of variation are presented in Fig. 8. The expected values, standard deviations and the output coefficients of variation are shown with respect to (w.r.t.) only perturbation parameter e (while the input coefficient of variation g is set as 0.10 and 0.25, separately) in Figs. 9–14.

Expectations: 2,4,6,8,10th orders

0.011

−3.5

−4

10

0.0105

Expectation

log10 (Errord)

0.0115

8 6

0.01

4

0.0095 2

−4.5

0.009 0.0085 0.8

−5 −0.5

0

0.9

0.5

ε

log10 (h)

6

−1

5

−1.2

4

−1.4

0.15

0.2

0.25

0.3

γ

2

−1.8

1

−2

0 0.8

0

0.5

log10 (h)

x 10−3

Standard deviations: 2,4,6th orders 6 4

3

−1.6

−2.2 −0.5

0.05

0.1

Fig. 6. Expected values for 2, 4, 6, 8, 10th orders. FEM−Q4−4GP(r=1.48) CSFEM −4SD(r=1.50) CSFEM −8SD(r=1.51) CSFEM−16SD(r=1.54)

stdev

log10 (Errore)

−0.8

1.1 0

−0.4 −0.6

1

2

0.9

ε

1 1.1 1.2

Fig. 5. Comparison of convergence rate between SFEM and FEM: (a) displacement norm; (b) energy norm.

0

0.05

0.1

0.15 γ

0.2

Fig. 7. Standard deviations for 2, 4, 6th orders.

0.25

0.3

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

Coefficients of variation: 2,4,6 orders 0.7

11.5

0.6

11

4

0.3 0.2

2

0.1 0 0.8 0.9

ε

1 1.1 0

0.05

0.1

0.15 γ

0.25

0.2

8.5 8

5

Standard deviations

E 8.98 8.97 8.96 8.95 0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

Fig. 9. Expected values; 2, 4, 6, 8, 10th orders; g: 0.10.

x

10.2

10−3

Expected values: 2,4,6,8,10th orders, γ =0.25

9.8

9.6

9.4

0.9

1 ε

1.05

1.1

1.15

1.2

2nd order 4th order 6th order

3.5 3 2.5 2

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

Fig. 12. Standard deviations; 2, 4, 6th orders; g:0.25.

2nd order 4th order 6th order 8th order 10th order

0.85

0.95

4

1.5 0.8

10

9.2 0.8

0.9

Standard deviations: 2,4,6th orders, γ =0.25

x 10−3

4.5

0.8

0.85

Fig. 11. Standard deviations; 2, 4, 6th orders; g:0.10.

8.99

E

9

6.5 0.8

9

10.4

9.5

7

2nd order 4th order 6th order 8th order 10th order

9.01

10

Expected values: 2,4,6,8,10th orders, γ =0.1

x 10−3

9.02

2nd order 4th order 6th order

7.5

0.3

Fig. 8. Coefficients of variation for 2, 4, 6th orders.

9.03

Standard deviations: 2,4,6th orders, γ=0.1

x 10−4

10.5

6

0.4

Standard deviations

coef ε and γ

0.5

57

0.95

1 ε

1.05

1.1

1.15

1.2

Fig. 10. Expected values; 2, 4, 6, 8, 10th orders; g: 0.25.

From Figs. 6–8, it is observed that the expected values and standard deviations increase nonlinearly along with the increase of the perturbation parameter e or the input coefficient of

variation g. Besides, it is apparent that the convergence of the GS_CS-FEM depends on the input g. The 2nd Taylor expansion perturbation is suitable for the situation when the input gðEÞ is no more than 0.10, and otherwise, it requires a higher perturbation order approximation. In addition, note that the input coefficient of variation plays more crucial impact on these probabilistic characteristics than the perturbation parameter, especially as g r0.10. Comparing the expected values, standard deviations and the output coefficients of variation in case of g set as 0.10 and 0.25 given in Figs. 9–14, we can easily observe that for a smaller value of input coefficient of variation, such as 0.10, the 4th order perturbation approximation can satisfy the accuracy/convergence requirement. However, even 10th order perturbation is not enough for the case using g ¼0.25, especially for standard deviations and the coefficients of variation. So the 2nd expansion perturbation technique performs quite well in precision when the input g is less than 0.10 and can be effective for 0.10r g r0.15, but may lead to low accuracy for larger g. In Figs. 15 and 16, the expected values for 2, 4, 6, 8, 10th order GS_CS-FEM approximation are compared to the results of Monte

58

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

Coefficients of variation: 2,4,6th orders, γ =0.1

0.13

4.5 2nd order 4th order 6th order

0.12

x 10−3

2nd order 4th order 6th order MC simulation

4

Standard deviations

3.5

COV

0.11

0.1

0.09

Standard deviations: 2,4,6th orders, ε=1

3 2.5 2 1.5 1

0.08 0.5 0

0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

Coefficients of variation: 2,4,6th orders, γ =0.25

0.4

COV

0.35 0.3

0.2

0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

Fig. 14. Coefficients of variation; 2, 4, 6th orders; g:0.25.

x 10−3

Expected values: 2,4,6,8,10th orders 2nd order 4th order 6th order 8th order 10th order MC simulation

10.2 10

E

9.8 9.6 9.4 9.2 9 8.8

0.15 γ

0.2

0.25

0.3

5.2. Infinite plate with circular hole

0.25

10.4

0.1

Carlo simulation (MCs). It is evident that using higher order perturbations in GS_CS-FEM provides higher accuracy of solution compared to the MCs results. Referring to computational efficiency, the CPU time for the GS_CS-FEM is only 15.5 s, while it requires more than thousands times of time for MCs calculation in a reasonable accuracy as a reference solution (even only 2000 times picked as the number of analysis for MCs, more than 7760 s required). If other randomized variables are considered, one can also apply similar implementations using the GS_CS-FEM approach.

2nd order 4th order 6th order

0.45

0.05

Fig. 16. Standard deviations; GS_CS-FEM vs. MCs.

Fig. 13. Coefficients of variation; 2, 4, 6th orders; g:0.10.

0.5

0

0

0.05

0.1

0.15 γ

0.2

Fig. 15. Expected values; GS_CS-FEM vs. MCs.

0.25

0.3

Fig. 17 represents a plate with a central circular hole subjected to a unidirectional tensile load of s ¼ 1.0 N/m at infinity in the x direction. Since the stress concentration around the hole is highly localized and decays very rapidly, essentially disappearing when the distance to the center is greater than 5a, only a finite plate with l¼ 5a need to be modeled. Due to its symmetry, only the upper right quadrant of the plate is chosen and discretized as 144 elements. The case is considered as a plane strain problem and Young’s modulus E ¼ 1:0  103 N/m2, Poisson ratio n ¼0.3. The inner edge of the hole is traction free and symmetric conditions are set on the left and bottomed edges. On the right (x ¼5.0 m) and top (y¼5.0 m) edges, traction boundary conditions are imposed according to the exact solution [24]. The analytical solution for displacement components is [3,24]

a r a a3 d1 ¼ þ 2 ð1 þ kÞcos y þ cos 3yÞ2 3 cos 3y 0 8m a r r

 a r a a3 d2 ¼ ðk3Þ2 ðk1Þ sin y þsin 3yÞ2 3 sin 3y ð46Þ 0 8m a r r where m0 ¼ ðE=2ð1þ nÞÞ and k ¼3 4n for plane strain cases, (r,y) are the polar coordinates and y is measured counterclockwise from the positive x-axis. In order to avoid the redundancy, it is unnecessary to repeat the plots for convergence rate of displacement and energy, as these performances are quite similar to the cantilever case. Fig. 18 describes 4-node quadrilateral elements mesh used to discretize the problem domain. Similar to plane stress example, the modulus e is adopted as the input random variable, the perturbation parameter e is chosen as the interval of [0.8,1.2], and the input coefficient of

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

59

Fig. 17. Infinite plate with a circular hole subjected to x-directional tension and a symmetric geometry.

x 10−3

Standard deviations: 2,4,6th orders

2

6

stdev

1.5

4

1

0.5

2

0 0.8 0.9

ε

1 1.1 1.2

0

0.05

0.1

0.15 γ

0.2

0.25

0.3

Fig. 20. Standard deviations for 2, 4, 6th orders.

Fig. 18. Domain discretization using 4-node quadrilateral elements of the infinite plate with a circular hole.

Coefficients of variation: 2,4,6 orders

0.7 0.6

x 10−3

Expectations: 2,4,6,8,10th orders

0.5 coef ε and γ

3.5 3.4 10

Expectation

3.3 3.2 3.1 3

0.2

6

0.1

2

0 0.8

2.9 2

2.8

4

0.3

8

4

6

0.4

0.9

ε

2.7 0.8

1 1.1 0

0.9

ε

1 1.1 1.2

0

0.05

0.1

0.15 γ

0.2

0.25

0.3

Fig. 19. Expected values for 2, 4, 6, 8, 10th orders.

variation of the randomized modulus is set as g(E)A[0.0,0.3]. The horizontal displacement at point (1,0) is adopted as the output probabilistic variable. The expected values, standard deviations and output coefficients of variation are collected, respectively in Figs. 19–21. The expected values, standard deviations and output coefficients of variation are plotted w.r.t. only perturbation parameter e (while the input g is set as 0.10 and 0.25, separately) in Figs. 22–27. In Figs. 28 and 29, the expected values and standard

0.05

0.1

0.15 γ

0.2

0.25

0.3

Fig. 21. Coefficients of variation for 2, 4, 6th orders.

deviations obtained from both GS_CS-FEM are compared to the corresponding outcome from MCs. It is shown again those using higher order perturbations in GS_CS-FEM yield higher accuracy approximation since the consequences may match MCs results better. Several other properties of GS_CS-FEM in this problem are also quite similar to those in the previous problem.

6. Conclusion In this paper, a generalized stochastic cell-based smoothed finite element method (GS_CS-FEM) is proposed. The method is

60

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

2.76

Expected values: 2,4,6,8,10th orders, γ =0.1

x 10−3

13

2nd order 4th order 6th order 8th order 10th order

2.755

2nd order 4th order 6th order

12 11 Standard deviations

2.75

Standard deviations: 2,4,6th orders, γ =0.25

x 10−4

E

2.745 2.74 2.735

10 9 8 7

2.73

6 2.725 0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

5 0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

Fig. 22. Expected values; 2, 4, 6, 8, 10th orders; g:0.10. Fig. 25. Standard deviations; 2, 4, 6th orders; g:0.25.

3.1

Expected values: 2,4,6,8,10th orders, γ =0.25

x 10−3

0.13

2nd order 4th order 6th order 8th order 10th order

3.05 3

0.12

Coefficients of variation: 2,4,6th orders, γ =0.1 2nd order 4th order 6th order

0.11 E

Variances

2.95 2.9

0.1

2.85

0.09

2.8

0.08

2.75 0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

1.2

0.8

Fig. 23. Expected values; 2, 4, 6, 8, 10th orders; g:0.25.

3.5

0.5

2nd order 4th order 6th order

0.45

0.95

1 ε

1.05

1.1

1.15

1.2

Coefficients of variation: 2,4,6th orders, γ=0.25 2nd order 4th order 6th order

0.4

3 Variances

Standard deviations

0.9

Fig. 26. Coefficients of variation; 2, 4, 6th orders; g:0.10.

Standard deviations: 2,4,6th orders, γ=0.1

x 10−4

0.85

0.35 0.3

2.5 0.25 0.2

2 0.8

0.85

0.9

0.95

1 ε

1.05

1.1

Fig. 24. Standard deviations; 2, 4, 6th orders; g:0.10.

1.15

1.2

0.8

0.85

0.9

0.95

1 ε

1.05

1.1

1.15

Fig. 27. Coefficients of variation; 2, 4, 6th orders; g:0.25.

1.2

G.R. Liu et al. / Finite Elements in Analysis and Design 63 (2013) 51–61

3.15

x 10−3

approximation of the strain field. Therefore, as long as the stochastic modeling is the same, the benefit of CS-FEM will be delivered. It is expected that there is no technical difficulty in further extending the application of GS_CS-FEM into other input random variables or application to other probability distributions.

Expected values: 2,4,6,8,10th orders 2nd order 4th order 6th order 8th order 10th order MC simulation

3.1 3.05 3

61

E

2.95 2.9

References

2.85 2.8 2.75 2.7

0

0.05

0.1

0.15 γ

0.2

0.25

0.3

Fig. 28. Expected values; GS_CS-FEM vs. MCs.

Fig. 29. Standard deviations; GS_CS-FEM vs. MCs.

applied to both plane stress and plane strain examples. The efficiency and accuracy of the proposed GS_CS-FEM are validated by comparison with results of MCs. The effects of input coefficients of variation and perturbation parameters are investigated, and some conclusions can be drawn as follows:

 The 2nd order approximation may be sufficient for a small



input coefficient of variation g such as less than 0.1, but higher g needs higher order perturbation approximation. For instance, at least 10th order perturbation are required for g ¼0.3. In deterministic computation, it has already proved that the cell-based smoothed finite element method (CS-FEM) can significantly improve accuracy and convergence compared to the FEM [3,25]. The same level of improvement can be achieved for the probabilistic solutions, as expected. This is because the improvement by CS-FEM is achieved by spatial

[1] G.R. Liu, S.S. Quek, The Finite Element Method: A Practical Course, Butterworth-Heinemann, Oxford, 2003. [2] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin meshfree method, Int. J. Numer. Methods Eng. 50 (2001) 435–466. [3] G.R. Liu, T.T. Nguyen, K.Y. Dai, K.Y. Lam, Theoretical aspects of the smoothed finite element method (SFEM), Int. J. Numer. Methods Eng. 71 (2007) 902–930. [4] G.R Liu, A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods, Int. J. Comput. Methods 5 (2) (2008) 199–236. [5] G.R. Liu, A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: Part I theory and Part II applications to solid mechanics problems, Int. J. Numer. Methods Eng. 81 (2010) 1093–1126. [6] G.R. Liu, G.R. Zhang, A normed G space and weakened weak (W2) formulation of a cell-based smoothed point interpolation method, Int. J. Comput. Methods 6 (1) (2009) 147–179. [7] G.R. Liu, G.Y. Zhang, K.Y. Dai, Y.Y. Wang, Z.H. Zhong, G.Y. Li, X. Han, A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems, Int. J. Comput. Methods 2 (4) (2005) 645–665. [8] G.R. Liu, G.R. Zhang, Edge-based smoothed point interpolation methods, Int. J. Comput. Methods 5 (4) (2008) 621–646. [9] G.R. Liu, G.R. Zhang, A normed G space and weakened weak (W2) formulation of a cell-based smoothed point interpolation method, Int. J. Comput. Methods 6 (1) (2009) 147–179. [10] Z.Q. Zhang, G.R. Liu, Upper and lower bounds for natural frequencies: a property of the smoothed finite element methods, Int. J. Numer. Methods Eng. 84 (2) (2010) 149–178. [11] G.R. Liu, T.T. Nguyen, Smoothed Finite Element Methods, CRC Press, Boca Raton, FL, 2010. [12] A. Haldar, S. Mahadevan, Reliability Assessment Using Stochastic Finite Element Analysis, John Wiley & Sons, Inc, New York, 2000. [13] G. Stefanou, The stochastic finite element method: past, present and future, Comput. Meth. Appl. Mech. Eng. 198 (2009) 1031–1051. [14] M. Kamin´ski, On generalized stochastic perturbation-based finite element method, Int. J. Numer. Methods Eng. Biomed. Eng. 22 (2006) 23–31. [15] M. Kamin´ski, Generalized stochastic perturbation technique in engineering computations, Math. Comput. Modell. 51 (2010) 272–285. [16] M. Kleiber, T.D. Hien, The Stochastic Finite Element Method, John Wiley & Sons, Ltd, West Sussex, England, 1992. ¨. C [17] O - avdar, A. Bayraktar, A. C - avdar, Effects of random material and geometrical properties on structural safety of steel-concrete composite systems, Int. J. Numer. Methods Eng. Biomed. Eng. 27 (2011) 1473–1492. [18] A. Haldar, S. Mahadevan, Reliability Assessment Using Stochastic Finite Element Analysis, John Wiley & Sons, Inc, New York, 2000. [19] G. Stefanou, The stochastic finite element method: past, present and future, Comput. Meth. Appl. Mech. Eng. 198 (2009) 1031–1051. [20] K.J. Bathe, Finite Element Procedures, MIT Press, Prentice-Hall, Cambridge, MA, Englewood Cliffs, NJ, 1996. [21] K.Y. Dai, G.R. Liu, T.T. Nguyen, An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics, Finite Elem. Anal. Des. 43 (2007) 847–860. [22] G.R. Liu, T. Nguyen-Thoi, K.Y Lam, An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids, J. Sound Vib. 320 (2009) 1100–1130. [23] M. Kamin´ski, Generalized perturbation-based stochastic finite element method in elastostatics, Comput. Struct. 85 (2007) 586–594. [24] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, New York, 1970. [25] H. Nguyen-Xuan, S. Bordas, H. Nguyen-Dang, Smooth finite element methods: convergence, accuracy and properties, Int. J. Numer. Methods Eng. 74 (2008) 175–208.