Faster computation of the Karhunen–Loève expansion using its domain independence property

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145 www.elsevier.com/locate/cma

Faster computation of the Karhunen–Lo`eve expansion using its domain independence property Srikara Pranesh, Debraj Ghosh ∗ Department of Civil Engineering, Indian Institute of Science, Bangalore 560012, India Received 14 June 2014; received in revised form 27 September 2014; accepted 31 October 2014 Available online 13 November 2014

Abstract The goal of this work is to reduce the cost of computing the coefficients in the Karhunen–Lo`eve (KL) expansion. The KL expansion serves as a useful and efficient tool for discretizing second-order stochastic processes with known covariance function. Its applications in engineering mechanics include discretizing random field models for elastic moduli, fluid properties, and structural response. The main computational cost of finding the coefficients of this expansion arises from numerically solving an integral eigenvalue problem with the covariance function as the integration kernel. Mathematically this is a homogeneous Fredholm equation of second type. One widely used method for solving this integral eigenvalue problem is to use finite element (FE) bases for discretizing the eigenfunctions, followed by a Galerkin projection. This method is computationally expensive. In the current work it is first shown that the shape of the physical domain in a random field does not affect the realizations of the field estimated using KL expansion, although the individual KL terms are affected. Based on this domain independence property, a numerical integration based scheme accompanied by a modification of the domain, is proposed. In addition to presenting mathematical arguments to establish the domain independence, numerical studies are also conducted to demonstrate and test the proposed method. Numerically it is demonstrated that compared to the Galerkin method the computational speed gain in the proposed method is of three to four orders of magnitude for a two dimensional example, and of one to two orders of magnitude for a three dimensional example, while retaining the same level of accuracy. It is also shown that for separable covariance kernels a further cost reduction of three to four orders of magnitude can be achieved. Both normal and lognormal fields are considered in the numerical studies. c 2014 Elsevier B.V. All rights reserved. ⃝ Keywords: Stochastic finite element; Random field; KL expansion; Integral equation; Numerical integration

1. Introduction While newer mathematical and computational tools are regularly being developed for studying engineering systems of increasing complexity and size, reliability of the predictions from these tools can be improved by systematically considering the underlying uncertainties. This consideration has led to the development of the area of uncertainty quantification (UQ). In a probability theory based UQ problem the uncertain parameters are modeled as random quantities such as variables, vectors, and processes or fields. For structural mechanics problems such parameters ∗ Corresponding author. Tel.: +91 080 2293 2818.

E-mail address: [email protected] (D. Ghosh). http://dx.doi.org/10.1016/j.cma.2014.10.053 c 2014 Elsevier B.V. All rights reserved. 0045-7825/⃝

126

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

are loads such as earthquake or winds [1,2], material properties such as Young’s modulus and density [3–5], support conditions and geometric properties [6]. These uncertainties are then propagated through methods such as Monte Carlo simulation [7], perturbation [8], and spectral stochastic finite element method [9], to name a few. For computational purpose often the random process (or field) models need to be discretized. Accordingly, a tensor product space is first defined on the cartesian product of the underlying probability space (Ω , F, P) and the physical domain D. For spatial processes, D is a closed interval in the space Rd , where d is 1, 2, or 3. Then a process κ(x, θ ) is approximated using a set of bases in this product space, where x ∈ D and θ ∈ Ω . This discretization or approximation is finite dimensional, and a lower dimensionality helps in optimizing the cost of uncertainty propagation. In the literature a number of discretization techniques such as midpoint [10,11], local (spatial) averaging [12,13], series expansions such as polynomial chaos (PC) [9], Karhunen–Lo`eve (KL) expansion [14], optimal linear estimation [15], and weighted integral [16] have been proposed. An attempt to unify some of these techniques has been made in [17]. When the mean and the covariance of the process are known, the KL expansion can be used to discretize the process, and this representation is optimal in the sense of dimensionality of the discretized process [15]. When the mean κ(x, ¯ θ ) and the continuous covariance kernel C(x1 , x2 ) are known, a square-integrable process κ(x, θ ) can be discretized in its KL expansion as κ(x, θ ) = ω(x, θ ) ≡ κ(x, ¯ θ) +

∞  

λi φi (x)ξi (θ ).

(1)

i=1

Here ξi (θ ) are zero-mean unit-variance uncorrelated random variables, λi and φi (x) are eigenvalues and eigenfunctions, respectively, of the covariance kernel satisfying the integral eigenvalue problem (IEVP)  C(x1 , x2 )φi (x1 )dx1 = λi φi (x2 ). (2) D

This integral equation is a homogeneous Fredholm equation of second type. Symmetry and positive semi-definiteness of a covariance kernel ensure that the eigenvalues are non-negative and the eigenfunctions are real-valued. Eq. (1) is convergent in mean-square sense. For computational purpose a finite-dimensional truncated version of the KL expansion is adopted as κ(x, ˆ θ ) = ω(x, ˆ θ ) ≡ κ(x, ¯ θ) +

L  

λi φi (x)ξi (θ ).

(3)

i=1

The truncation level L depends upon the correlation length of the process; a higher correlation length leads to a faster decay in the eigenvalues, thus a lower L. In many practical cases an exponentially decaying covariance kernel is used, which often needs to be modified [18] to remove the non-differentiability at the origin and thus improving the smoothness of the kernel. Although the KL expansion is used mostly for Gaussian processes, a few non-Gaussian distributions such as lognormal or beta distribution also have been modeled successfully [19–21]. Exact solutions of Eq. (2) are available only for a few specific types of covariance kernels over simple domains [14]. For most practical cases a numerical method is required to solve it. The current work is regarding this numerical method. Two successful methods for solving this IEVP are the finite element based Galerkin method (FE) [9] and numerical quadrature such as Nystr¨om method [22]. In the structural mechanics community FE is the most popular choice since the FE mesh developed for solving the mechanics problem can be re-used to solve the IEVP [9,23]. The FE approximation followed by a Galerkin projection leads to a generalized algebraic eigenvalue problem. A number of works have been conducted in this regard for accelerating the computation. In [24] an adaptive FE scheme was followed for stationary processes. In [25] the aforementioned algebraic eigenvalue problem was solved using a generalized fast multipole method. Thus the cubic complexity of the computational cost was reduced to log-linear. In [26] the eigenvalue problem was solved using a hierarchical sparse matrix approximation. Other modifications on the KL expansion have also been proposed. In [27] a Fourier–Karhunen–Lo`eve representation is followed where a truncated Fourier approximation of a wide sense stationary process is followed by a principal component analysis (PCA). Although the convenience of re-using the mechanics mesh in solving the IEVP is an advantage, one important concern naturally arises here. Should the optimal mesh for the mechanics problem be accurate and optimal for discretizing the random field as well? This issue has already been raised in the literature [10,17,23], where it was found

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

