A stochastic B-spline wavelet on the interval finite element method for problems in elasto-statics

Probabilistic Engineering Mechanics 58 (2019) 102996

Contents lists available at ScienceDirect

Probabilistic Engineering Mechanics journal homepage: www.elsevier.com/locate/probengmech

A stochastic B-spline wavelet on the interval finite element method for problems in elasto-statics Shashank Vadlamani, Arun C.O. ∗ Department of Aerospace Engineering, Indian Institute of Space Science and Technology, Thiruvananthapuram Kerala, 695547, India

ARTICLE

INFO

Keywords: B-spline wavelet on the interval Scaling function Multiresolution analysis Random field Auto-covariance function Perturbation method Monte Carlo simulation

ABSTRACT The current paper presents the formulation of stochastic B-spline wavelet on the interval (BSWI) based wavelet finite element method (WFEM) for one dimensional (1D) and two dimensional (2D) problems in elasto-statics wherein, the spatial variation of modulus of elasticity is modelled as a homogeneous random field. Random field discretization is done using BSWI scaling functions. Perturbation approach is used to calculate the response statistics. The results are compared with that obtained from Monte Carlo simulation (MCS) and a parametric study is carried out by considering different coefficient of variation values by varying the standard deviation. For 1D problem in particular, results from proposed stochastic WFEM method are compared with those found using stochastic finite element method (FEM) wherein random field discretization is done using Lagrange shape functions. A comparative study of the computational time needed for the execution of perturbation approach and MCS based on WFEM and FEM is also performed.

1. Introduction During the design of engineering structural systems, a deterministic approach with a factor of safety does not provide adequate information to maximize the safety of the system. On the other hand, a stochastic approach takes into account the inherent uncertainties during the design and analysis stage of the system. It provides additional statistical information about the system which leads to a more comprehensive description and optimum design of the system, thereby improving its reliability [1,2]. On the downside, a stochastic approach increases the model complexity and requires a higher computational effort to obtain the solution of the system response as compared to a deterministic approach. However, the advent of powerful computational resources during the past few decades has lead to an extensive research in the field of stochastic mechanics [3]. Nonetheless, due to the growing complexity in the system design, there is a need for the development of more efficient stochastic based numerical methods. Finite element method (FEM), being a popular numerical tool, becomes the favourite choice for solving stochastic partial differential equations. Extensive research has been done in stochastic finite element methods (SFEM), which is used to discretize the input random fields and do a response analysis [4–13]. However, the high mesh dependency of FEM creates difficulties in mapping the random field discretization to response discretization [14,15]. Further, to alleviate the mesh dependency of FEM, meshfree methods have been introduced into stochastic analysis [16–19]. Another attractive alternate tool, which can reduce ∗

the mesh dependency and re-meshing issues considerably, is wavelet finite element method (WFEM). Over the past few decades, there has been a widespread research into the development of wavelet based numerical methods in many areas of scientific importance [20]. Wavelets are mathematical functions that can be used to approximate other functions at different levels of resolution. The multiresolution analysis (MRA) [21–23] property that wavelets possess leads to the development of a hierarchy of solutions during the approximation process. Wavelets have a compact support and are localized in space, which leads to a refinement of solution locally in the regions of high gradient [24]. Therefore, issues such as slow convergence in the vicinity of high gradients and re-meshing that are encountered in conventional FEM can be avoided using wavelet based numerical methods. As a result, the computational effort that is needed at pre-processing stage in adaptive mesh sensitive problems can be reduced by making use of the properties of wavelets. There exist different wavelet based numerical methods in the literature [20]. Among these, B-spline wavelet on the interval (BSWI) based WFEM has gained wide prominence due to the underlying properties of BSWI [22,25] that make it more efficient as compared to other wavelet based numerical methods. In view of the same reasons that make BSWI more efficient, it is selected for usage in the current paper. One of the first papers to implement spline wavelets expansions as shape functions in the framework of finite element analysis was by Chen and Wu [26,27]. Xiang et al. [28] proposed the construction of

Corresponding author. E-mail addresses: [email protected] (S. Vadlamani), [email protected] (Arun C.O.).

https://doi.org/10.1016/j.probengmech.2019.102996 Received 30 April 2019; Accepted 13 August 2019 Available online 19 August 2019 0266-8920/© 2019 Published by Elsevier Ltd.

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

of the random field is established with the help of a convergence study. Correlation length parameters are varied and its effects on the response statistics are also analysed in the 1D case. Since, BSWI scaling functions form the basis for the current study, in the next section, for the benefit of the reader, a brief description of BSWI is given.

one-dimensional (1D) 𝐶0 type and 𝐶1 type BSWI elements for structural analysis using BSWI WFEM. A class of 𝐶0 type plate elements for plane elasto-mechanics and moderately thick plate problems were constructed by Xiang et al. [29] based on two-dimensional (2D) tensor product of BSWI scaling functions. Further, Xiang et al. [30] using BSWI WFEM, proposed a new approach for the analysis of stress intensity factors (SIFs) for cracked plane plate. More recently, BSWI plane elastomechanical plate element was used to develop a numerical method to calculate the elastic band structures of 2D phononic crystals (2DPCs), composed of square lattices of solid cuboids in a solid matrix [31]. The usage of wavelets in the framework of stochastic processes can be traced back to Wornell [32], who used orthonormal wavelet bases to construct for scaling processes from a set of uncorrelated random variables. Zeldin and Spanos [33] have shown that Daubechies wavelet bases can be used effectively to model random fields. However, Daubechies wavelets lack explicit function expression and a separate algorithm is needed for the evaluation of integrals of products of wavelet basis functions, with or without derivatives, known as connection coefficients [34,35]; which increases computational effort. Further, it is noted that connection coefficients derivation can only be obtained for integration in global coordinates and it fails when the integrand involves variant Jacobian [26,27,29]. Phoon et al. [36] and Proppe [37,38] implemented a wavelet Galerkin approach using basis functions of Haar wavelets to solve the Fredholm integral equation involved in Karhunen–Loève (KL) expansion technique. Further, Phoon et al. [39] showed that use of a wavelet-Galerkin scheme for determination of Eigen-solutions in KL expansion is computationally equivalent to using wavelets directly for stochastic expansion of a Gaussian process. However, the drawback of using KL expansion is that the values of true variance are underpredicted [40]. Angulo and Ruiz-Medina [41] addressed an inverse problem wherein an input random field is derived from the observation of an output random field using the orthogonal expansion of secondorder random fields in terms of a wavelet basis. Le Maitre et al. [42] constructed an uncertainty quantification scheme based on Polynomial Chaos (PC) representations, in which a Haar or a Legendre basis was used for the orthogonal projection of uncertain data and solution variables. Later, Spanos et al. [43] estimated the power spectrum of non-stationary stochastic processes using dyadic, generalized, and filtered harmonic wavelets. Spline wavelet-based FEM along with Monte Carlo simulation (MCS) was used for bending analysis of plates by Han et al. [44] when parameters are random variables. Andreev [45], investigated the numerical aspects of the stochastic Galerkin method when applied to linear stochastic PDEs. The spatial variable was discretized using the piecewise linear and quadratic spline wavelet-Galerkin bases which lead to well-conditioned system matrices. Zhang et al. [46] used BSWI based finite element method by combining with MCS for structural analysis. It is noticed that the existing literature does not discuss about modelling of random field using BSWI WFEM. Even though BSWI WFEM has been used with MCS for structural analysis as mentioned in the preceding paragraph, MCS makes the whole process computationally expensive and is not a viable option. Hence, in the present study, a formulation of stochastic BSWI WFEM for 1D and 2D problems in elasto-statics using perturbation method is shown. The spatial variation of modulus of elasticity is modelled as a homogeneous random field. For the discretization of both random field and response, it is proposed to use BSWI scaling functions, whose MRA properties can be used to refine the discretization without much computational effort. Numerical examples based on 1D and 2D formulations are solved and compared with solutions from MCS. A parametric study is also done to estimate the effect of coefficient of variation (CV) of the input field on results of perturbation. For 1D problem, a comparative study is done with stochastic FEM wherein the random field is modelled using Lagrange shape functions. The size of random variable vector for discretization

2. B-spline wavelet on the interval [0, 1] The theory of spline wavelets for whole square integrable real space 𝐿2 (𝐑) was developed by Chui and Wang [47–49]. Wavelets defined on 𝐿2 (𝐑) cannot be directly used as interpolating functions, as it would result in numerical instability [50]. This issue was addressed by Chui and Quak [25] who constructed wavelet bases for the bounded interval [0, 1], known as BSWI. The decomposition and reconstruction algorithms for the spline wavelets on the bounded interval were developed by Quak and Weyrich [51]. Further, BSWI scaling and wavelet functions were also used by Goswami et al. [52] for solving first kind integral equations. Spline wavelets are semi-orthogonal. These wavelets retain interscale orthogonality and there is no necessity for the basis functions to be orthogonal to its translates within the same resolution level. Splines can readily adapt to the case of the bounded interval [0, 1] by introducing multiple knots at the endpoints. As a result, no truncation is needed when the function on 𝐿2 (𝐑) is restricted and by way of suitable adaptation at the endpoints, MRA of 𝐿2 (𝐑) can be implemented over to [0, 1]. However, at points where multiple knots exist (at 0 and 1 in the case of BSWI), the smoothness decreases due to knot coalescence, but remains unaffected elsewhere. The continuity of B-splines depends on the selected order 𝑚 in such a way that B-splines with order 𝑚 are in 𝐶𝑚−2 continuity [28]. The analytical expressions for the BSWI scaling functions 𝜙 and wavelet functions 𝜓 for given order 𝑚 and resolution 𝑗 = 0 can be found in the paper by Goswami et al. [52] and the expressions for order 𝑚 and any resolution 𝑗 were given in the paper by Xiang et al. [28] as, 𝜙𝑙𝑚,𝑘 (2𝑗−𝑙 𝜉), 𝑘 = −𝑚 + 1, … , −1, ⎫ ⎪ (0 boundary scaling functions) ⎪ 𝑙 𝑗−𝑙 𝑗 𝑗 𝜙 (1 − 2 𝜉), 𝑘 = 2 − 𝑚 + 1, … , 2 − 1,⎪ 𝜙𝑗𝑚,𝑘 (𝜉) = 𝑚,2𝑗 −𝑚−𝑘 ⎬ (1 boundary scaling functions) ⎪ ⎪ 𝜙𝑙𝑚,0 (2𝑗−𝑙 𝜉 − 2−𝑙 𝑘), 𝑘 = 0, … , 2𝑗 − 𝑚, ⎪ (inner scaling functions) ⎭

(1)

⎫ ⎪ ⎪ 𝑗 𝜓𝑚,𝑘 (𝜉) = ⎬ (1 boundary wavelets) ⎪ 𝑙 𝑗−𝑙 −𝑙 𝑗 𝜓𝑚,0 (2 𝜉 − 2 𝑘), 𝑘 = 0, … , 2 − 2𝑚 + 1, (inner wavelets)⎪ ⎭ 𝑙 (2𝑗−𝑙 𝜉), 𝑘 = −𝑚 + 1, … , −1, (0 boundary wavelets) 𝜓𝑚,𝑘 𝑙 𝑗−𝑙 𝜉), 𝑘 = 2𝑗 − 2𝑚 + 2, … , 2𝑗 − 𝑚 𝜓𝑚,2 𝑗 −2𝑚−𝑘+1 (1 − 2

(2) The compactly supported intervals of wavelets are, [0, (2𝑚 − 1 + 𝑘)2−𝑗 ], (0 boundary wavelets)⎫ ⎪ 𝑗 supp 𝜓𝑚,𝑘 (𝜉) = [𝑘2−𝑗 , 1], (1 boundary wavelets) ⎬ −𝑗 −𝑗 [𝑘2 , (2𝑚 − 1 + 𝑘)2 ], (inner wavelets) ⎪ ⎭

(3)

BSWI scaling functions are categorized as the boundary scaling functions that exist at boundary points 0 and 1 on the domain and inner scaling functions that are dilations and translations of cardinal B-splines as shown in Eq. (1), (2) and (3). The corresponding wavelets can be constructed by utilizing the scaling functions eventually. For illustration, the scaling functions of BSWI using different order at resolution 𝑗 = 3 are shown in Fig. 1. In the next section, the formulation of 𝐶0 continuity BSWI WFEM elements for 1D and 2D cases is shown. 2

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 1. BSWI scaling functions using different order at resolution 𝑗 = 3.

and resolution of 𝑗, the 1D element solving domain 𝛺𝑒 is discretized using 𝑛 + 1 nodes where, 𝑛 = 2𝑗 + 𝑚 − 2 segments. The generalized functional of potential energy for a bar element spanning between the global co-ordinates 𝑥𝑎 and 𝑥𝑏 can be written as, { } 𝑛𝑑 𝑥𝑏 𝑥𝑏 ∏ ∑ ( ) 𝑑𝑢0 2 1 = 𝐸𝐴 𝑑𝑥 − 𝑓 (𝑥)𝑢0 𝑑𝑥 − 𝑄 𝑖 𝑢0 𝑥 𝑖 (7) ∫𝑥𝑎 2 ∫𝑥𝑎 𝑑𝑥

Fig. 2. Distribution of nodes in a single 1D element with 𝐶0 continuity.

𝑖=1

(4)

Here, 𝑓 (𝑥) is the distributed load, 𝑄𝑖 are the generalized lumped nodal forces, 𝑛𝑑 is the total number of nodes per BSWI wavelet finite bar element, 𝑢0 is the axial displacement and 𝐸 is the modulus of elasticity. The formulation of a BSWI WFEM bar element was given by Xiang et al. [28]. The unknown axial displacement field of a bar element in the element solving domain 𝜉 is approximated in terms of wavelet scaling functions as,

(5)

𝑢0 (𝜉) =

3. Formulation of BSWI element for elasto-statics The equilibrium equations and boundary conditions for small displacements in linear elastic solids are, 𝛁 ⋅ 𝝈 + 𝒃 = 0 in 𝛺 𝝈 ⋅ 𝒏 = 𝒕 on 𝛤𝑡 (natural boundary conditions)

𝑗 −1 2∑

𝑎𝑗𝑚,𝑘 𝜙𝑗𝑚,𝑘 (𝜉) = 𝝋𝒂𝒆

(8)

𝑘=−𝑚+1

𝒖 = 𝒖 on 𝛤𝑢 (essential boundary conditions)

where, 𝝋 = {𝜙𝑗𝑚,−𝑚+1 (𝜉) … 𝜙𝑗𝑚,2𝑗 −1 (𝜉)} is the row vector of BSWI scaling

(6)

functions and 𝒂𝒆 = {𝑎𝑗𝑚,−𝑚+1 𝑎𝑗𝑚,−𝑚+2 … 𝑎𝑗𝑚,2𝑗 −1 }𝑇 is the column vector of wavelet coefficients that needs to be determined. The 𝐶0 element type transformation matrix and physical DOF are obtained as [28],

Here 𝝈 = 𝑫𝜺 is the stress vector, D is the material property matrix, 𝜺 = 𝛁𝑠 𝒖 is the strain vector, 𝒖 is the displacement vector, 𝒃 is the body force vector, 𝒕 and 𝒖 are the vectors of prescribed surface tractions and displacements respectively, 𝒏 is a unit normal to the domain 𝛺. 𝛤𝑡 and 𝛤𝑢 are the portions of boundary { 𝛤 , where tractions and displacements } are prescribed respectively. 𝛁 = 𝜕𝑥𝜕 , 𝜕𝑥𝜕 … is the vector of gradient 1 2 operators and 𝛁𝑠 𝒖 is the symmetric part of 𝛁𝒖.

𝒖𝒆 = 𝑹𝒆 𝒂𝒆 , ⎫ 𝑹𝒆 = [𝝋𝑇 (𝜉1 )𝝋𝑇 (𝜉2 ) … 𝝋𝑇 (𝜉𝑛+1 )]𝑇 ,⎪ ⎪ { } 𝑇 𝒖𝒆 = 𝑢1 𝑢2 … 𝑢𝑛+1 , ⎬ ⎪ 𝑻 𝑒 = (𝑹𝒆 )−1 , ⎪ 𝑵 𝑒 = 𝝋𝑻 𝑒 ⎭

3.1. 1D BSWI WFEM bar element

(9)

The elemental transformation matrix transforms the stiffness matrix from wavelet space into physical space. The transformation matrix maintains the continuity and compatibility within the element and at the interface between the neighbouring elements; it is maintained by using an assembly matrix. By substituting Eq. (9) into Eq. (8), the unknown displacement field is obtained as, ( )−1 𝒆 𝑢0 (𝜉) = 𝝋 𝑹𝒆 𝒖 = 𝝋𝑻 𝑒 𝒖𝒆 = 𝑵 𝑒 𝒖𝒆 (10)

In BSWI WFEM, the 1D problem domain 𝛺 is divided into subdomains 𝛺𝑖 (𝑖 = 1, 2. . . ). Each 𝛺𝑖 is then mapped into the standard element solving domain 𝛺𝑒 = {𝜉∣ 𝜉 ∈ [0, 1]} as shown in Fig. 2, where instead of using the traditional polynomial interpolation, scaling or scaling and wavelet functions of BSWI are used to form the shape functions over the elements 𝛺𝑒 . Here 𝜉 is the local co-ordinate used for solving 1D BSWI on [0, 1] along the axis of the element. For a BSWI scaling function of order 𝑚

Upon substituting the displacement field of Eq. (10) into the weak form and invoking the stationary condition for variation of admissible 3

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

where, 𝝋 is obtained by taking 2D tensor product of BSWI scaling functions in the element solving domain 𝛺𝑒 as, ⎫ 𝝋 = 𝝋1 ⊗ 𝝋2 ⎪ 𝝋1 = {𝜙𝑗𝑚,−𝑚+1 (𝜉)𝜙𝑗𝑚,−𝑚+2 (𝜉) … 𝜙𝑗𝑚,2𝑗 −1 (𝜉)} ⎬ 𝑗 𝑗 𝑗 𝝋2 = {𝜙𝑚,−𝑚+1 (𝜂)𝜙𝑚,−𝑚+2 (𝜂) … 𝜙𝑚,2𝑗 −1 (𝜂)}⎪ ⎭ 𝒄 𝑒 = {𝑐1 , 𝑐2 … 𝑐(𝑛+1)2 }𝑇 𝒅 𝑒 = {𝑑1 , 𝑑2 … 𝑑(𝑛+1)2 }𝑇

(15)

} (16)

Here, 𝝋1 , 𝝋2 are the rows vectors combined by the scaling functions for order 𝑚 and resolution 𝑗 and 𝒄 𝑒 , 𝒅 𝑒 are the column vector of wavelet coefficients to be determined. The column vector of physical degrees of freedom (DOF) at nodes is defined as, 𝒖𝑒 = {𝑢(𝜉1 , 𝜂1 ) … 𝑢(𝜉1 , 𝜂𝑛+1 ) … 𝑢(𝜉𝑛+1 , 𝜂𝑛+1 )}𝑇 𝒗𝑒 = {𝑣(𝜉1 , 𝜂1 ) … 𝑣(𝜉1 , 𝜂𝑛+1 ) … 𝑣(𝜉𝑛+1 , 𝜂𝑛+1 )}𝑇

} (17)

By substituting Eq. (14) into Eq. (17), the 2D elemental transformation Fig. 3. Distribution of nodes in a single 2D element with 𝐶0 continuity.

matrix 𝑹𝑒 can be obtained as, 𝒖𝑒 = 𝑹𝑒 𝒄 𝑒 , 𝒗𝑒 = 𝑹𝑒 𝒅 𝑒

displacements, the element solving equations for BSWI WFEM bar element can be obtained as,

} (18)

such that, 𝑲 𝑒 𝒖𝒆 = 𝑭 𝑒

(11)

⎫ 𝑹𝑒 = 𝑻 𝑒1 ⊗ 𝑻 𝑒2 , ⎪ 𝑻 𝑒1 = {𝝋𝑇1 (𝜉1 )𝝋𝑇1 (𝜉2 ) … … 𝝋𝑇1 (𝜉𝑛+1 )}𝑇 ,⎬ 𝑒 𝑇 𝑇 𝑇 𝑇 𝑻 2 = {𝝋2 (𝜂1 )𝝋2 (𝜂2 ) … … 𝝋2 (𝜂𝑛+1 )} ⎪ ⎭

where, (

)𝑇 (

)

⎫ 𝑑𝝋 𝑑𝝋 𝐴 ⎪ 𝐸 (𝑻 𝑒 )𝑇 𝑲𝑒 = (𝑻 𝑒 ) 𝑑𝜉, ∫0 𝑙𝑒𝑥 𝑑𝜉 𝑑𝜉 ⎪ ⎬ 𝑛𝑑 1 ∑ ( ) 𝑇 𝑇 𝑭𝑒 = 𝑙𝑒𝑥 𝑓 (𝜉) (𝑻 𝑒 ) (𝝋)𝑇 𝑑𝜉 + 𝑄𝑒𝑖 (𝑻 𝑒 ) (𝝋)𝑇 𝜉𝑖 ⎪ ⎪ ∫0 𝑖=1 ⎭ 1

(19)

Hence, the equations for displacement field in the standard element

(12)

solving domain can be replaced as, 𝑢0 (𝜉, 𝜂) = 𝝋𝑇 (𝑹𝑒 )−1 𝒖𝑒 = 𝑵 1 𝒖𝑒 , 𝑣0 (𝜉, 𝜂) = 𝝋𝑇 (𝑹𝑒 )−1 𝒗𝑒 = 𝑵 1 𝒗𝑒