127

that such optimality does not hold. Similar mesh resolution issue appears in other discretization techniques as well [15,28]. A few alternate ways have been explored and reported in the literature. In [29] it was found numerically that if the eigenvectors are known to have a periodic behavior, then a Rayleigh–Ritz method with trigonometric basis set is better than a collocation scheme. In [30] a numerical parametric study was conducted on the factors affecting the convergence of the KL expansion, and the requirement of an efficient numerical scheme to solve Eq. (2) was highlighted. In this context another related question can be asked as follows. While specifying the characteristics of the random process only the mean, variance, and probability distribution are mentioned, the exact shape of the physical domain D is not required. However, solutions of the eigenvalue problem stated in Eq. (2) would be dependent upon D. Therefore, should the KL approximated process depend upon D? Motivated by these questions, here we proceed as follows. First, using Mercer’s theorem it is mathematically argued that the KL expansion for a process is physical domain independent. That is, although a modification in the domain D changes the eigenvalues and eigenfunctions of the covariance kernel, it does not change the realizations of the process estimated using the KL expansion. It is also shown that the truncated KL expansion should include the most dominant terms to ensure this invariance. Based on this domain independence property, a new method is proposed here. Accordingly, a complicated domain is replaced by a bounding box such as a rectangle in two-dimensions and a box in three-dimensions, over which the IEVP is solved using a numerical quadrature. Through two numerical studies – one two- and one three-dimensional – it is demonstrated that the proposed method reduces the computational cost significantly compared to the FE approach. Both Gaussian and lognormal processes are modeled. It is further demonstrated that for separable kernels the cost can be further reduced by using the separability. A set of preliminary numerical studies on the proposed approach applied to a two-dimensional problem was reported in [31]. The approach of embedding a complicated domain into a simpler domain is also used in [32], based on the finite cell method in mechanics [33]. In the current work, the domain-independence of the KL-approximated random field is proved (in Proposition 1), and thus it is established that any numerical method can be used on a simpler domain, without the need for a local mesh refinement as in [32]. The paper is organized as follows. In Section 2 a mathematical proof for the domain independence of the KL expansion for a Gaussian process is presented. The role of truncation is also discussed in this context. Based on the results in this section the proposed method is presented in Section 3. The domain independence is numerically demonstrated and the proposed method is applied and tested on two numerical problems in Section 4. A comparative study of the proposed method along with the conventional Galerkin method is performed in this section. Finally the concluding remarks are made. 2. Domain independence of the KL expansion Following the motivations mentioned in the previous section, now it will be demonstrated that if the domains D and D′ are such that D∩D′ ̸= φ, then the mean and second moments of the process remain unaltered. That is, at a common spatial location of two domains x ∈ D∩D′ , the mean and variance remain the same. Furthermore, the covariance of the process between two spatial locations x, y ∈ D∩D′ also remains unaltered. This domain independence property, when proved, would imply that if the process is completely characterized by first two moments – such as Gaussian processes – the marginal and joint probability density functions (pdfs) remain the same on the common spatial locations of D and D′ . The fundamental tools used to prove these domain independence properties are (i) the fact that the KL expansion is a second order representation of a random process, and (ii) Mercer’s theorem. A few definitions will be stated first, followed by the statement of Mercer’s theorem and KL theorem for the completeness of the discussion. Definition 1. Let D ⊂ Rd be a bounded domain and let k : D × D → R.  If |k(x, y)|2 dxdy < ∞

(4)

D

then k is referred to as a Hilbert–Schmidt kernel. Definition 2. Let k(x, y) be a Hilbert–Schmidt kernel and u(y) ∈ L 2 (D). An integral operator K : L 2 (D) → L 2 (D) defined as  [K u](x) = k(x, y)u(y)dy (5) D

is called as a Hilbert–Schmidt operator.

128

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Mercer’s Theorem. Let k : D × D → R be a continuous Hilbert–Schmidt kernel. Let the corresponding Hilbert–Schmidt operator K : L 2 (D) → L 2 (D) be positive definite, and self-adjoint with the ith eigenvalue and eigenfunction denoted as λi and φi , respectively. Then k(x, y) =

∞ 

λi φi (x)φi (y) ∀x, y ∈ D.

(6)

i=1