3.2. BSWI WFEM plane elasto-statics element

} (20)

By using the principle of minimum potential, the element stiffness In BSWI WFEM, the problem in 2D domain 𝛺 is divided into subdomains 𝛺𝑖 (𝑖 = 1, 2. . . ) where, each 𝛺𝑖 is then mapped into the standard 2D element solving domain 𝛺𝑒 = {𝜉, 𝜂∣ 𝜉, 𝜂 ∈ [0, 1]}. Here, 𝜉, 𝜂 are the local co-ordinates used for solving 2D BSWI on [0, 1] as shown in Fig. 3. For a BSWI scaling function of order 𝑚 and resolution of 𝑗, the 2D element solving domain 𝛺𝑒 is discretized into (𝑛 + 1) × (𝑛 + 1) nodes where, 𝑛 = 2𝑗 + 𝑚 − 2 segments on each side. Instead of using the traditional polynomial interpolation, 2D tensor product of scaling or scaling and wavelet functions of BSWI are used to form the shape functions over the elements 𝛺𝑒 in 2D domain.

matrix and force vector (for a unit thickness) can be obtained as [29],

The generalized functional of potential energy for a sub-domain 𝛺𝑖 in a 2D plane stress problem can be written as,

𝑷 e𝑎

𝛱𝑝 (𝑢) =

∫𝛺𝑖

𝑲 𝒆𝑼 𝒆 = 𝑭 𝒆

(21)

where, 𝑲𝒆 =

[

] 𝑲 e,1 𝑲 e,2 , 𝑲 e,3 𝑲 e,4

𝑼𝒆 =

{ e} 𝒖 , 𝒗e

𝑭𝒆 =

{

𝑷 e𝒂 𝑷 e𝒃

}} (22)

such that, (

⎫ 𝑏𝑥 𝝋𝑇 d𝛺 ⎪ ⎪ 𝑛𝑒 ⎪ ∑ 𝑒 −1 𝑇 𝑇 ⎪ + ((𝑹 ) ) 𝝋 (𝜉𝑖 , 𝜂𝑖 )𝐹𝑖𝑥 , ⎪ 𝑖=1 ( ⎬ 𝑷 e𝑏 = ((𝑹𝑒 )−1 )𝑇 𝑝𝑦 𝝋𝑇 ds + 𝑏𝑦 𝝋𝑇 d𝛺 ⎪ ⎪ ∫𝑆𝜎 ∫ )𝛺𝑒 ⎪ 𝑛𝑒 ∑ ⎪ + ((𝑹𝑒 )−1 )𝑇 𝝋𝑇 (𝜉𝑖 , 𝜂𝑖 )𝐹𝑖𝑦 , ⎪ ⎭ 𝑖=1

𝑛𝑒 ∑ 1 𝑇 𝜺 𝑫 𝜺 𝑡𝑑𝑥𝑑𝑦 − 𝒖𝑇 𝒃 𝑡𝑑𝑥𝑑𝑦 − 𝒖𝑇 𝒑 𝑡𝑑𝛤 − 𝒖𝑇𝑖 𝑭 𝑖 (13) ∫𝛺𝑖 ∫𝛤𝑡 2 𝑖=1

such that, 𝒃 is the body force column vector, 𝒖 is the displacement vector, 𝑡 is the uniform thickness and 𝒑 is the column vector of surface tractions in x and y directions respectively. 𝑭 is column vector of point forces, 𝒖𝑖 is the column vector of point displacements at the point of application of forces and 𝑛𝑒 is the total number of nodes per 2D BSWI plane elasto-statics element. The formulation of a BSWI WFEM plane elasto-statics element was given by Xiang et al. [29].

𝑒 −1 𝑇

= ((𝑹 ) )

∫𝑆𝜎

𝑝𝑥 𝝋𝑇 ds +

∫ )𝛺𝑒

(23)

and ⎫ 1 − 𝜈 00 𝐸 (𝐀11 ⊗ 𝐀00 + 𝐀1 ⊗ 𝐀11 ), ⎪ 2 2 2 (1 − 𝜈 2 ) 1 ⎪ ⎪ 𝐸 1 − 𝜈 10 01 01 10 = (𝜈𝐀1 ⊗ 𝐀2 + 𝐀1 ⊗ 𝐀2 ), ⎬ 2 (1 − 𝜈 2 ) ⎪ 𝐸 1 − 𝜈 11 = (𝑲 e,2 )𝑇 , 𝑲 e,4 = (𝐀00 ⊗ 𝐀11 + 𝐀1 ⊗ 𝐀00 ),⎪ 2 1 2 ⎪ 2 (1 − 𝜈 2 ) ⎭

𝑲 e,1 =

The displacement field function in the local co-ordinate system can be approximated as, } 𝑢0 (𝜉, 𝜂) = 𝝋 𝒄 𝑒 (14) 𝑣0 (𝜉, 𝜂) = 𝝋 𝒅 𝑒

𝑲 e,2 𝑲 e,3

4

(24)

S. Vadlamani and Arun C.O.

with

(

Probabilistic Engineering Mechanics 58 (2019) 102996

𝐀00 = ((𝑻 𝑒1 )−1 )𝑇 1 𝐀01 = ((𝑻 𝑒1 )−1 )𝑇 1 𝐀10 = ((𝑻 𝑒1 )−1 )𝑇 1 𝐀11 = ((𝑻 𝑒1 )−1 )𝑇 1

is necessary. But 𝛼(𝐱) being a random field does not possess an explicit expression and hence requires an approximation. Generally, this is done by approximating a function over a set of random variables distributed in the domain obtained by discretization of the field and is known as random field discretization. In the current work, it is proposed to use a shape function method for random field modelling. Shape function method is found to give better approximation when compared to KL expansion technique. Shape function method using Lagrange interpolation and moving least square shape functions has been employed in SFEM [40] and stochastic meshless methods [18] respectively. In the present study, it is proposed to use wavelet functions to model both the random field and response. A detailed discussion on the modelling the random field in 1D and 2D domain and stochastic BSWI WFEM formulation is given in the following section.

)

⎫ ( )𝑇 ( ) ⎪ 𝝋1 d𝜉 (𝑻 𝑒1 )−1 , 𝑙𝑒𝑥 𝝋1 ∫0 ⎪ ( ⎪ ( ) ) 1( )𝑇 𝑑𝝋1 ⎪ 𝝋1 d𝜉 (𝑻 𝑒1 )−1 , ⎪ ∫0 𝑑𝜉 ⎪ ( ) ⎬ 1 ( 𝑑𝝋 )𝑇 ( ) ⎪ 1 𝝋1 d𝜉 (𝑻 𝑒1 )−1 , ⎪ ∫0 𝑑𝜉 ⎪ ( ) ( ) ( ) ⎪ 𝑇 1 𝑑𝝋 𝑑𝝋 1 1 1 d𝜉 (𝑻 𝑒1 )−1 ,⎪ ⎪ 𝑙𝑒𝑥 ∫0 𝑑𝜉 𝑑𝜉 ⎭ 1

(25)

Here, 𝝋1 , 𝝋2 are the vectors of associated wavelet scaling functions. Further, 𝐀𝑝𝑞 (𝑝, 𝑞 = 0, 1) is analogous to 𝐀𝑝𝑞 (𝑝, 𝑞 = 0, 1) when 𝑙𝑒𝑥 , 𝑑𝜉, 𝑻 𝑒1 2 1 𝑒 and 𝝋1 are replaced with 𝑙𝑒𝑦 , 𝑑𝜂, 𝑻 2 and 𝝋2 respectively. The essential boundary conditions are implemented by setting the corresponding DOF to the respective value and eliminating them from the system of equations. The element stiffness matrices formulated as given in Eq. (21) for all the sub-domains can be assembled to form the algebraic set of global equilibrium equations as,

4.1. 1D stochastic BSWI WFEM bar element On similar lines, as the displacement field is approximated in Eq. (8), the unknown Gaussian random field can be approximated in the element solving domain in terms of wavelet scaling functions as,

(26)

𝑲𝑼 = 𝑭 4. Formulation of stochastic BSWI element for elasto-statics

𝛼(𝜉) =

𝑗

𝑟

(33)

𝑟

𝑗

where, 𝝋R = {𝜙𝑚𝑟

𝑗 (𝜉)} 𝑚𝑟 ,2𝑗𝑟 −1 𝑗𝑟 𝑗𝑟 {𝑏𝑚 ,−𝑚 +1 𝑏𝑚 ,−𝑚 +2 𝑟 𝑟 𝑟 𝑟

𝑟 ,−𝑚𝑟 +1

𝑒

(𝜉) … 𝜙 𝑟

is the row vector of BSWI 𝑗

scaling functions, 𝒃 = … 𝑏 𝑟 𝑗𝑟 }𝑇 is the column 𝑚𝑟 ,2 −1 vector of wavelet coefficients that needs to be determined and 𝑚𝑟 , 𝑗𝑟 are the order and resolution of BSWI scaling function chosen for the discretization of random field. The subscript 𝑟 used here denotes that the function or variable used is associated with random field. It can be noted from Eqs. (8) and (33) that the order and resolution that is used for the discretization of the displacement field and random field can be different from each other. The unknown random field function can be expressed in terms of 𝐶0 element type transformation matrix as, ( )−1 𝑒 𝛼(𝜉) = 𝝋R 𝑹𝑒R 𝜶 R = 𝝋R 𝑻 𝑒R 𝜶 𝑒R = 𝑵 𝑒R 𝜶 𝑒R (34) { }𝑇 𝑒 where, 𝜶 R = 𝛼1R 𝛼2R … 𝛼(𝑛+1)R is the set random variables distributed over the domain of the element. Thus, element stiffness coefficients and hence the element displacements will become functions of random variables 𝜶 𝑒R and in general Eq. (11) can be written as,

(27)

where, 𝛼(𝐱) is the zero mean Gaussian field with exponential autocovariance kernel given as, [ ( | | )] 𝛥𝑖 𝛤𝛼 = 𝜎𝛼2 . exp − 𝛴𝑖𝑛 | | . (28) 𝑐𝑖

( ) ( ) 𝑲 𝑒 𝜶 𝑒R 𝒖𝑒 𝜶 𝑒R = 𝑭 𝑒

where, 𝑛 ∈ R𝑛 , 𝛥𝑖 is the distance between two points 𝑥𝑎 , 𝑥𝑏 along 𝑖, 𝑐𝑖 is the correlation length parameter which determines the statistical correlation of field variable in the domain and 𝜎𝛼 is the standard deviation of the zero mean random field which can be expressed in terms of standard deviation 𝜎𝐸 of 𝐸(𝐱) as,

(35)

When 𝐸(𝐱) is modelled as a homogeneous Gaussian field as given in Eq. (27), the 𝑲 𝑒 in Eq. (35) can be written as, 𝑲𝑒 =

(29)

𝜎𝛼 = 𝜎𝐸 ∕𝜇𝐸 .

𝑗

𝑏𝑚𝑟 ,𝑘 𝜙𝑚𝑟 ,𝑘 (𝜉) = 𝝋R 𝒃𝑒

𝑘=−𝑚𝑟 +1

In the current work, the Young’s modulus is considered as a spatially varying homogeneous random field. Consequently, the equilibrium equations given in Eq. (4), and the generalized functional of total potential for 1D and 2D problems as given in Eqs. (7) and (13) will also become stochastic in nature and hence the response too. This necessitates for modelling of random fields present in the physical system. In order to simplify the algebraic equations required to estimate the statistics, it is a general practise to transform the field to a Gaussian field before modelling [18]. In the current study two types random fields are considered to model 𝐸(𝐱), homogeneous Gaussian field and homogeneous lognormal field. If 𝐸(𝐱) is Gaussian random field with mean 𝜇𝐸 , it can be written as [18], 𝐸(𝐱) = 𝜇𝐸 (1 + 𝛼(𝐱)),

𝑗𝑟 −1 2∑

𝜇𝐸 𝐴 1 𝑒 𝑇 (𝑻 ) 𝑙𝑒𝑥 ∫0

(

𝑑𝝋 𝑑𝜉 (

) 𝜇𝐸 𝐴 1 ( 𝝋R 𝑻 𝑒R 𝜶 𝑒R (𝑻 𝑒 )𝑇 𝑙𝑒𝑥 ∫0

)𝑇 (

𝑑𝝋 𝑑𝜉

𝑑𝝋 𝑑𝜉 )𝑇 (

)

⎫ ⎪ ⎪ ⎬ ) 𝑑𝝋 ⎪ (𝑻 𝑒 ) 𝑑𝜉,⎪ 𝑑𝜉 ⎭ (𝑻 𝑒 ) 𝑑𝜉

(36)

However, when 𝐸(𝐱) is a homogeneous lognormal field with mean𝜇𝐸𝑙 and standard deviation 𝜎𝐸𝑙 it can be expressed in terms of 𝛼(𝐱) as,

+

𝐸(𝐱) = 𝐶𝑙 𝑒𝛼(𝐱)

Here, 𝑲 𝑒 is the elemental stochastic stiffness matrix. When 𝐸(𝐱) is modelled as a homogeneous lognormal field as given in Eq. (30), the 𝑲 𝑒 in Eq. (35) can be written as,

(30)

with 2 𝜇𝐸 𝑙 𝐶𝑙 = √ 2 + 𝜎2 𝜇𝐸 𝐸 𝑙

(31)

𝑙

𝑲𝑒 =

The auto-covariance kernel for 𝛼(𝐱) can be written as [53], 𝛤𝛼𝑙

[ ( | | )] ⎛ 𝜎𝐸2 ⎞ 𝛥𝑖 = ln ⎜1 + 𝑙 ⎟ exp − 𝛴𝑖𝑛 | | 2 ⎟ ⎜ 𝑐𝑖 𝜇𝐸 ⎝ 𝑙 ⎠

√ 𝑙𝑒𝑥

( 1+

𝜎𝐸

)2 ∫0 𝑙

𝜇𝐸

𝑙

(32)

(

1

𝜇𝐸𝑙 𝐴

exp

𝝋R 𝑻 𝑒R 𝜶 𝑒R

)

(𝑻 𝑒 )𝑇

(

𝑑𝝋 𝑑𝜉

)𝑇 (

𝑑𝝋 𝑑𝜉

)

⎫ ⎪ ⎪ (𝑻 𝑒 ) 𝑑𝜉 ⎬ ⎪ ⎪ ⎭ (37)

The formulation of 2D 𝐶0 continuity stochastic BSWI WFEM plane stress element is shown in the following section.

To evaluate the integrals defined in Eqs. (12) and (25), an explicit expression for 𝐸(𝐱), and hence 𝛼(𝐱) as a function of spatial coordinates 5

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

{ }𝑁 small [55]. Let 𝜦 = 𝛼𝑖 𝑖=1 denote the vector of 𝑁 zero mean random variables representing the random field in the global domain 𝛺. The Taylor series expansion of the global stochastic stiffness matrix 𝑲 and the displacement vector 𝑼 can be obtained as, } 𝑁 𝑁 𝑁 ∑ 1 ∑ ∑ 𝐼𝐼 𝑲 = 𝑲0 + 𝑲 𝑖𝑗 𝜶 𝑖 𝜶 𝑗 + ⋯ , (47) 𝑲 𝐼𝑖 𝜶 𝑖 + 2 𝑖=1 𝑗=1 𝑖=1 } 𝑁 𝑁 𝑁 ∑ 1 ∑ ∑ 𝐼𝐼 𝐼 𝑼 = 𝑼0 + 𝑼 𝑖 𝜶𝑖 + 𝑼 𝜶 𝜶 + ⋯, (48) 2 𝑖=1 𝑗=1 𝑖𝑗 𝑖 𝑗 𝑖=1

4.2. Stochastic BSWI WFEM plane elasto-statics element The unknown Gaussian random field can be approximated by taking 2D tensor product of BSWI scaling functions on the same lines as the displacement field is approximated in Eq. (14) as, 𝛼(𝜉, 𝜂) = 𝝋R 𝒈𝑒

(38)

where, ⎫

𝝋R = 𝝋1R ⊗ 𝝋2R 𝑗 𝑗 𝝋1R = {𝜙𝑚𝑟 ,−𝑚 +1 (𝜉)𝜙𝑚𝑟 𝑗

𝑟

𝝋2R = {𝜙𝑚𝑟

𝑟

𝑟 ,−𝑚𝑟

𝑗

𝑟 ,−𝑚𝑟 +2

(𝜂)𝜙𝑚𝑟 +1

𝑟 ,−𝑚𝑟

𝒈𝑒 = {𝑔1 , 𝑔2 … 𝑔(𝑛+1)2 }𝑇

(𝜉) …

(𝜂) … +2

⎪ 𝑗 𝜙 𝑟 𝑗𝑟 (𝜉)} ⎪ 𝑚𝑟 ,2 −1 𝑗 ⎬ 𝜙 𝑟 𝑗𝑟 (𝜂)}⎪ 𝑚𝑟 ,2 −1

Upon substituting Eqs. (47) and (48) in the equilibrium equations and equating the terms with the same order we obtain,

(39)

⎪ ⎭

⎫ 𝜕 2 𝑲 || , ⎪ | 𝜕𝜶 𝑖 𝜕𝜶 𝑗 ||𝜦=0 ⎪ ( 𝐼 ) ⎬ 𝑼 0 = 𝑲 −1 𝑭 0 , 𝑼 𝐼𝑖 = −𝑲 −1 𝑲𝑖 𝑼 0 , 0 0 ( )⎪ ⎪ −1 𝑼 𝐼𝐼 𝑲 𝐼𝑖 𝑼 𝐼𝑗 + 𝑲 𝐼𝑗 𝑼 𝐼𝑖 + 𝑲 𝐼𝐼 𝑖𝑗 = −𝑲 0 𝑖𝑗 𝑼 0 ⎭ 𝑲 𝐼𝑖 =

The unknown random field function can be expressed in terms of 𝐶0 element type 2D transformation matrix as, } 𝛼(𝜉, 𝜂) = 𝝋𝑇R (𝑹𝑒R )−1 𝜶 𝑒R = 𝑵 2 𝜶 𝑒R (40)

First order approximation 𝝁𝑈 = 𝑼 0 , 𝜸𝑈 =

𝑲 𝑒,2 𝛽 𝑲 𝑒,3 𝛽 𝑲 𝑒,4 𝛽

(

𝑼 𝐼𝑖 𝑼 𝐼𝑗

)𝑇

⎫ ⎪ ⎬ 𝜞 𝜶 (𝜶 𝑖 , 𝜶 𝑗 )⎪ ⎭

(50)

Second order approximation 𝑁 𝑁 1 ∑ ∑ 𝐼𝐼 𝑼 𝜞 (𝜶 , 𝜶 ), 2 𝑖=1 𝑗=1 𝑖𝑗 𝜶 𝑖 𝑗

⎫ ⎪ ⎪ ⎪ 𝑁 ∑ 𝑁 ( )𝑇 ∑ ⎪ 𝑼 𝐼𝑖 𝑼 𝐼𝑗 𝜞 𝜶 (𝜶 𝑖 , 𝜶 𝑗 )+ 𝜸𝑈 = ⎪ ⎬ 𝑖=1 𝑗=1 ⎪ 𝑁 𝑁 𝑁 𝑁 ⎪ 1 ∑ ∑ ∑ ∑ 𝐼𝐼 ( 𝐼𝐼 )𝑇 { 𝑼 𝑖𝑗 𝑼 𝑘𝑙 𝜞 𝜶 (𝜶 𝑖 , 𝜶 𝑙 )𝜞 𝜶 (𝜶 𝑗 , 𝜶 𝑘 ) ⎪ 4 𝑖=1 𝑗=1 𝑘=1 𝑙=1 ⎪ } ⎪ ⎭ + 𝜞 𝜶 (𝜶 𝑖 , 𝜶 𝑘 )𝜞 𝜶 (𝜶 𝑗 , 𝜶 𝑙 ) 𝝁𝑈 = 𝑼 0 +

(43)

(51)

Similarly statistics of other response functions of interest, like stresses and strains can also be found out. Further, the study can also be extended to a reliability analysis on the same lines as shown in Rahman and Rao [56]. In the next section, a few 1D and 2D numerical examples are solved based on the preceding formulations and the results are analysed.

such that, ) ( ) ( ) ( ) 𝑙 ( 1 − 𝜈 ) 𝑑𝑵 1 𝑇 𝑑𝑵 1 ⎫ 𝑑𝑵 1 𝑇 𝑑𝑵 1 + 𝑒𝑥 𝑙𝑒𝑥 𝑑𝜉 𝑑𝜉 𝑙𝑒𝑦 2 𝑑𝜂 𝑑𝜂 ⎪ ⎪ ( )𝑇 ( ) ( ( ) ( ) 𝑇 ) 𝑑𝑵 ⎪ 𝑑𝑵 1 𝑑𝑵 1 𝑑𝑵 1 1−𝜈 1 ⎪ =𝜈 + 𝑑𝜉 𝑑𝜂 2 𝑑𝜂 𝑑𝜉 ⎪ ⎬ (45) ( ) ( ) ( ) ( 𝑑𝑵 )𝑇 ( 𝑑𝑵 ) 𝑑𝑵 1 𝑇 𝑑𝑵 1 1−𝜈 ⎪ 1 1 =𝜈 + ⎪ 𝑑𝜂 𝑑𝜉 2 𝑑𝜉 𝑑𝜂 ( ) ( ) ( ) ( )⎪ 𝑙𝑒𝑦 ( 1 − 𝜈 ) 𝑑𝑵 1 𝑇 𝑑𝑵 1 ⎪ 𝑙𝑒𝑥 𝑑𝑵 1 𝑇 𝑑𝑵 1 = + ⎪ 𝑙𝑒𝑦 𝑑𝜂 𝑑𝜂 𝑙𝑒𝑥 2 𝑑𝜉 𝑑𝜉 ⎭ 𝑙𝑒𝑦