The convergence of this series is absolute and uniform. Karhunen–Lo`eve Theorem. Consider a centered, mean square continuous stochastic process κ : D × Ω → R, which is also square integrable, that is, κ ∈ L 2 (D × Ω ). Then there exist a basis {φi (x)} of L 2 (D) and random variables ξi (θ ) such that n  

L 2 (Ω )

λi ξi (θ )φi (x) −−−→ κ

(7)

i=1

where the convergence is in L 2 (Ω ) sense, that is, lim ∥κ(x, θ ) −

n→∞

Here ξi (θ ) =

n  

λi ξi (θ )φi (x)∥ L 2 (Ω ) → 0

∀ x ∈ D.

(8)

i=1 √1 λi



D

κ(x, θ )φi (x)dx, holding the properties

1. E[ξi (θ )] = 0 2. E[ξi (θ )ξ j (θ )] = δi j 3. V ar [ξi (θ )] = 1 where E[·] denotes the expectation operator with respect to the probability space (Ω , F, P) and V ar (·) denotes the variance. Upon using Mercer’s theorem in the Karhunen–Lo`eve theorem, it is found that λi -s and φi -s are eigenvalues and corresponding eigenfunctions of the covariance kernel of the random process, that is, are solutions of Eq. (2). This leads to the KL expansion stated in Eq. (1). Based on these definitions and theorems the domain independence property is stated and proved next. Proposition 1. Let κ(x, θ ) ∈ L 2 (D × Ω ) be a mean-square continuous random field defined on D ⊂ Rd , specified up to second order moments with κ(x) ¯ as the mean and C(x, y) as the covariance function. Consider another domain D′ ⊂ Rd such that D ∩ D′ is not a null set, and define another mean-square continuous random field κ ′ (x, θ ) ∈ L 2 (D′ × Ω ) with the same mean κ(x) ¯ and covariance function C(x, y). Let (λi , φi ) be the ith eigenpair of IEVP equation (2) on D and (λi′ , φi′ ) be the ith eigenpair on D′ . Let ω(x, θ ) and ω′ (x, θ ) denote the complete (that is, without truncation) KL expansions of κ(x, θ ) and κ ′ (x, θ ), respectively. Then for any x ∈ D ∩ D′ 1. E[ω(x, θ )] = E[ω′ (x, θ )] ′ 2. E[ω2 (x, θ )] = E[ω 2 (x, θ )] and for any x, y ∈ D ∩ D′ 3. E[ω(x, θ )ω(y, θ )] = E[ω′ (x, θ )ω′ (y, θ )]. Proof. Since the covariance C(x, y) between any two points x, y ∈ D ∩ D′ is the same for both the domains, by Mercer’s theorem we can conclude that C(x, y) =

∞ 

λi φi (x)φi (y) =

i=1

∞ 

λi′ φi′ (x)φi′ (y)

∀ x, y ∈ D ∩ D′ .

i=1

The KL expansions ω(x, θ ) and ω′ (x, θ ) of random fields κ(x, θ ) and κ ′ (x, θ ), respectively, are given by ω(x, θ ) = κ(x) ¯ + ω′ (x, θ ) = κ(x) ¯ +

∞  

λi φi (x)ξi (θ )

i=1 ∞   i=1

λi′ φi′ (x)ξi (θ ).

(9)

129

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Each of the conclusions will be shown separately. 1. 

∞  

E[ω(x, θ )] = E κ(x) ¯ +

 λi φi (x)ξi (θ )

i=1 ∞  

= κ(x) ¯ +

λi φi (x)E[ξi (θ )]

i=1

= κ(x) ¯

following property 1 of the Karhunen–Lo`eve theorem.

Similarly E[ω′ (x, θ )] = κ(x). ¯ Hence E[ω(x, θ )] = E[ω′ (x, θ )]. 2. 

∞  ∞    E[ω (x, θ )] = E λi λ j φi (x)φ j (x)ξi (θ )ξ j (θ )



2

i=1 j=1 ∞  ∞    λi λ j φi (x)φ j (x)E[ξi (θ )ξ j (θ )]. = i=1 j=1

Using the property 2 of ξi -s from the Karhunen–Lo`eve theorem E[ω2 (x, θ )] =

∞ 

λi φi2 (x).

i=1 ′2

Similarly E[ω (x, θ )] =

′ ′2 i=1 λi φi (x)

∞

′

and by substituting x = y in Eq. (9), E[ω2 (x, θ )] = E[ω 2 (x, θ )]

3. 

∞  ∞    λi λ j φi (x)φ j (y)ξi (θ )ξ j (θ ) E[ω(x, θ )ω(y, θ )] = E



i=1 j=1 ∞  ∞    = λi λ j φi (x)φ j (y)E[ξi (θ )ξ j (θ )]. i=1 j=1