𝑁 ∑ 𝑁 ∑ 𝑖=1 𝑗=1

Using the properties of Kronecker product, the stiffness matrix components 𝑲 e,1 , 𝑲 e,2 , 𝑲 e,3 and 𝑲 e,4 can be obtained, which are same as given in Eqs. (24) and (25). Further, 𝑲 𝑒𝛽 can be simplified as, [ 𝑒,1 𝑒,2 ] 𝑲𝛽 𝑲𝛽 𝑒 𝑲𝛽 = (44) 𝑲 𝑒,3 𝑲 𝑒,4 𝛽 𝛽

𝑲 𝑒,1 = 𝛽

(49)

(41)

Thus, element stiffness coefficients and hence the element displacements will become functions of random variables 𝜶 𝑒R . When 𝐸(𝐱) is modelled as a homogeneous Gaussian field, the 𝑲 𝒆 in Eq. (21) can be written as, ( ) ( ) 𝑲 𝒆 𝜶 𝒆𝑹 𝑼 𝒆 𝜶 𝒆𝑹 = 𝑭 𝒆 (42) where, [ e,1 e,2 ] 1 1( ) 𝜇𝑬 𝑲 𝑲 𝑵 2 𝜶 𝒆R 𝑲 𝒆𝛽 d𝜉d𝜂, 𝑲𝒆 = e,3 e,4 + 2 𝑲 𝑲 1 − 𝜈 ∫0 ∫0

𝑲 𝐼𝐼 𝑖𝑗 =

By applying the expectation and variance operators on the first order or second order approximation of Eq. (48), the first and second order statistics of displacement can be obtained as,

where, ⎫ 𝑹𝑒R = 𝑻 𝑒1R ⊗ 𝑻 𝑒2R ⎪ 𝑻 𝑒1R = {𝝋𝑇1R (𝜉1 )𝝋𝑇1R (𝜉2 ) … … 𝝋𝑇1R (𝜉𝑛R +1 )}𝑇 ⎬ 𝑒 𝑇 𝑇 𝑇 𝑇 𝑻 2R = {𝝋2R (𝜂1 )𝝋2R (𝜂2 ) … … 𝝋2R (𝜂𝑛R +1 )} ⎪ ⎭

𝜕𝑲 || , 𝜕𝜶 𝑖 ||𝜦=0

(

5. Numerical examples Three numerical examples are solved based on the proposed stochastic BSWI WFEM formulations for 1D and 2D problems. Once the response statistics are calculated via perturbation method, the results are compared with the statistics obtained from MCS [57,58].

Following the random field discretization and the element stiffness matrices formulated as given in Eq. (42) for all the sub-domains can be assembled to form the algebraic set of global equilibrium equations as, ( ) ( ) 𝑲 𝜶𝑹 𝑼 𝜶𝑹 = 𝑭 (46)

5.1. One-dimensional bar with a linear body force In this example, a 1D bar with a linear body force 𝑏(𝑥) = 𝑥 as shown in Fig. 4 is considered for analysis. Two cases are studied in which the modulus of elasticity is modelled as a homogeneous random field with a Gaussian distribution as given by Eq. (27); and with a lognormal distribution as given in Eq. (30). For both the cases, the results from stochastic WFEM wherein both the input random field and response are discretized using BSWI scaling functions are compared with those obtained from stochastic FEM wherein both random field and response are discretized using Lagrange shape functions. The mean value of Young’s modulus is taken as 𝜇𝐸 = 1 × 105 MPa for both the distributions, with 𝐿 = 100 mm and 𝑎𝑟𝑒𝑎 = 1 mm2 . The entire domain

From Eq. (46), to calculate the second moment characteristics of displacements the present work proposes to use the perturbation method which is discussed in the next section. 4.3. Perturbation method Perturbation method uses the expansion of the global stiffness matrix and the displacement vector via Taylor series [3,18,54]. It is based on the assumption that the variance of the random field should be 6

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996 Table 1 Normalized computational times for perturbation method and MCS using linear and quadratic WFEM and FEM shape functions, for Gaussian distribution of 𝐸(𝐱). Order of shape functions

Fig. 4. One end fixed bar subjected to linear body force 𝑏(𝑥).

Linear Quadratic a 5000

of the bar is modelled using one BSWI 𝐶0 type element. Initially, a convergence study is done to determine the order and resolution of scaling functions used for displacement approximation. Figs. 5a and 5b shows the displacement and stress obtained from deterministic WFEM analysis for different order and resolution. It is observed that converged results in displacement and axial stress can be obtained when linear (order 𝑚 = 2) and quadratic (order 𝑚 = 3) BSWI wavelet scaling functions with resolution 𝑗 = 3 are used. Hence, for further analysis the displacement field is approximated using the order 𝑚 = {2, 3} and resolution 𝑗 = 3. The minimum resolution needed to accurately capture the spatially varying Young’s modulus while approximating with the wavelet based shape functions is established by examining the mean and standard deviation of generated sample values along the bar from a MCS. From a convergence study, based on the calculation of relative percentage error in 𝐿2 norm of mean and standard deviation values of Young’s modulus for various MCS sample size; it is noted that an error of less than 1% is obtained when the MCS sample size is 5000. Hence, a sample size of 5000 is selected to be used during the process. Further, this study is extended to Lagrange shape functions also to identify the number of nodes required for random field discretization. The process is repeated for both Gaussian and non-Gaussian approximations and is explained in the following discussions.

WFEM

FEM

Perturbation

MCSa

Perturbation

MCS

4.57 11.47

53.65 129.77

1 2.46

51.46 75.58

samples.

compared with the same obtained from MCS. The values are calculated and plotted for different values of CV, obtained by varying the standard deviation of Young’s modulus 𝐸(𝐱). The correlation length parameter considered during the process is equal to 50. It can be observed that when linear shape functions are used, at a CV of 20%, FOP results underestimate the mean values of displacement when compared with SOP and MCS results by around 11%. When quadratic shape functions are used, the deviation between FOP and the rest increases to 15% at a CV of 20%. Further, it can be seen from Fig. 9 that at a CV of 20%, the standard deviation of displacements obtained from FOP and SOP deviates by around 6% in comparison with the MCS results. Hence, this observation reiterates the assumption that the perturbation approach used in Eqs. (50) and (51) is suitable for small variances of 𝐸(𝐱) for which CV is kept under 20%. For a CV of 10%, Fig. 10 shows variation of the mean values calculated at 𝑥 = 100 mm for different correlation length parameters by employing linear and quadratic shape functions for both response and random field modelling. It is observed from Fig. 10 that regardless of correlation length parameters, the mean values of displacement evaluated using FOP and SOP accord well with the MCS values for both WFEM and FEM using linear or quadratic shape functions. Fig. 11 shows the variation of the standard deviation of displacement evaluated using linear and quadratic shape functions against correlation length parameter. When the correlation length parameter is decreased to a small value of 0.1, it is observed that the deviation between the results obtained from WFEM and FEM is less than 1% when linear shape functions are used for modelling. However, the deviation between WFEM and FEM values increase to 3% when quadratic shape functions are used in the analysis. It is also observed that when the correlation length parameter is increased to 100, the deviation between the results obtained from the perturbation method and MCS (for both WFEM and FEM) increases to 1%–2%. The dependency of size of the stochastic mesh on correlation length parameter for FEM shape functions is well documented in the literature [40,59–61]. In the present study, the deviation that is observed while using stochastic FEM can be attributed to improper mapping between response mesh and random field mesh due to non-selection of appropriate mesh size for the given correlation length parameters. This reinstates the fact that using FEM shape functions will put a limitation on the selection of the random field mesh size. On the contrary, no such limitation exists when BSWI scaling functions are used and a coarse random field discretization can also be used to accurately capture the results irrespective of the correlation length parameter. Besides the evaluation of mean and standard deviation of displacement as discussed so far, the normalized computational time required for implementing the perturbation method (FOP and SOP combined) and MCS using WFEM and FEM based shape functions is also calculated and shown in Table 1. The configuration of CPU (Central Processing Unit) is kept the same throughout the computational process. While evaluating the computational time, the correlation length parameter is kept equal to 50 and CV of 10% is assumed. The sample size used for MCS is equal to 5000. It can be observed from Table 1 that perturbation approach requires less computational effort when compared with MCS. It is also noted that FEM based perturbation approach takes less time when compared with WFEM. The higher computational times needed by

Case 1: Random field 𝐸(𝐱) follows Gaussian distribution The mean values and standard deviation of a Young’s modulus with a Gaussian approximation as defined in Eq. (27) is reproduced from the MCS samples when random field 𝛼 is modelled using Eq. (33). The same has been calculated when 𝛼 is modelled with the help of linear and quadratic Lagrange shape functions too. Variation of mean values and standard deviation of Young’s modulus at 𝑥 = 100 mm for the discretization of the random field using different number of random variables is shown in Figs. 6a and 6b respectively. It can be observed that a resolution of 𝑗 = 1 with linear and quadratic scaling functions resulting in 3 and 4 nodes respectively captures the mean and standard deviation values of Young’s modulus accurately. Further, in the case of FEM based shape function, 2 linear elements and 2 quadratic elements resulting in 3 and 5 nodes respectively is sufficient to model the random field accurately. Using MCS, the values of the first row of covariance matrix are found out for various correlation length parameters using linear and quadratic WFEM and FEM based shape functions and compared with the exact values as shown in Fig. 7. A CV of 10% is considered for the study. A resolution of 𝑗 = 1 is used for the discretization of a random field using WFEM. It can be observed that the covariance matrix values calculated using WFEM and FEM are in good agreement with the exact values for all the correlation length parameters. This also proves that the shape function methods can accurately reproduce the second order statistics of a random field regardless of the correlation length parameter. Based on the above studies, for further analysis to solve the stochastic boundary value problem as described in Case 1, the displacement field and random field are approximated by BSWI scaling function having resolution 𝑗 = 3 and 𝑗 = 1 respectively. In the same way, for discretization of displacement field and random field using FEM, 8 and 2 linear elements respectively and 5 and 2 quadratic elements respectively are used. The mean values of displacement field at 𝑥 = 100 mm obtained by using the first order perturbation (FOP) approximation as in Eq. (50) and second order perturbation approximation (SOP) as given in Eq. (51) respectively are shown in Fig. 8. The values are 7

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 5. Variation of axial displacement (𝛿) and axial stress (𝜎𝑥𝑥 ) along the length of the rod obtained from a deterministic study while linear and quadratic BSWI scaling functions are used to model displacement.

Fig. 6. Variation of mean and standard deviation values of Young’s modulus at 𝑥 = 100 mm with number of nodes (random variables) while linear and quadratic BSWI scaling functions and FEM based shape functions are used to model Gaussian field.

Fig. 7. Comparison of values of first row of covariance matrix obtained using linear and quadratic BSWI WFEM and FEM based shape functions for Gaussian field, with exact values for various correlation length parameters.

WFEM could be attributed to the fact that the scaling functions of

further reduced by way of better programming practices and optimized implementation of algorithms.

WFEM are formed using B-splines, which are piecewise polynomials

Case 2: Random field 𝐸(𝐱) with lognormal distribution Here, the modulus of elasticity is modelled as a homogeneous random field with a lognormal distribution as given in Eq. (30). A similar procedure analogous to the one discussed in the context of a Gaussian distribution is followed and convergence is established for the discretization of random field. A MCS sample size of 5000 is selected. A resolution of 𝑗 = 1 using linear and quadratic wavelet scaling functions resulting in 3 and 4 nodes respectively and with FEM shape functions 2 linear elements and 2 quadratic elements resulting in 3 and 5 nodes respectively lead to converged solutions. Hence, these parameters are used in further calculations. The mean values of displacement field

and its explicit expression are obtained at Gauss points during the evaluation of system of equations. As a result, when the number of Gauss points increases, the number of function calls for obtaining the explicit expression of B-splines and thereby evaluating the scaling functions also increases. This increases the computational overhead resulting in a higher computational time. On the contrary, FEM shape functions are not piecewise polynomials due to which, this issue is not encountered. The issue of higher computational overhead as encountered in the case of WFEM can be possibly addressed and the computational times can be 8

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 8. Mean value of displacement at 𝑥 = 100 mm on the bar for different values of CV by varying standard deviation of 𝐸(𝐱) with Gaussian random field approximation, obtained by using linear and quadratic BSWI WFEM and FEM based shape functions for both response and random field modelling.

Fig. 9. Standard deviation of displacement at 𝑥 = 100 mm on the bar for different values of CV by varying standard deviation of 𝐸(𝐱) with Gaussian random field approximation, obtained by using linear and quadratic BSWI WFEM and FEM based shape functions for both response and random field modelling.

Fig. 10. Mean value of displacement at 𝑥 = 100 mm on the bar for varying correlation length parameter for 𝐸(𝐱) with Gaussian random field approximation, obtained by using linear and quadratic BSWI WFEM and FEM based shape functions for both response and random field modelling.

Fig. 11. Standard deviation of displacement at 𝑥 = 100 mm on the bar for varying correlation length parameter for 𝐸(𝐱) with Gaussian random field approximation, obtained by using linear and quadratic BSWI WFEM and FEM based shape functions for both response and random field modelling.

9

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 12. Mean value of displacement at 𝑥 = 100 mm on the bar for different values of CV by varying standard deviation of 𝐸(𝐱) with lognormal random field approximation, obtained by using linear BSWI WFEM and FEM based shape functions for both random field and response modelling.

Fig. 13. Standard deviation of displacement at 𝑥 = 100 mm on the bar for different values of CV by varying standard deviation of 𝐸(𝐱) with lognormal random field approximation, obtained by using linear BSWI WFEM and FEM based shape functions for both random field and response modelling. Table 2 Normalized computational times for Perturbation method and MCS using linear and quadratic WFEM and FEM shape functions, for lognormal distribution of 𝐸(𝐱).

obtained using FOP, SOP and MCS is shown in Fig. 12. The values are calculated at 𝑥 = 100 mm, for different values of CV obtained by varying the standard deviation of Young’s modulus 𝐸(𝐱). It can be observed from Fig. 12 that at a CV of 20%, the mean values of displacement field using FOP deviate away from SOP and MCS values by more than 4.5%. Moreover, there is a difference of 13%–15% between the MCS values and the SOP values for CV beyond 20% when WFEM (linear and quadratic shape functions) and FEM (linear shape functions) are used. Hence, the perturbation approach that is used here is suitable for CV of less than 20%. The standard deviation of displacement field at 𝑥 = 100 mm obtained for different values of CV is shown in Fig. 13. It can be observed that at a CV of 20%, SOP values deviate away from MCS values by 2.5%–3.5%. When the CV is increased beyond 20%, the deviation also increases to 7.5%–9%. The mean values of displacement field calculated for varying correlation length parameter is shown in Fig. 14. It can be seen that the values obtained from SOP deviate away from MCS values by at the most 1.5%–2% for different correlation length parameters. The variation of standard deviation of displacement against correlation length parameters for linear and quadratic shape functions is shown in Fig. 15. When linear shape functions are used for modelling the random field, the values obtained from WFEM and FEM concur well with each other. In the case of quadratic shape functions, it can be observed that there is a deviation of 2%–3% between WFEM and FEM values at small correlation length parameters. Analogous to Gaussian distribution, the normalized computational time required for implementing the perturbation method and MCS based on WFEM and FEM shape functions is also calculated for lognormal distribution and given in Table 2. The parameters such as correlation length parameter, CV of 𝐸(𝐱) and MCS sample size are all kept the same as, in the case of Gaussian distribution. From Table 2 it can be noted that MCS takes more computational effort when compared to a perturbation approach. The perturbation method based on

Order of shape functions

Linear Quadratic a

WFEM

FEM

Perturbation

MCSa

Perturbation

MCS

4.87 10.32

51.76 115.51

1 2.08

44.19 62.52

5000 samples.

FEM shape functions has the least computational time when compared with other methods, which could be attributed to the reasons already discussed in the case of Gaussian distribution. In conclusion, perturbation approach using wavelet scaling functions is found to be accurate for both Gaussian and non-Gaussian fields, when compared with MCS and FEM based evaluations. The computational times obtained from WFEM based perturbation approach are on a higher side in comparison with FEM. However, it is important to note that if randomness associated with geometric parameters are involved, re-meshing and convergence studies needs to be done; in which case WFEM due to its inherent properties like MRA could become more attractive with respect to computational time too. Further, the effect of small correlation length on convergence associated with FEM discretization also needs to be accounted for in the overall computational time. Thus, it is clear from this example that BSWI WFEM results are comparable in different aspects to FEM, for both Gaussian and non-Gaussian input fields. Hence, for further numerical examples only Gaussian input field and WFEM is used. 5.2. Two-dimensional plane stress problem In this example, a plane stress problem under uni-axial loading is studied. Due to symmetry, only a quarter portion of the plate is 10

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 14. Mean value of displacement at 𝑥 = 100 mm on the bar for varying correlation length parameter for 𝐸(𝐱) with lognormal random field approximation, obtained by using linear BSWI WFEM and FEM based shape functions for both random field and response modelling.

Fig. 15. Standard deviation of displacement at 𝑥 = 100 mm on the bar for varying correlation length parameter for 𝐸(𝐱) with lognormal random field approximation, obtained by using linear BSWI WFEM and FEM based shape functions for both random field and response modelling.

considered for analysis. Roller supports are employed at the left edge and bottom edge as shown in Fig. 16. The modulus of elasticity is modelled as a homogeneous Gaussian random field as given by Eq. (27). The mean value of Young’s modulus is taken as 𝜇𝐸 = 2 × 105 MPa, with 𝜈 = 0.3, 𝐿 = 5 mm, 𝑏 = 6 mm and 𝐹 = 1000 N∕mm2 . The analytical solution obtained from deterministic analysis is given in Sadd [62]. The entire domain is modelled using one 2D BSWI 𝐶0 type element. In this problem, the displacement and random field are approximated using linear and quadratic shape functions and the results are studied for each case separately. For the current study to approximate the displacement field, scaling functions with 𝑚 = 2 and 𝑚 = 3 with a resolution 𝑗 = 3 is considered. Further, random field is approximated with scaling function of same order that is used for displacement approximation but having a coarse resolution of 𝑗 = 1 to reduce the computational effort. For illustration, the nodal distribution used for displacement and random field approximation using linear wavelet scaling functions is shown in Figs. 17(a) and 17(b) respectively. The mean and standard deviation values of displacement in 𝑦 direction (𝜇𝑣 and 𝜎𝑣 ), along the path A to B (𝑥 = 2.5 mm, 𝑦 = [0 ∶ 6] mm) as seen in Fig. 16, is calculated using FOP, SOP and MCS for both linear and quadratic WFEM scaling functions. During this analysis, the spatially varying Young’s modulus 𝐸(𝐱) is considered having a variance of magnitude 2 × 108 which leads to a CV of 7.07%. The correlation length parameters 𝑐1 , 𝑐2 in 𝑥 and 𝑦 directions are considered equal to 2.5 and 3 respectively. Table 3 shows the mean and standard deviation values of displacement obtained when linear scaling functions and quadratic scaling functions are employed. It can be observed from Table 3 that the mean and standard deviation values obtained using FOP and SOP concur will with the MCS values for linear as well as quadratic shape functions. The variation of mean values of displacement in 𝑦 direction at position C (𝑥 = 5 mm, 𝑦 = 6 mm) for different CV is shown in Fig. 18a. It can be observed from Fig. 18 that mean value obtained from perturbation approach are in good agreement with MCS values, for the

Fig. 16. Uni-axial loading of a plate under plane stress.

range of CV considered. The variation of the standard deviation values of displacement in 𝑦 direction at position C are shown in Fig. 18b. The difference in the standard deviation of perturbation results from MCS is within 5% as observed for different values of CV. Thus, it can be concluded that the mean values can be accurately captured using the WFEM based perturbation approach for a 2D plane stress problem. It is also observed that when the CV is increased beyond 12%, the standard deviation values obtained from WFEM deviate by not more than 5% in comparison with MCS values. In conclusion, 11

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Fig. 17. Nodal distribution used for displacement and random field approximation when the domain is discretized using single 2D BSWI 𝐶0 type element constructed using linear WFEM scaling functions.

Fig. 18. Mean value and standard deviation of displacement in 𝑦 direction at position C, for different values of CV by varying standard deviation of 𝐸(𝐱), obtained using linear and quadratic BSWI WFEM based scaling functions.

Table 3 Mean and standard deviation values of displacement in 𝑦 direction (𝜇𝑣 and 𝜎𝑣 ) along path A to B obtained by using linear and quadratic WFEM scaling functions. Linear/Quadratic

Distance along y (mm)

𝜇𝑉

𝜎𝑉

FOP

SOP

MCS

FOP

SOP

MCS

Linear 𝑚 = 2

0 1.50 3.00 4.50 6.00

0 0.00750 0.01500 0.02250 0.03000

0 0.00751 0.01503 0.02255 0.03008

0 0.00752 0.01503 0.02254 0.03006

0 0.000365 0.000682 0.001000 0.001286

0 0.000366 0.000683 0.001002 0.001289

0 0.00036 0.00069 0.00101 0.00131

Quadratic 𝑚 = 3

0 1.33 2.00 3.33 4.66 6.00

0 0.00665 0.00992 0.01668 0.02335 0.03000

0 0.00671 0.00999 0.01675 0.02343 0.03016

0 0.00672 0.01000 0.01675 0.02343 0.03015

0 0.00060 0.00072 0.00088 0.00106 0.00153

0 0.00061 0.00072 0.00088 0.00106 0.00154

0 0.00062 0.00073 0.00089 0.00106 0.00149

of Young’s modulus is taken as 𝜇𝐸 = 2 × 105 MPa, Poisson’s ratio is considered to be deterministic with a value, 𝜈 = 0.3 and a far-field stress of 𝐹 = 1000 N∕mm2 is applied. The 2D irregular domain of this problem is modelled using two 2D BSWI 𝐶0 type elements as it is the minimum number of quadrilateral elements needed to discretize the geometry. The analytical solution of the deterministic analysis is given in Sadd [62]. The results of convergence study using a deterministic approach showing the stress concentration factor along 𝑦 = 0 for various order and resolutions of wavelet scaling functions are shown in [63] and it is observed that 𝑚 = 3, 𝑗 = 3 leads to a converged solution. Hence, in the current stochastic approach, the displacement field is approximated using 𝑚 = 3, 𝑗 = 3; but random field is approximated

perturbation results obtained from WFEM concur well with the results of MCS for a 2D plane stress problem with only a small margin of error at higher CV values. 5.3. Plate with a hole under uni-axial far field tension In this example, a plate with a hole under uni-axial far field tensile loading is studied. Due to symmetry, a quarter portion of the plate with dimensions, 𝐿 = 5 mm, 𝑏 = 6 mm, and radius of the hole, 𝑎 = 0.1 mm having symmetry boundary conditions at the left edge and bottom edge as shown in Fig. 19 is considered for analysis. The modulus of elasticity is modelled as a homogeneous Gaussian random field. The mean value 12

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

used for displacement and random field approximation is shown in Figs. 20(a) and 20(b) respectively. A variance of magnitude 2 × 108 which leads to a CV of 7.07% is considered. Further, in the present study the correlation length parameters 𝑐1 , 𝑐2 in 𝑥 and 𝑦 directions are considered equal to 2.5 and 3 respectively. The mean and standard deviation values of displacement in 𝑥 direction, 𝜇𝑢 and 𝜎𝑢 respectively, along 𝑦 = 0 are plotted in Fig. 21. It can be seen that the values obtained from perturbation approach are in good agreement with MCS values. Further, for different values of CV (5% to 14.14%), mean and standard deviation of displacement are calculated at various positions A and B on the problem domain. Table 4 shows the perturbation and MCS results of the mean and standard deviation of displacement at position A (𝑥 = 0.1 mm, 𝑦 = 0 mm). It can be observed from Table 4 that the perturbation and MCS results match well even for the highest CV considered in the study. The difference between perturbation and MCS is observed to be less than 2%. A similar observation can be made for the mean and standard deviation of displacement at position B (𝑥 = 0 mm, 𝑦 = 0.1 mm) which is tabulated in Table 4. In conclusion, it can be observed that accurate results are obtained when BSWI scaling functions at lower order and resolution, resulting in less number of random variables, are used for modelling the random field. Hence, with the proposed stochastic WFEM, a good accuracy can be achieved.

Fig. 19. Far field tensile loading of a plate with a hole.

using a coarser discretization of 𝑚 = 2, 𝑗 = 1 for both the elements to reduce the computational effort. For illustration, the nodal distribution

Fig. 20. Nodal distribution used for displacement and random field approximation when the domain is discretized using two 2D BSWI 𝐶0 type elements.

Fig. 21. Mean (𝜇𝑢 ) and standard deviation (𝜎𝑢 ) of displacement in 𝑥 direction, along path A to D. 13

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996

Table 4 Mean and standard deviation values of displacement in 𝑥 direction (𝜇𝑢 and 𝜎𝑢 ) at position A and in 𝑦 direction (𝜇𝑣 and 𝜎𝑣 ) at position B, for different values of CV obtained by varying standard deviation of 𝐸(𝐱). Position

CV %

𝜇𝑈 FOP

𝜎𝑈

SOP

A

5.000 7.070 10.00 14.14

−3.822 −3.822 −3.822 −3.822

B

5.000 7.070 10.00 14.14

1.195 1.195 1.195 1.195

× × × ×

10−4 10−4 10−4 10−4

MCS

−3.822 −3.822 −3.821 −3.821

× × × ×

10−4 10−4 10−4 10−4

FOP

−3.822 −3.821 −3.820 −3.818

× × × ×

10−4 10−4 10−4 10−4

SOP × × × ×

10−6 10−6 10−5 10−5

6.94 9.83 1.39 1.97

× × × ×

1.12 1.58 2.24 3.17

× × × ×

10−5 10−5 10−5 10−5

1.12 1.58 2.25 3.19

× × × ×

𝜇𝑉 × × × ×

10−3 10−3 10−3 10−3

1.195 1.196 1.196 1.198

× × × ×

MCS

6.93 9.81 1.38 1.96

10−6 10−6 10−5 10−5

6.94 9.83 1.40 2.01

× × × ×

10−6 10−6 10−5 10−5

1.12 1.58 2.26 3.26

× × × ×

10−5 10−5 10−5 10−5

𝜎𝑉

10−3 10−3 10−3 10−3

1.195 1.196 1.197 1.199

6. Conclusion

× × × ×

10−3 10−3 10−3 10−3

10−5 10−5 10−5 10−5

[7] F. Yamazaki, M. Shinozuka, Safety analysis of stochastic finite element systems by Monte Carlo simulation, Doboku Gakkai Ronbunshu (1988) 109–119, http: //dx.doi.org/10.2208/jscej.1988.398_109. [8] M. Shinozuka, G. Deodatis, Response variability of stochastic finite element systems, J. Eng. Mech. 114 (1988) 499–519, http://dx.doi.org/10.1061/(ASCE) 0733-9399(1988)114:3(499). [9] P.D. Spanos, R. Ghanem, Stochastic finite element expansion for random media, J. Eng. Mech. 115 (1989) 1035–1053, http://dx.doi.org/10.1061/(ASCE)07339399(1989)115:5(1035). [10] T. Takada, Weighted integral method in stochastic finite element analysis, Probab. Eng. Mech. 5 (1990) 146–156, http://dx.doi.org/10.1016/02668920(90)90006-6. [11] T. Takada, Weighted integral method in multi-dimensional stochastic finite element analysis, Probab. Eng. Mech. 5 (1990) 158–166, http://dx.doi.org/10. 1016/0266-8920(90)90016-D. [12] G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probab. Eng. Mech. 25 (2010) (2010) 183–197, http://dx.doi.org/10.1016/j.probengmech.2009.10.003. [13] M. Manjuprasad, C.S. Manohar, Adaptive random field mesh refinements in stochastic finite element reliability analysis of structures, Comput. Model. Eng. Sci. 19 (2007) 23–54. [14] G. Deodatis, Weighted integral method. i: Stochastic stiffness matrix, J. Eng. Mech. 117 (1991) 1851–1864, http://dx.doi.org/10.1061/(ASCE)07339399(1991)117:8(1851). [15] G. Deodatis, M. Shinozuka, Weighted.Integral. Method, Weighted integral method II: Response variability and reliability, J. Eng. Mech. 117 (1991) 1865–1877, http://dx.doi.org/10.1061/(ASCE)0733-9399(1991)117:8(1865). [16] J. Dolbow, T. Belytschko, An introduction to programming the meshless element free Galerkin method, Arch. Comput. Methods Eng. 5 (1998) 207–241, http: //dx.doi.org/10.1007/BF02897874. [17] G.R. Liu, Y.T. Gu, An Introduction To Meshfree Methods and their Programming, Springer-Verlag, Berlin/Heidelberg, 2005, http://dx.doi.org/10.1007/14020-3468-7. [18] S. Rahman, B.N. Rao, A perturbation method for stochastic meshless analysis in elastostatics, Internat. J. Numer. Methods Engrg. 50 (2001) 1969–1991, http://dx.doi.org/10.1002/nme.106. [19] S. Rahman, H. Xu, A meshless method for computational stochastic mechanics, Int. J. Comput. Methods Eng. Sci. Mech. 6 (2005) 41–58, http://dx.doi.org/10. 1080/15502280590888649. [20] B. Li, X. Chen, Wavelet-based numerical analysis: A review and classification, Finite Elem. Anal. Des. 81 (2014) 14–31, http://dx.doi.org/10.1016/j.finel.2013. 11.001. [21] I. Daubechies, Orthonormal bases of compactly supported wavelets, Comm. Pure Appl. Math. 41 (1988) 909–996, http://dx.doi.org/10.1002/cpa.3160410705. [22] C.K. Chui, An Introduction To Wavelets - Wavelet Analysis and Its Applications, Academic press, 1992. [23] I. Daubechies, Ten lectures on wavelets, in: IEEE Symp. Comput. Med. Syst. Society for Industrial and Applied Mathematics, 1992, p. 357. http://dx.doi.org/ 10.1137/1.9781611970104. [24] K. Amaratunga, J.R. Williams, S. Qian, J. Weiss, Wavelet-Galerkin solutions for one-dimensional partial differential equations, Internat. J. Numer. Methods Engrg. 37 (1994) 2703–2716, http://dx.doi.org/10.1002/nme.1620371602. [25] C.K. Chui, E. Quak, Wavelets on a Bounded Interval, in: Numer. Methods Approx. Theory, vol. 9, Birkhäuser Basel, Basel, 1992, pp. 53–75, http://dx.doi.org/10. 1007/978-3-0348-8619-2_4. [26] W.H. Chen, C.W. Wu, A spline wavelets element method for frame structures vibration, Comput. Mech. 16 (1995) 11–21, http://dx.doi.org/10.1007/ BF00369881. [27] C.W. Wu, W.H. Chen, Extension of spline wavelets element method to membrane vibration analysis, Comput. Mech. 18 (1996) 46–54, http://dx.doi.org/10.1007/ BF00384175.

In the current paper, the formulation of stochastic BSWI WFEM for 1D and 2D problems in elasto-statics is proposed in which the spatial variation of modulus of elasticity is modelled as a homogeneous random field. The present work suggests using BSWI scaling functions to discretize random field. The response statistics of mean and standard deviation are calculated using the perturbation approach. Numerical examples based on 1D and 2D elasto-statics are solved and the results obtained from WFEM based perturbation approach are compared with the results of MCS. Further, in order to validate the proposed methodology perturbation and MCS results based on FEM Lagrange shape functions for 1D problem are also obtained and a comparative study between the results of WFEM and FEM is done. It can be concluded from the results of 1D and 2D numerical examples that the proposed BSWI WFEM based perturbation approach (both FOP and SOP) accurately captures the response statistics of the displacement field for CV values of up to 15%–20%. MRA properties of BSWI wavelets help the physical problem to be discretized using less number of elements; say 1 or 2, as witnessed in the numerical examples and thus can reduce the computational cost required for mesh assembly when compared to FEM. A single BSWI WFEM element with coarse nodal distribution for random field modelling, perform well by accurately capturing the standard deviations even at extremely small or large correlation length parameters as observed in 1D problem. On the other hand, FEM demands an appropriate response and random field mesh size in accordance with the correlation length of the input field, which otherwise would result in an erroneous result. Normalized computational times are also shown for 1D problems, from which it can be concluded that FEM based perturbation times are less in comparison with WFEM. However, it is to be noted that those computational times are estimated without accounting for the time associated with re-meshing and convergence studies. The 2D plane stress problems demonstrate that high accuracy can be maintained even when the random field is modelled with wavelet scaling functions of lower order and lower resolution, thereby reducing the computational effort needed. References [1] S.-K. Choi, R.V. Grandhi, R.A. Canfield, Reliability-Based Structural Design, Springer London, 2007, http://dx.doi.org/10.1007/978-1-84628-445-8. [2] A. Haldar, S. Mahadevan, Probability, Reliability, and Statistical Methods in Engineering Design, John Wiley & Sons, 1999. [3] G. Stefanou, The stochastic finite element method: Past, present and future, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1031–1051, http://dx.doi.org/ 10.1016/j.cma.2008.11.007. [4] E. Vanmarcke, M. Grigoriu, Stochastic finite element analysis of simple beams, J. Eng. Mech. 109 (1983) 1203–1214, http://dx.doi.org/10.1061/(ASCE)07339399(1983)109:5(1203). [5] W.K. Liu, T. Belytschko, A. Mani, Random field finite elements, Internat. J. Numer. Methods Engrg. 23 (1986) 1831–1845, http://dx.doi.org/10.1002/nme. 1620231004. [6] F. Yamazaki, A. Member, M. Shinozuka, G. Dasgupta, Neumann expansion for stochastic finite element analysis, J. Eng. Mech. 114 (1988) 1335–1354, http://dx.doi.org/10.1061/(ASCE)0733-9399(1988)114:8(1335). 14

S. Vadlamani and Arun C.O.

Probabilistic Engineering Mechanics 58 (2019) 102996 [45] R. Andreev, Implementation of Sparse Wavelet-Galerkin FEM for Stochastic PDEs, 2009. [46] X. Zhang, X. Chen, Z. Yang, B. Li, Z. He, A stochastic wavelet finite element method for 1D and 2D structures analysis, Shock Vib. 2014 (2014) 1–15, http://dx.doi.org/10.1155/2014/104347. [47] C.K. Chui, J.-Z. Wang, On compactly supported spline wavelets and a duality principle, Trans. Amer. Math. Soc. 330 (1992) 903–915, http://dx.doi.org/10. 1090/S0002-9947-1992-1076613-3. [48] C.K. Chui, J. Wang, A cardinal spline approach to wavelets, Proc. Amer. Math. Soc. 113 (1991) 785, http://dx.doi.org/10.1090/S0002-9939-1991-1077784-X. [49] C.K. Chui, J.Z. Wang, An analysis of cardinal spline-wavelets, J. Approx. Theory. 72 (1993) 54–68, http://dx.doi.org/10.1006/jath.1993.1006. [50] S. Bertoluzza, G. Naldi, J.C. Ravel, Wavelet methods for the numerical solution of boundary value problems on the interval, in: C.K. Chui, et al. (Eds.), Wavelets Theory, Algorithms, Appl., 1994, pp. 425–448, http://dx.doi.org/10.1016/B9780-08-052084-1.50024-7. [51] E. Quak, N. Weyrich, Decomposition and reconstruction algorithms for spline wavelets on a bounded interval, Appl. Comput. Harmon. Anal. 1 (1994) 217–231, http://dx.doi.org/10.1006/acha.1994.1009. [52] J.C. Goswami, A.K. Chan, C.K. Chui, On solving first-kind integral equations using wavelets on a bounded interval, IEEE Trans. (1995) http://dx.doi.org/10. 1109/8.387178. [53] H. Xu, S. Rahman, Decomposition methods for structural reliability analysis, Probab. Eng. Mech. 20 (2005) 239–250, http://dx.doi.org/10.1016/j. probengmech.2005.05.005. [54] A. Gupta, C.O. Arun, Stochastic meshfree method for elastic buckling analysis of columns, Comput. Struct. 194 (2018) 32–47, http://dx.doi.org/10.1016/j. compstruc.2017.08.014. [55] C.O. Arun, B.N. Rao, S.M. Srinivasan, Stochastic meshfree method for elastoplastic damage analysis, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2590–2606, http://dx.doi.org/10.1016/j.cma.2010.04.009. [56] S. Rahman, B.N. Rao, An element-free Galerkin method for probabilistic mechanics and reliability, Int. J. Solids Struct. 38 (2001) 9313–9330, http://dx.doi.org/ 10.1016/S0020-7683(01)00193-7. [57] R.Y. Rubinstein, D.P. Kroese, Simulation and the Monte Carlo Method, second ed., Wiley, 2007. [58] P. Glasserman, Monte Carlo Methods in Financial Engineering, Springer, New York, 2003, http://dx.doi.org/10.1007/978-0-387-21617-1. [59] A. Der Kiureghian, J.-B. Ke, The stochastic finite element method in structural reliability, Probab. Eng. Mech. 3 (1988) 83–91, http://dx.doi.org/10.1016/02668920(88)90019-7. [60] C. Li, A. Der Kiureghian, Optimal discretization of random fields, J. Eng. Mech. 119 (1993) 1136–1154, http://dx.doi.org/10.1061/(ASCE)0733-9399(1993)119: 6(1136). [61] B.A. Zeldin, P.D. Spanos, On random field discretization in stochastic finite elements, J. Appl. Mech. 65 (1998) 320, http://dx.doi.org/10.1115/1.2789057. [62] M.H. Sadd, Elasicity - Theory, Applications, and Numerics, Elsevier Academic Press, 2005. [63] S. Vadlamani, Arun. C.O., A background cell-based numerical integration for B-spline wavelet on the interval finite element method, Eng. Comput. 36 (2019) 569–598, http://dx.doi.org/10.1108/EC-07-2018-0315.

[28] J.W. Xiang, X.F. Chen, Z.J. He, H.B. Dong, The construction of 1D wavelet finite elements for structural analysis, Comput. Mech. 40 (2007) 325–339, http: //dx.doi.org/10.1007/s00466-006-0102-5. [29] J. Xiang, X. Chen, Y. He, Z. He, The construction of plane elastomechanics and mindlin plate elements of B-spline wavelet on the interval, Finite Elem. Anal. Des. 42 (2006) 1269–1280, http://dx.doi.org/10.1016/j.finel.2006.06.006. [30] J.-W. Xiang, M. Liang, Y.-T. Zhong, Computation of stress intensity factors using wavelet-based element, J. Mech. 32 (2016) N1–N6, http://dx.doi.org/10.1017/ jmech.2016.2. [31] M. Liu, J. Xiang, Y. Zhong, Band structures analysis method of two-dimensional phononic crystals using wavelet-based elements, Crystals 7 (2017) 328, http: //dx.doi.org/10.3390/cryst7110328. [32] G.W. Wornell, A Karhunen–Loeve-like expansion for 1/f processes via wavelets, IEEE Trans. Inf. Theory. 36 (1990) 859–861, http://dx.doi.org/10.1109/18. 53745. [33] B.A. Zeldin, P.D. Spanos, Random field representation and synthesis using wavelet bases, J. Appl. Mech. 63 (1996) 946, http://dx.doi.org/10.1115/1. 2787251. [34] A. Latto, H. Resnikoff, E. Tenenbaum, The evaluation of connection coefficients of compactly supported wavelets, in: Proc. French-USA Work. Wavelets Turbul. Princet. Univ., 1991, Springer-Verlag, 1991, pp. 76–89. [35] J. Ma, J. Xue, S. Yang, Z. He, A study of the construction and application of a Daubechies wavelet-based beam element, Finite Elem. Anal. Des. 39 (2003) 965–975, http://dx.doi.org/10.1016/S0168-874X(02)00141-5. [36] K.K. Phoon, S.P. Huang, S.T. Quek, Implementation of Karhunen–Loeve expansion for simulation using a wavelet-Galerkin scheme, Probab. Eng. Mech. 17 (2002) 293–303, http://dx.doi.org/10.1016/S0266-8920(02)00013-9. [37] C. Proppe, Reliability analysis with stochastic finite elements and wavelet approximations, in: 8th. World Congr. Comput. Mech. Venice, Italy, 2008. [38] C. Proppe, Multiresolution analysis for stochastic finite element problems with wavelet-based Karhunen-Loève expansion, Math. Probl. Eng. 2012 (2012) 1–15, http://dx.doi.org/10.1155/2012/215109. [39] K.K. Phoon, H.W. Huang, S.T. Quek, Comparison between Karhunen–Loeve and wavelet expansions for simulation of Gaussian processes, Comput. Struct. 82 (2004) 985–991, http://dx.doi.org/10.1016/j.compstruc.2004.03.008. [40] B. Sudret, A. Der Kiureghian, Stochastic Finite element methods and reliability: A State-of-the-Art Report, 2000, www.ethz.ch/content/dam/ethz/specialinterest/baug/ibk/risk-safety-and-uncertainty-dam/publications/reports/SFEreport-Sudret.pdf. [41] J.M. Angulo, M.D. Ruiz-Medina, Multi-resolution approximation to the stochastic inverse problem, Adv. Appl. Probab. 31 (1999) 1039–1057, http://www.jstor. org/stable/1428343. [42] O.P. Le Maıtˆre, O.M. Knio, H.N. Najm, R.G. Ghanem, Uncertainty propagation using Wiener–Haar expansions, J. Comput. Phys. 197 (2004) (2004) 28–57, http://dx.doi.org/10.1016/j.jcp.2003.11.033. [43] P.D. Spanos, J. Tezcan, P. Tratskas, Stochastic processes evolutionary spectrum estimation via harmonic wavelets, Comput. Methods Appl. Mech. Engrg. 194 (2005) 1367–1383, http://dx.doi.org/10.1016/j.cma.2004.06.039. [44] J.-G. Han, W.-X. Ren, Y. Huang, A wavelet-based stochastic finite element method of thin plate bending, Appl. Math. Model. 31 (2007) 181–193, http: //dx.doi.org/10.1016/j.apm.2005.08.020.

15