Using the property 2 of ξi -s from the Karhunen–Lo`eve theorem E[ω(x, θ )ω(y, θ )] =

∞ 

λi φi (x)φi (y),

i=1

similarly E[ω′ (x, θ )ω′ (y, θ )] =

∞ 

λi′ φi′ (x)φi′ (y).

i=1

Using Eq. (9) E[ω(x, θ )ω(y, θ )] = E[ω′ (x, θ )ω′ (y, θ )].  Remark 1. Proposition 1 implies that the first and second order moments of a random process generated by the KL expansion are invariant to a change in the physical domain. However, no invariance is claimed for the individual KL coefficients. Remark 2. If the random field is Gaussian, the mean and second order moments completely specify the random field. Proposition 1 implies that a Gaussian random field is invariant to the change of the physical domain provided functional form of the covariance kernel is retained the same. For an example in mechanics, consider a beam where the modulus of elasticity is modeled as a Gaussian random field. Proposition 1 implies that making a hole in the beam will not alter the probability distribution of the modulus of elasticity in the remaining part of the beam.

130

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Note that the domain independence property is proved here for the complete KL expansion. However, as mentioned earlier, for computational purpose the expansion is truncated after a few terms, as expressed in Eq. (3). Now the effect of this truncation on the domain independence property is explored and an error bound is derived. Note that from Proposition 1, the KL expansions ω(x, θ ) and ω′ (x, θ ) lead to the same process. Therefore, now onwards only the notation κ(x, θ ) is used and the notation κ ′ (x, θ ) is discarded. Proposition 2. Let ω, ˆ ωˆ ′ denote the truncated versions of ω and ω′ , respectively. Let ω(x, θ ) − ω(x, ˆ θ ) = ϵ(x, θ ) ω′ (x, θ ) − ωˆ ′ (x, θ ) = ϵ ′ (x, θ ). Then ∥ω(x, ˆ θ ) − ωˆ ′ (x, θ )∥ L 2 (Ω ) ≤ (∥ϵ ′ (x, θ )∥ L 2 (Ω ) + ∥ϵ(x, θ )∥ L 2 (Ω ) )

(10)

where ∥ · ∥ L 2 (Ω ) denotes the L 2 (Ω ) norm. Proof. For brevity the explicit dependence on the arguments x and θ will not be mentioned. Now ∥ω − ω∥ ˆ L 2 (Ω ) = ∥ϵ∥ L 2 (Ω )

(11)

∥ω′ − ωˆ ′ ∥ L 2 (Ω ) = ∥ϵ ′ ∥ L 2 (Ω ) .

(12)

An error bound for ∥ωˆ − ωˆ ′ ∥ L 2 (Ω )

(13)

is now derived. Add and subtract ω to (ωˆ − ωˆ ′ ) and then use the triangle inequality in Eq. (13), to get ∥ωˆ − ωˆ ′ ∥ L 2 (Ω ) = ∥ωˆ − ω + ω − ωˆ ′ ∥ L 2 (Ω ) ≤ ∥ω − ω∥ ˆ L 2 (Ω ) + ∥ω − ωˆ ′ ∥ L 2 (Ω ) = ∥ϵ∥ L 2 (Ω ) + ∥ω − ωˆ ′ ∥ L 2 (Ω ) . Now consider ∥ω

− ω′ ∥

L 2 (Ω ) .

Add and subtract κ to ω

(14) − ω′ ,

take the norm, and use the triangle inequality to get

∥ω − ω′ ∥ L 2 (Ω ) = ∥ω − κ + κ − ω′ ∥ L 2 (Ω ) ≤ ∥κ − ω∥ L 2 (Ω ) + ∥κ − ω′ ∥ L 2 (Ω ) .

(15)

Following the convergence of the KL theorem as mentioned in Eq. (8), ∥κ − ω∥ L 2 (Ω ) = 0 and ∥κ − ω′ ∥ L 2 (Ω ) = 0 ∀x, y ∈ D ∩ D′ . Therefore using the non-negative property of norms in Eq. (15), it can be concluded that ∥ω − ω′ ∥ L 2 (Ω ) = 0. Next consider the last term on the right hand side of inequality (14) ∥ω − ωˆ ′ ∥ L 2 (Ω ) = ∥ω − ω′ + ω′ − ωˆ ′ ∥ L 2 (Ω ) ≤ ∥ω − ω′ ∥ L 2 (Ω ) + ∥ω′ − ωˆ ′ ∥ L 2 (Ω ) = ∥ω′ − ωˆ ′ ∥ L 2 (Ω )

(using triangle inequality)

(using ∥ω − ω′ ∥ L 2 (Ω ) = 0)

This inequality, combined with Eq. (12), leads to ∥ω − ωˆ ′ ∥ L 2 (Ω ) ≤ ∥ϵ ′ ∥ L 2 (Ω ) .

(16)

Substituting inequality (16) into inequality (14) we obtain ∥ωˆ − ωˆ ′ ∥ L 2 (Ω ) ≤ (∥ϵ ′ ∥ L 2 (Ω ) + ∥ϵ∥ L 2 (Ω ) ).



Note that for a desired level of accuracy the number of required terms L in the truncated KL expansions in two domains may vary. Proposition 2 implies that if the KL expansions in both the domains are truncated after sufficient number of terms, then the random processes generated by them differ only by a very small amount. Inclusion of more terms will reduce this difference further. Based on these observations the proposed method is presented next.

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

131

Fig. 1. Example 1: 2D problem, a plate with a hole. An FE mesh is also shown.

Fig. 2. Example 2: 3D problem, a structural steel angle with bolting holes. An FE mesh is also shown.

3. Proposed method using domain independence The domain independence property of the KL expansion of a random process will be used in the following way. The first step is to simplify the spatial domain. That is, the original physical domain D is replaced by a new domain D′ over which the numerical solution of the IEVP will be easier and faster. Here D′ is chosen as the minimal bounding box around D, such that D ⊂ D′ . For demonstration of the proposed idea consider two examples. First, consider a two-dimensional domain shown in Fig. 1. The mechanics problem here is a plane stress problem on a rectangular domain with a circular hole at the center—call this domain as D. For solving the KL eigenvalue problem this domain will be replaced simply by a rectangle whose sides coincide with the sides of the outer boundary of the original domain D—call this new domain as D′ . Then, consider a three-dimensional domain shown in Fig. 2, this is a structural steel connector angle with bolting holes. For this problem, the modified domain D′ is chosen as a (solid) cuboid as shown in Fig. 3, whose sides are of minimal dimensions to fit around the object. In other words, D′ is the smallest cuboid which completely contains the angle. The next step is to use the Nystr¨om method for solving the integral equation (2) over the new domain D′ . Accordingly, a numerical quadrature rule is chosen to evaluate the integral on domain D′ . Due to the rectangular (in two dimensions) or cuboid (in three dimensions) nature of the domain D′ , finding a quadrature rule over D′ will be far simpler than to find a quadrature rule over the complex shaped domain D. Furthermore, in this approach the accuracy and cost of solving the IEVP will no longer be affected by the FE mesh chosen for the mechanics problem. Using an N point quadrature rule the IEVP becomes N  j=1

( j)

( j)

w ( j) C(x1 , x2 )φi (x1 ) = λi φi (x2 )

(17)

132

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 3. Bounding box corresponding to the steel angle.

Fig. 4. Eigenvalue decay in the 2D problem.

Fig. 5. Comparison of pdfs for plate with and without hole using exponential covariance.

133

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 6. Comparison of pdfs for plate with and without hole using triangular covariance.

Fig. 7. Comparison of pdfs for the lognormal process in the 2D problem. ( j)

where {x1 } Nj=1 denotes the set of quadrature points and {w( j) } Nj=1 denotes the set of corresponding weights. Dis(k)

cretizing the above equation completely over N quadrature points by replacing x2 = x2 , k = 1, 2, . . . , N , we obtain N simultaneous equations N 

( j)

(k)

( j)

(k)

w ( j) C(x1 , x2 )φi (x1 ) = λi φi (x2 ),

k = 1, 2, . . . , N .

(18)

j=1

The above system of equations can be written in the form of a matrix eigenvalue problem as Kφ i = λi φ i

(19) ( j)

(k)

(k)

where K ∈ R N ×N with the (k, j)th element as w( j) C(x1 , x2 ), and φ i ∈ R N with the kth element as φi (x2 ). This symmetric, positive definite algebraic eigenvalue problem can be solved using a standard solver. Note that only a few largest eigenvalues and associated eigenvectors are sufficient for the KL approximation. Therefore methods such as the power iterations or Krylov-subspace based solvers might be efficient choices.

134

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 8. Decay of eigenvalues in the 3D problem.

Fig. 9. Comparison of pdfs for steel angle and the cuboid in the 3D problem.

• • • • •

The proposed method thus can be summarized as follows. Step 1: Enclose the given domain in the minimal bounding box. Step 2: Choose an appropriate quadrature rule for numerical integration over this box and calculate the corresponding weights. ( j) (k) Step 3: Construct the matrix K by evaluating w( j) C(x1 , x2 ) for k = 1, 2, . . . , N , j = 1, 2, . . . , N . Step 4: Solve the algebraic eigenvalue problem defined in Eq. (19). Step 5: Post processing: When required, interpolate the eigenfunctions as required following the quadrature rule used. For a mechanics problem the eigenfunctions need to be evaluated at the Gauss points of the FE mesh.

Because of the regular shape of the domain D′ a uniform numerical mesh can be used, which makes the implementation simple. Quadrature rules such as Simpsons rule, trapezoidal rule, and mid-point rule can be used. In the following numerical study the mid-point rule has been used. 4. Numerical studies Two stochastic mechanics problems are considered for numerical studies. The first one is a plane stress problem, where a thin plate with a circular hole at center – as shown in Fig. 1 – is subjected to an in-plane uniaxial tensile

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

135

Fig. 10. Match of pdfs using correlation length as 1/12 of the width and six KL expansion terms in the 2D problem.

Fig. 11. Match of pdfs using correlation length as 1/12 of the width and fifty KL expansion terms in the 2D problem.

force. The second one is a three-dimensional one: a connector angle with bolting holes – as shown in Fig. 2 – is subjected to a uniaxial tensile force. Young’s modulus is assumed to be a random field in these two examples, with a known covariance function. For the two-dimensional (2D) problem two probability distributions are considered for the field: Gaussian and lognormal, whereas for the three-dimensional (3D) problem only Gaussian distribution is considered. For the mechanics problem, in both cases the presence of holes leads to stress concentration around them, thus demands a refined FE mesh in these areas. For instance, see the typical meshes in Figs. 1 and 2. Whereas the random field should not require any such localized refinement. Using these two numerical examples three different numerical studies are conducted here. First, the domain independence of the KL approximation of the random field is numerically demonstrated and the effect of the correlation length on the truncation of the KL expansion is studied. Next, a computational cost comparison between the proposed method and the existing FE-discretization based Galerkin method is carried out. Both CPU time and memory requirement are compared for this purpose. Accuracy of the proposed method is also verified. Finally, for a separable kernel it is demonstrated that further cost saving is possible by a corresponding separation of the eigenvalue problem.

136

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 12. Convergence of stress with mesh refinement in the 2D problem.

Fig. 13. Convergence of variance of the simulated process with mesh refinement in the 2D problem using the Galerkin method. Table 1 Percentage relative error of eigenvalues using the proposed method: the 2D example. Mode No

Nystr¨om method

Analytical value

Percentage error

1 2 3 4

6.59 0.51 0.26 0.14

6.63 0.49 0.26 0.13

0.6 3.9 0 7.1

4.1. Numerical demonstration of domain independence For demonstrating the domain independence of the KL approximations, the widely used Galerkin approach in an FE-discretization [9,30] is used to solve the IEVP. Accordingly, a set of FE bases – usually the one used for the mechanics problem – is used to approximate the eigenfunctions. To find the unknown coefficients, this approximation is substituted in the IEVP—that is, in Eq. (2), and a Galerkin projection of the residual is taken on the same FE basis set. These steps lead to a generalized algebraic eigenvalue problem, which is then solved by a standard solver.

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

137

Fig. 14. Percentage relative error with mesh refinement for mode 1 in 2D problem using the proposed method.

Fig. 15. Convergence of variance of the simulated process with mesh refinement in the 2D problem using the proposed method.

First the 2D problem is considered. The thin plate is chosen to have an aspect ratio of two with the diameter of the hole being 0.5 times the width of the plate. For a Gaussian random field model of Young’s modulus the mean is assumed to be 200 MPa and standard deviation as 10 MPa. Two covariance kernels are considered: one exponential with the form C(x, y) = e−

|x1 −x2 | α

e−

|y1 −y2 | α

and one triangular with the form    |x1 − x2 | |y1 − y2 | C(x, y) = 1 − 1− α α

(20)

(21)

in two dimensions [30], where (x1 , y1 ) and (x2 , y2 ) denote the coordinates of two points. In both the cases the correlation length α is chosen to be five times the width of the plate. Constant strain triangle (CST) plane stress elements were used in the FE modeling. From a mesh convergence study for the mechanics problem – details of which are given in the next subsection – it is found that about 1000 elements can capture the stress field accurately. Therefore this converged mesh is used for the Galerkin method of solution of the IEVP. As mentioned earlier, the

138

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 16. Accuracy comparison between the proposed and Galerkin methods, the 2D example.

Fig. 17. Convergence of stress with mesh refinement in the 3D problem.

minimal bounding box here is a thin rectangular plate of same dimensions, but without the hole. For brevity here these two domains are referred to as with hole and without hole. The FE mesh generated for the modified domain (without hole) is independent from the mesh used in the original domain (with hole). Upon solving the eigenvalue problems for these two domains, decay of the eigenvalues of the covariance function with unit variance is plotted in Fig. 4. From this figure it is noted that both the domains have similar trend of decay, however, the numerical values of the eigenvalues vary. Based on these decay patterns six terms were retained in constructing the KL approximation. For an arbitrarily chosen spatial location the pdfs of Young’s modulus for these two domains are then estimated using the truncated KL expansion. These pdfs are plotted in Fig. 5 for the exponential covariance, and in Fig. 6 for the triangular covariance, respectively. In both of these figures an excellent match between the pdfs in two domains is observed, which was verified using a 5% Kolmogorov–Smirnov test. Often the log-normal distribution is suggested as a more realistic model for elastic properties [4]. To study the domain independence property, a lognormal process is next considered. The statistical moments for this lognormal field model of Young’s modulus are taken from [4] and an exponential covariance function is considered. Upon finding the KL approximations in two domains, the pdfs of Young’s modulus at the same spatial location are plotted in Fig. 7. Again an excellent agreement is found and this is verified using a 5% Kolmogorov–Smirnov test.

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

139

Fig. 18. Convergence of variance with mesh refinement in the 3D problem using the Galerkin method.

Fig. 19. Convergence of variance with mesh refinement in the 3D problem using the proposed method.

The 3D example is considered next. As mentioned earlier, the minimal bounding box here is a solid cuboid of same outer dimensions. Constant strain tetrahedron elements were used for the FE modeling. The modulus of elasticity of the angle was modeled as a Gaussian random field with a mean of 200 MPa and standard deviation of 10 MPa, with the covariance function being of exponential type. Decay of the eigenvalues for two domains is plotted in Fig. 8, and the pdfs of the KL approximated process evaluated at an arbitrarily chosen spatial location are plotted in Fig. 9. The observations are similar to the 2D case. That is, although the individual eigenvalues may vary with a modification in the domain, pdfs of the KL-approximated field remain domain independent. The match in the pdfs is again verified using a 5% Kolmogorov–Smirnov test. Recall that in Section 2 the importance of including sufficient number of terms was highlighted for this domain independence. This sufficiency requirement is now numerically studied for the 2D problem. In general, with a reduction in the correlation length, the eigenvalue decay becomes slower and thus the KL approximation requires 1 × width of the plate. The pdfs from two more terms to be included. Here the correlation length is changed to 12 different domains are compared with different levels of truncation. In Fig. 10 six terms were included, that is, L = 6 in Eq. (3). It is observed that the pdfs have different variance. Next the level of truncation L is increased to 50 and the pdfs are plotted in Fig. 11. Now there is an excellent agreement between the generated pdfs, which highlights the importance of including more number of terms.

140

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 20. Comparison of memory used for the 2D problem.

Therefore by numerical experiments it is found that although the individual  L √ terms λi , φi (x) in the KL expansion may vary with a change in the domain, the overall KL approximation i=1 λi φi (x)ξi (θ ) — that is, the summation over a sufficient number of KL modes L — does not depend upon the domain. 4.2. Cost and accuracy analysis of the proposed method Next the accuracy and computational cost of the proposed method is studied. In these studies the Galerkin method is applied to the original domains and the proposed method is applied to the modified domains. First the 2D example is considered. Since the accuracy is of concern, first the convergence of meshes in both the methods is studied. Convergence of the FE mesh is studied for both the mechanics problem and the random field discretization problem. Accordingly, a sequence of meshes with successive refinement is considered. The axial stress at an arbitrarily chosen point in the stress concentration region is computed for these meshes and the trend is plotted in Fig. 12. The stress has converged in a mesh with around 1000 elements. Note that this converged mesh was used in the previous section. The variance of the process estimated using the Galerkin method is plotted against the mesh size in Fig. 13. It is noted that the variance converges at around 200 elements, thus much earlier than the mechanics problem. To study the convergence of the proposed method the number of quadrature points is increased successively. Let λ N ys denote the eigenvalue found using the Nystr¨om method and λ Anal denote the eigenvalue computed analytically. For the highest |λ −λ | eigenvalue the relative error, defined as N ysλ AnalAnal ∗ 100%, is plotted against the refinement in the quadrature mesh in Fig. 14. Here the word element refers to the rectangular region within four quadrature points, this terminology is introduced only to compare with the FE meshes. For a mesh with 100 elements the relative errors in highest four eigenvalues are tabulated in Table 1. The convergence of the variance of the process with respect to the refinement in the quadrature mesh is plotted in Fig. 15. It is found that around 1400 elements yield a converged estimate. To verify the accuracy of the proposed method the pdf of Young’s modulus estimated using the proposed method (on the modified domain) and the Galerkin method (on the original domain) is plotted in Fig. 16, and an excellent match is found. This match shows that the proposed method is accurate. The same set of experiments is carried out for the 3D example. Using FE mesh the convergences of the axial stress and the variance of the process are plotted in Figs. 17 and 18, respectively. Again the variance converges with a much coarser mesh compared to the stress. This observation for both the examples suggests that the optimal meshes for the mechanics problem and the random field discretization need not be the same. Convergence of the variance using the proposed method is plotted in Fig. 19. The computational cost is studied next. For the 2D example the memory requirement and the CPU time on a single processor with 2.8 GHz speed are plotted in Figs. 20 and 21, respectively. A few important data points in these plots are tabulated in Table 2. From Fig. 20 it is noted that for any fixed mesh size the Galerkin method requires more memory compared to the proposed method. However, recall that the optimal mesh sizes for these two methods are not same. Therefore, the cost comparison should be made for the optimal meshes. To this end, consider three points in this

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

141

Fig. 21. Comparison of CPU times for the 2D problem.

Table 2 Time and memory comparison for various optimal meshes for the 2D problem.

Optimal Mechanics mesh Optimal mesh for Galerkin method Optimal mesh for proposed method

Number of elements

Time (s)

Memory (Mb)

1000 200 1400

1.2 × 104 250 0.8

29 4.8 16.2

figure marked as Q1, Q2, and Q3. They correspond to the optimal mesh for the proposed method, the optimal mesh for the mechanics problem in FE—referred to here as the optimal mechanics mesh, and the optimal mesh for the random field discretization using the FE mesh—referred to here as the optimal stochastic mesh, respectively. From these three points it is observed that the memory requirement of the proposed method is lower than the optimal mechanics mesh but higher than the optimal stochastic mesh in the Galerkin method. Next consider the CPU time comparison that is plotted in the semi-log scale in Fig. 21. Here it is noted that the proposed method is much faster compared to the Galerkin method, and the difference is of three to four orders of magnitude. The exact figures are given in Table 2. Both the optimal meshes in the Galerkin method are far more expensive than the optimal mesh in the proposed method. Note that the optimal mesh in the proposed method is finer than the stochastic optimal mesh in the Galerkin method. Thus the higher memory requirement in the proposed method can be attributed to higher mesh density. However, the cost saving can be explained using a computational complexity analysis for 2D problem. Solving Eq. (2) using the Galerkin method with FE basis results in a generalized eigenvalue problem Ax = λBx.

(22)

Here the cost of formulation is O(e2 n 2e +e2 n 2e F +en 2e G) where e is the number of elements, n is the number of nodes, n e is the number of nodes per element, F is the cost of computation for each entry in the elemental A matrix, and G is the cost of computation for each entry in the elemental B matrix. Evaluation of each entry involves numerical evaluation of the integral whose cost depends on the chosen integration rule. On the other hand, in the proposed method let the rectangular domain be divided into N × N rectangular elements. If E is the cost of computation of each entry, then the total computational cost would be O(N 2 E). It must be observed that the proposed method leads to an ordinary algebraic eigenvalue problem Eq. (19), hence requires only computation of a single matrix. Further E is much lower than F or G, as it is only the cost of evaluation of the covariance kernel. These investigations in terms of a computational complexity analysis reveal that the computation of individual elements in the matrices is far more expensive in the Galerkin method compared to the proposed method. Thus the

142

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 22. Comparison of memory used for the 3D problem.

Fig. 23. Comparison of time taken with mesh refinement in the proposed method and the Galerkin method in the 3D problem.

stage of assembling of matrices contributes dominantly to the total cost in the Galerkin method. This relative cost saving makes the proposed method computationally inexpensive. Similar observations are made for the 3D example as well. For this example the memory requirement and CPU times are plotted in Figs. 22 and 23, respectively. In Fig. 22 it is noted that the memory requirement is very close in both the methods. However, from Fig. 23 it is noted that the proposed method is about one to two orders of magnitude faster than the Galerkin method. The point Q2 is not generated here as it required far more computational time, and its contribution to the conclusions here may not be significant. 4.3. Further cost reduction for separable covariance functions Many of the commonly used covariance functions are separable [30,18]. For these functions the proposed method can be further enriched to exploit this separability, and thus to reduce the computation further. For these kernels the eigenvalues and eigenfunctions can be obtained as products of solutions of one-dimensional problems [25]. However, in this case the domain must be of a rectangular shape in two-dimensions and cuboid shape in three-dimensions. Thus, irrespective of the numerical method used to solve the integral eigenvalue problem, modification of the domain from D to D′ is necessary and the domain independence property plays a key role. In the current work, although the chosen

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

143

Fig. 24. Effect of using separability of the covariance kernel: CPU time comparison for the 2D example.

Fig. 25. Effect of using separability of the covariance kernel: CPU time comparison for the 3D example.

kernels are separable, the separability was not used in the numerical results presented in the previous subsection. This choice was made to test the proposed method for any general kernel. This separability is next used for both 2D and 3D examples, and the computational cost is presented in Figs. 24–27. All results are generated using the Nystr¨om method. In both of these plots it is noted that the separability reduces the cost by a few orders of magnitude. For both the 1400 elements mesh in 2D and the 5000 elements mesh 3D the cost saving is of four orders of magnitude. Furthermore, in both the problems the cost growth rate with respect to the mesh size is far lower while using the separability. 5. Conclusion By using the domain independence property of the KL expansion along with the Nystr¨om method, the computational cost of solving the integral eigenvalue problem is reduced by up to four orders of magnitude compared to the finite element discretization. For separable kernels the cost is reduced even further. Since the mechanics problem and the IEVP are discretized using different discretizations, optimality of the meshes in terms of cost and accuracy for these two problems can be achieved independently. The proposed method, in its current form, is simple to implement. Thus it could serve as a very efficient alternative to currently used methods, especially to the FE based

144

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

Fig. 26. Effect of using separability of the covariance kernel: Memory requirement comparison for the 2D example.

Fig. 27. Effect of using separability of the covariance kernel: Memory requirement comparison for the 3D example.

Galerkin method. It is found to be successful for a non-Gaussian processes as well. Further potential developments of this method should include exploring various quadrature rules, finding a more suitable algebraic eigensolver, and parallelizing. Acknowledgments The authors would like to acknowledge financial support from the Board of Research in Nuclear Sciences (BRNS) Grant No. 2011/36/41-BRNS/1977 and the ISRO-IISc Space Technology Cell Grant No. ISTC/CCE/DG/273. References [1] K.R. Gurley, M.A. Tognarelli, A. Kareem, Analysis and simulation tools for wind engineering, Probab. Eng. Mech. 12 (1) (1997) 9–31. [2] G.I. Schu¨eller, H.J. Pradlwarter, C.A. Schenk, Non-stationary response of large linear FE models under stochastic loading, Comput. Struct. 81 (8) (2003) 937–947. [3] G. Stefanou, M. Papadrakakis, Stochastic finite element analysis of shells with combined random material and geometric properties, Comput. Methods Appl. Mech. Engrg. 193 (1) (2004) 139–160.

S. Pranesh, D. Ghosh / Comput. Methods Appl. Mech. Engrg. 285 (2015) 125–145

145

[4] JCSS-probabilistic model code. part III resistance models — steel. Joint Committee on Structural Safety, 2001. [5] M. Ostoja-Starzewski, Random field models of heterogeneous materials, Int. J. Solids Struct. 35 (19) (1998) 2429–2455. [6] A. Nouy, A. Clement, F. Schoefs, N. Mo¨es, An extended stochastic finite element method for solving stochastic partial differential equations on random domains, Comput. Methods Appl. Mech. Engrg. 197 (51) (2008) 4663–4682. [7] C.A. Schenk, G.I. Schu¨eller, Uncertainty Assessment of Large Finite Element Systems, Springer, Berlin/Heidelberg, 2005. [8] M. Kleiber, T.D. Hien, The Stochastic Finite Element Method: Basic Perturbation Technique and Computer Implementation, Wiley, 1992. [9] R. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, revised ed., Dover Publications, 2003. [10] A. Der Kiureghian, J.-B. Ke, The stochastic finite element method in structural reliability, Probab. Eng. Mech. 3 (2) (1988) 83–91. [11] S. Mahadevan, A. Haldar, Practical random field discretization in stochastic finite element analysis, Struct. Saf. 9 (4) (1991) 283–304. [12] E. Vanmarcke, M. Grigoriu, Stochastic finite element analysis of simple beams, J. Eng. Mech. 109 (1983) 1203–1214. [13] G.A. Fenton, E.H. Vanmarcke, Simulation of random fields via local average subdivision, J. Eng. Mech. 116 (8) (1990) 1733–1749. [14] H.L. Van Trees, Detection, Estimation and Modulation Theory, Part I, Wiley, 2001. [15] C.-C. Li, A. Der Kiureghian, Optimal discretization of random fields, J. Eng. Mech. 119 (6) (1993) 1136–1154. [16] G. Deodatis, M. Shinozuka, The weighted integral method. II: response variability and reliability, J. Eng. Mech. 117 (8) (1991) 1865–1877. [17] B.A. Zeldin, P.D. Spanos, On random field discretization in stochastic finite elements, Trans. ASME J. Appl. Mech. 65 (1998) 320–327. [18] P.D. Spanos, M. Beer, J. Red-Horse, Karhunen–Lo`eve expansion of stochastic processes with a modified exponential covariance kernel, J. Eng. Mech. 133 (7) (2007) 773–779. [19] R. Ghanem, The nonlinear Gaussian spectrum of log-normal stochastic processes and variables, J. Appl. Mech. 66 (4) (1999) 964–973. [20] K.K. Phoon, S.P. Huang, S.T. Quek, Simulation of second-order processes using Karhunen–Lo`eve expansion, Comput. Struct. 80 (12) (2002) 1049–1060. [21] K.K. Phoon, H.W. Huang, S.T. Quek, Simulation of strongly non-gaussian processes using Karhunen–Lo`eve expansion, Probab. Eng. Mech. 20 (2) (2005) 188–198. [22] K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, 1997. [23] G. Stefanou, The stochastic finite element method: past, present and future, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1031–1051. [24] D.L. Allaix, V.I. Carbone, Numerical discretization of stationary random processes, Probab. Eng. Mech. 25 (2010) 332–347. [25] C. Schwab, R.A. Todor, Karhunen–Lo`eve approximation of random fields by generalized fast multipole methods, J. Comput. Phys. 217 (1) (2006) 100–122. [26] B. Khoromskij, A. Litvinenko, H. Matthies, Application of hierarchical matrices for computing the Karhunen–Lo`eve expansion, Computing 84 (2009) 49–67. [27] C.F. Li, Y.T. Feng, D.R.J. Owen, D.F. Li, I.M. Davis, A Fourier–Karhunen-0Lo`eve discretization scheme for stationary random material properties in SFEM, Internat. J. Numer. Methods Engrg. 73 (2008) 1942–1965. [28] D.C. Charmpis, M. Papadrakakis, Improving the computational efficiency in finite element analysis of shells with uncertain properties, Comput. Methods Appl. Mech. Engrg. 194 (2005) 1447–1478. [29] R. Guti´errez, J.C. Ruiz, M.J. Valderrama, On the numerical expansion of a second order stochastic process, Appl. Stoch. Models Data Anal. 8 (2) (1992) 67–77. [30] S.P. Huang, S.T. Quek, K.K. Phoon, Convergence study of the truncated Karhunen–Lo`eve expansion for simulation of stochastic processes, Internat. J. Numer. Methods Engrg. 52 (9) (2001) 1029–1043. [31] S. Choudhary, Numerical methods for solving the eigenvalue problem involved in the Karhunen-Lo`eve decomposition (Ph.D. thesis), Indian Institute of Science, 2012. [32] W. Betz, I. Papaioannou, D. Straub, Numerical methods for the discretization of random fields by means of the Karhunen–Lo`eve expansion, Comput. Methods Appl. Mech. Engrg. 271 (2014) 109–129. [33] J. Parvizian, A. D¨uster, E. Rank, Finite cell method, Comput. Mech. 41 (1) (2007) 121–133.