A Modified Computational Scheme for the Stochastic Perturbation Finite Element Method

Abstract

A modified computational scheme of the stochastic perturbation finite element method (SPFEM) is developed for structures with low-level uncertainties. The proposed scheme can provide second-order estimates of the mean and variance without differentiating the system matrices with respect to the random variables. When the proposed scheme is used, it involves finite analyses of deterministic systems. In the case of one random variable with a symmetric probability density function, the proposed computational scheme can even provide a result with fifth-order accuracy. Compared with the traditional computational scheme of SPFEM, the proposed scheme is more convenient for numerical implementation. Four numerical examples demonstrate that the proposed scheme can be used in linear or nonlinear structures with correlated or uncorrelated random variables.

Keywords:
Stochastic finite element method; Uncertain structures; Perturbation; Computational scheme

1 INTRODUCTION

The responses of actual structures depend on many uncertain parameters, such as unpredictable external excitations, random material parameters, random geometric properties, etc. Analyses of the influences of these uncertain parameters on actual structures have become an important topic in the field of computational mechanics and have attracted the attention of many researchers.

Many methods have currently been developed for analyzing structures with uncertain parameters; (Stefanou 2009). (Falsone and Ferro 2007) developed a method that can provide explicit relationships between the system responses and the random variables defining the uncertain parameters. (Kamiński 2010) proposed a perturbation-based stochastic finite element method using the polynomial response function. (Kim and Inoue 2004) suggested a spectral stochastic element-free Galerkin method for problems involving the random material property. (Li and Chen 2004) presented a method based on the generalized probability density evolution. (Lin et al. 2001) proposed the pseudo-excitation method for the frequency domain analysis of structures with random excitations. (Deb et al. 2001) developed a variational formulation for the elliptic boundary-value problems with stochastic input data. Their method allows numerical treatment by the finite element method. (Ganapathysubramanian and Zabaras 2007) proposed the sparse grid method (SGM) for stochastic problems. (Rahman 2008) presented a new polynomial dimensional decomposition method for solving stochastic problems.

Among the proposed stochastic approaches, the Monte Carlo simulation (MCS) may be the most widely used method; (Papadrakakis and Kotsopulos 1999), (Szekely and Schuëller 2001), (Johnson et al. 2003), (Lei and Qiu 2000), (Boyaval 2012). MCS is a type of statistical approach based on the law of large numbers. It is applicable to any stochastic system. However, this method is often time-consuming because the errors of the MCS results depend on the size of the sample set. The quais-Monte Carlo method(QMC) is based on the discrepancy of point sets and is often more efficient than the MCS; (Holtz 2011). The spectral stochastic finite element method (SSFEM) is another type of stochastic method, in which the random field is usually discretized using the Karhunen-Loève expansion and the structure's nodal displacements are approximated using a polynomial chaos expansion; (Sakata et al. 2011), (Xiu and Karniadakis 2002) The computational cost of the SSFEM is smaller than that of the MCS. However, the computational effort increases exponentially with the order and number of uncertain quantities involved, which places some practical restrictions on this method. For a system with low-level uncertainties, the stochastic perturbation finite element method (SPFEM) was proposed, which can provide the second-order estimates of the mean and covariance; (Liu et al. 1986), (Van den Nieuwenhof and Coyette 2003), (Kamiński 2013). If the coefficient of variation is too large (e.g., larger than 20 percent), the results of the SPFEM are usually unacceptable. To overcome this defect, the n-th order stochastic perturbation technique was developed with the help of a symbolic computer program to extend the applied range of the SPFEM; (Kamiński 2007), (Kamiński 2006) For the case of a large coefficient of variation, a sub-interval perturbation technique was proposed; (Xia and Yu 2012). In order to further reduce the computational effort of the SPFEM, a transformation technique was combined with the SPFEM to analyze the structure with correlated random variables ; (Liu et al. 1986).

Although there are many limits to the SPFEM, it is still widely used for two reasons; (Kamiński 2010), (Kamiński 2013), (Kamiński 2011), (Wang and Yao 2010), (Sakata et al. 2011), (Kamiński and Szafran 2009), (Chang et al. 2008). The first reason is that the SPFEM is often more efficient than other stochastic approaches, and the second one is that for a real engineering structure with random material and geometric parameters, the coefficients of variation of these random parameters are usually small. However, compared with the MCS, the traditional computational scheme of the SPFEM is too complicated because it requires computing the derivatives of the system matrices (e.g., the stiffness and mass matrices) with respect to the random variables. Generally, the traditional computational scheme depends on the type of problem that one takes into account, and each type of problem has its own computational scheme.

In this paper, a unified computational scheme of the SPFEM, which can provide the second-order estimates of the mean and covariance matrices of the system responses, is developed. In the proposed computational scheme, the derivatives of the system matrices with respect to the random variables are not required. Compared with the original computational scheme, the modified scheme is more convenient for complicated problems, and the computational procedures are the same for different types of random problems. In the case of one random variable with a symmetric probability density function (PDF), the proposed computational scheme can even provide the result with fifth-order accuracy.

In the next section, the basic theory of the SPFEM is outlined. In Section 3, the proposed computational scheme of the SPFEM is addressed in detail. In Section 4, four numerical examples are presented. The numerical results are compared with the results obtained using the MCS and the original computational scheme of the SPFEM.

2 THE STOCHASTIC PERTURBATION FINITE ELEMENT METHOD

In this section, the basic theory of the SPFEM is briefly introduced. More information about the SPFEM can be found in works of (Liu et al. 1986), (Wu and Zhong 2013). Suppose that the model of the random static system produced by the finite element method can be written as

where ε = (ε ₁, ε ₂,... ε i,..., ε q) is a zero-mean random vector, in which εi is the random variable and q is the number of the random variables.K(ε) is the stochastic stiffness matrix, and u(ε) and f are the displacement and load vectors, respectively.

Using the SPFEM for Eq. (1) involves inevitably the partial derivations of random quantities with respect to random variables. For convenience, some notations are introduced first. If the random quantity x(E) is function of the zero-mean random vector ε, i.e., E(ε) = 0, then

where x₀ denotes x(ε) evaluated at ε = 0, x_1,i denotes the partial derivate of x(ε) with respect to εi evaluated at ε = 0, and x_2,ij denotes the mixed partial derivate of x(ε) with respect to ε i and ε j (in the case of i ≠ j) or half of the second-order partial derivate of x(ε) with respect to ε i (in the case of i = j) evaluated at ε = 0.

The Taylor series expansion of the stochastic stiffness matrix K(E) can be expressed as

where K₀, K 1,i and K ₂,_ij are defined by Eq. (2), and O() represents the truncated remainder which satisfies the following condition:

The displacement vector u can also be expressed as

where u₀, u1,i and u_2,ij are defined by Eq. (2). Substituting Eqs. (3)- (5) into Eq. (1) and equating equal terms gives

and

In terms of Eq. (5), we have

Using Eqs. (5) and (9), the mean vector and covariance matrices of the displacement vector can then be written as:

And

Where

cov(ε_i, ε_j) = σ_iσ_jρ_ij, σ = (σ₁, σ₂, ..., σ_i, ..., σ_q)

σi is the standard deviation of εi, and o is a vector consisting of σi. ||σ|| σ is the infinite norm of the vector σ. ρij is the correlation coefficient of the random variables εi and εj. Eq. (11) is the appropriate truncation of cov (u (ε), u (ε)) in the sense of second-order accuracy and has been introduced in many reports; (Stefanou 2009), (Liu et al. 1986), (Falsone and Impollonia 2002). If the random variables are uncorrelated, i.e., E(εi, εj) = 0 Eq. (10) and Eq. (11) can be simplified as

And

respectively. Equations (6)-(8) show that when the random variables are correlated, the computations of the u₀, u1,i and u2ij vectors involve img deterministic equations. However, Eqs. (12) and (13) show that when the random variables are uncorrelated, only u₀, u1,i and u2ij are needed. Therefore, only deterministic equations need to be computed in this case. For the case of correlated random variables, (Liu et al. 1986) proposed a technique that can transform correlated random variables into uncorrelated random variables.

Although the SPFEM is widely applied to systems with low-level uncertainties, the computational scheme of this method is too complex because it requires the derivatives of the system matrices with respect to random variables, e.g., K1,i and K2ij, which may be difficult to obtain for some complex systems.

3 THE MODIFIED COMPUTATIONAL SCHEME OF THE SPFEM

In the original computational scheme of the SPFEM, which involves computing the derivate matrices (i.e., K1,i and K2ij) of the system matrix K(ε) with respect to the random vector ε evaluated at E(ε), Eqs. (5)-(7) are used to compute the u₀, u1,i and u2ij vectors. However, the derivate matrices are difficult to compute in some cases. To avoid differentiating the system matrices, a modified computational scheme is developed in this section. The modified scheme uses another technique, which is based on the second-order estimations of the mean and covariance matrices derived from the SPFEM (i.e., Eqs. (10)- (11)), rather than Eqs. (6)-(8), to calculate the u1,iσi and u2ijσiσj vectors directly. The proposed computational scheme is developed with the help of the static FEM model, but it can be extended to other problems without any difficulty.

Because the case of uncorrelated random variables is often encountered, it is discussed first in Subsection 3.1. The special case in which the random variables are uncorrelated and have a symmetric joint PDF is addressed in Subsection 3.2. The case of correlated random variables is addressed in Subsection 3.3. The computational scheme presented in the previous three subsections is developed in terms of the static FEM model and extended to other problems in Subsection 3.4.

3.1. The case of Uncorrelated Random Variables

When the SPFEM is used, the random variables ε i are treated naturally as small quantities (i.e., |ε_i| << 1), and the displacement vector is expressed as a Taylor series, i.e., Eq. (5), which contains the terms u1,i ε i and u₂,ij ε i ε j. However, it is worth noting in terms of Eqs. (10)-(13) that u1,iσi and u2,ijσiσj are needed rather than u1,i ε i and u2,ij ε i ε j, and ε i is an uncertain variable. Hence, we can focus our attention on computing u1,iσi and u2,ijσiσj.

For the case of uncorrelated random variables, the third-order estimates of the mean and covariance matrices of the displacement vector can be written as

where

If we replace the random vector ε in Eq. (5) with a deterministic vector a_s which is defined as

then we have

Replacing the random vector ε in Eq. (5) with -a_s yields

Using Eqs. (17) and (18), it is easy to obtain the following two equations:

And

Where

Note that

The combination of Eqs. (19), (20) and (22) yields

And

Substituting Eqs. (19), (20), (23) and (24) into Eq. (14), the mean and covariance matrices of the displacement vector can be written as

And

Noting that z_s = O(), Eq. (26) can be rewritten as

Equations (25) and (27) show that the proposed estimates of the mean and covariance matrices of the displacement vector need only the vectors u₀, w_s and z_s, in which u₀ can be obtained using Eq. (6), and in terms of Eq. (21), the computations of the vectors w_s and z_s involve only the displacement vectors u(a_s) and u(-a_s). According to the form of a_s, i.e., Eq. (16), u(a_s) and u(-a_s) can be computed by the following equation:

The computations of the vectors u₀, u(a_s) and u(-a_s) require solving 2q+1 deterministic equations. The derivative matrices of the system matrix K(ε) with respect to the random vector ε, which are difficult to compute in some problems, do not need to be computed in the proposed method. Therefore, compared with the original scheme of SPFEM, the proposed scheme is more convenient and the computational efforts are the same.

It is worth noting that the proposed estimate of the covariance matrix given by Eq. (27) contains the third-order terms u_1sρsss and u_2,ssρsss. These third-order terms are not contained in the estimate given by the original computational scheme, i.e., Eq. (13). Hence, the proposed estimate is more accurate than the original one.

Moreover, if the random variables are not only uncorrelated but also independent, i.e., ρijk ≠ 0 holds if and only if i = j = k. In this case, Eq. (26) can be rewritten as

Therefore, if the random variables are not only uncorrelated but also independent, Eq. (29) gives a third-order estimate of the covariance matrix that is more accurate than the estimate given by the original computational scheme of the SPFEM.

When the mean and covariance matrices of the random variables are the only known information, the terms multiplied by ρsss in Eq. (27) can be omitted, and the covariance matrix of the displacement vector can be written as

It can be shown that all second-order terms are retained in Eq. (30), so the accuracy of Eq. (30) is the same as that of the estimate given by the original computational scheme of Eq. (13).

Sometimes, it is necessary to evaluate the statistical moments of some other response vectors, e.g., the stress or strain vector, which can be obtained by using the same approach proposed for the displacement vector. Generally, the other response vector is a function of the displacement u and can be denoted by

Therefore, F is also a function of ε and can be expanded via the Taylor series as

where F₀, F_1,i and F_2,ij are defined by Eq. (2).

The mean and covariance matrices of F(ε) can then be written as

And

Using the same procedure for the estimate of the displacement given in this subsection, we have

where

3.2. The Case that the Random Variables are Uncorrelated and Have a Symmetric Joint PDF

When the random variables are uncorrelated and have a symmetric joint PDF, more higher-order information can be retained in the estimates of the mean and covariance matrices. Using the properties of the symmetric joint PDF, we have

And

The fourth-order Taylor series expansion of the displacement vector can be written as

Combing Eqs. (37) and (38) with Eq. (39), the fourth-order estimates of the mean and covariance matrices of the displacement vector can be written as

Where

If we replace the random vector ε in Eq. (39) with a deterministic vector b_s, which is defined as

we have

Replacing ε in Eq. (39) with a deterministic vector -b_s gives

Adding Eq. (43) and Eq. (44) and subtracting Eq. (44) from Eq. (43), respectively, will yield the following equations

and

where

The following equations are then obtained in terms of Eqs. (45), (21) and (46)

Substituting Eqs. (45), (46) and (48) into Eq. (40), the mean and covariance matrices can be written as

and

Equations (49) and (50) show that the proposed estimates of the mean and covariance matrices of the displacement vector need only the vectors u₀, ws and zs, in which u₀ can be obtained by solving Eq. (6). In terms of Eq.(47), the computations of the vectors ws and zs involve only the displacement vectors u(± bs) which can be computed by the following equation

The computations of the vectors u₀, u(bs) and u(-bs) require solving 2q+1 deterministic equations. Hence, the computational effort of the proposed computational scheme is the same as that of the original scheme.

Equation (40) shows that the fourth-order estimates of the mean and covariance matrices of the displacement vector include the coupling terms and the uncoupling terms , , and . Comparing Eqs. (45)- (48) with Eqs. (49) and (50) shows that the estimates of the mean and covariance matrices of the displacement vector include all the fourth-order uncoupling terms. However, these terms are not included in the estimates given by the original computational scheme of the SPFEM. Hence, the modified computational scheme is more accurate than the original one.

As shown above, for the case of multi random variables, Eqs. (49) and (50) give third-order estimates of the mean and covariance matrices. However, if there is only one random variable, Eqs. (49) and (50) can provide more accurate estimates of the mean and covariance matrices. In this case, we have

Thus, Eq. (40) can be rewritten as

Substituting Eqs. (45), (46) and (48) into Eq. (53) yields

Equation (54) shows that if there is only one random variable, the fifth-order accuracy scheme can be given for the estimates of the mean and covariance matrices.

3.3. The Case of Ccorrelated Random Variables

In Subsections 3.1 and 3.2, the case of uncorrelated random variables is addressed in detail. When the components in the random vector e are correlated, the covariance matrix cov(ε, ε) is usually a symmetric positive definite matrix rather than a diagonal matrix. (Liu et al. 1986) once proposed a transformation technique to simplify the computational scheme of the SPFEM. By using the transformation technique, the correlated variables can be transformed into a set of uncorrelated variables. This transformation technique can also be combined with the methods proposed in the previous two subsections to address the problems involving correlated random variables.

The transformation technique involves the eigenvalues of cov(ε, ε), defined by the following equations

where Φ is the eigenvector matrix of cov(ε, ε), Λ is a diagonal matrix whose diagonal elements are the eigenvalues of cov(ε, ε), and I is a unit matrix. Let

substituting Eq. (56) into Eq. (55) gives

Equation (57) shows that the covariance matrix of b is diagonal, which implies that the random vector b can be viewed as an uncorrelated random vector.

Once the random vector b is obtained, the computational schemes presented in the previous two subsections can be applied. Generally, the computational scheme developed in Subsection 3.1 could be used in the case of random variables with a symmetric or asymmetric joint PDF. Let the joint PDF of ε and b be defined by p(ε) and g(b), respectively. If the joint PDF p(ε) is symmetric, the joint PDF g(b)= p(ΦTb) of the random vector b satisfies g(b) = g(-b). Therefore, the joint PDF of the random vector b is symmetric, and the computational scheme addressed in Subsection 3.2 can be used.

The expressions of ρsss and ρssss for the random vector b can be written as

where Λs is the s-th eigenvalue of the cov(ε, ε) matrix. One can use numerical integration approaches to calculate Eq. (58), e.g., the MCS. However, if the joint PDF of e is a Gaussian function, ρssss =3; if the joint PDF of ε is a uniform PDF, ρssss =9/5. By transforming the random vector ε into b, only 2q + 1 deterministic equations need to be solved. Furthermore, the variables in the vector b have small variance values can be omitted to further reduce the computational effort.

3.4. Analyses of Other Problems

In Subsection 3.1, the static linear system was discussed. However, the computational scheme presented there can be applied to other problems, such as the nonlinear problem, the eigenvalue problem, and so on. Suppose that the problem can be expressed as

where the operator L depends on the problem considered. For example, for the statics linear system discussed in Subsection 3.1, L is given by

It is worth noting that the proposed computational schemes, i.e., Eqs. (25) and (27) or Eqs. (49) and (50), are developed based on Eqs. (10)-(13), which are the mean and covariance matrices of the displacement vector obtained using the SPFEM, rather than the FEM model (i.e., Eq. (1)) of a random static system. In other words, Eqs. (25) and (27) or Eqs. (49) and (50) are independent of the problem. If the method proposed in this paper is employed to solve Eq. (59), the solution can be expanded as

This expansion is the same as Eq. (5), and the mean and covariance are also the same as Eqs. (12) and (13). Hence, the computational schemes presented in subsections 3.1 and 3.2 can also be used for the problem governed by Eq. (59). If the random variables are uncorrelated, Eqs. (25) and (27) can be used to compute the mean and covariance matrices of the u vector; if the uncorrelated random variables have a symmetric joint PDF, Eqs. (49) and (50) can be employed. In the former case, the term u(± as) = u(0,..., ±σs,...,0) can be obtained by solving the follow equation

In the latter case, the term u(±bs) = u(0,..., ±σs,...,0) is obtained by solving the follow equation

In summary, the modified computational scheme of SPFEM proposed in Section 3 can be applied to the problems that can be solved using the original scheme, whereas the modified scheme does not require the derivatives of the system matrices with respect to the random variables. In the original scheme of the SPFEM, the expressions of the mean and covariance matrices include the terms of u1,i and u2,ij, which can be obtained using the perturbation technique. For a stochastic problem, the original computational scheme of the SPFEM depends on the form of the perturbation solution of the problem considered. It is well known that using the perturbation technique for different types of problems, the forms (i.e., the terms of u1,i and u2,ij) of the solutions are always different. Hence, for the original scheme of the SPFEM, each type of problem has its own computational scheme that depends on the problem which one takes into account. For example, Eqs. (6)-(8) can be used for static problems, but they cannot be used to other problems, such as the eigenvalue problem. However, the proposed computational scheme does not require computing the terms of u1,i and u2,ij, which means that an unified computational procedure is developed for random systems.

4 NUMERICAL EXAMPLES

4.1 The One-dimensional Elastic Bar

Consider the simple structural mechanics problem which has been discussed by (Kleiber and Hien 1992). Fig. 1 shows a cantilever bar of length l, cross-sectional area A and Young's modulus E subjected to axial load Q.

Figure 1
Cantilever bar subjected to axial force.

The Young's modulus is selected as the uncertain parameter, which can be defined by

The problem is to evaluate the first two statistical moments for the displacement at the free end of the bar. The mean and variance of the displacement q(ε) obtained using the original stochastic perturbation method have been found in the work of (Kleiber and Hien 1992) and may be shown as

where the subscript denotes the stochastic perturbation method. Assuming that the random variable ε ∈ [-a,a] is a random variable with uniform probability density function, then

Substituting Eq. (66) into Eq. (65) gives

The exact mean and variance of the displacement q(ε) may be written as

where the subscript denotes the exact solutions. Using the proposed method, the mean and variance of the displacement can be written as

where the subscript denotes the modified stochastic perturbation method.

The percentage errors of the results obtained using the original stochastic perturbation method and the proposed one are defined by

Figures 2(a) and 2(b) show the percentage errors of the means and variances, respectively. It can be shown that compared with the original stochastic perturbation method, the proposed modified one produces more precise evaluations for both means and variances. While a > 0.3, e_c,spm > 10%. Hence the variances computed using the original stochastic perturbation method are unacceptable in the case of a > 0.3. However the modified stochastic perturbation method performs well even in the case of a = 0.5.

Figure 2
Percentage errors for means (a) and variance (b) of the displacement.

4.2 The Plane Steel Frame

The eigenvalue problem of a plane steel frame is taken into account. The plane frame of a 20 story building is depicted in Fig. 3. All the stories of the building have equal masses, m = 168750kg. The uncertain story stiffness for each story is considered, which is defined by

where εi are assumed to be the independent zero-mean random variables with the same truncated Gaussian distribution functions, defined by

Figure 3
Plane frame.

The MCS, the sparse grid method (SGM) and the original and modified computational schemes of the SPFEM are employed in this example for σ = 0.1, 0.2 and 0.25. For convience, the original and modified computational schemes of the SPFEM are denoted as OSPFEM and MSPFEM, respectively. The solutions obtained using the MCS with 10,000 samples are considered as the reference solutions. The percentage errors e_max of the frequencies are defined by

in which y_i,_MCS is the mean or standard deviation of the i^th frequency computed using the MCS and y_i is the mean or standard deviation of the i^th frequency computed using the OSPFEM or MSPFEM. Compared with MCS, the proposed method involves only 41 calculations of deterministic eigenvalue problems. For the SGM, the Gauss-Hermite quadrature formula with two integration points is used. Hence, the SGM also requires 41 calculations of deterministic eigenvalue problems.

Tables 1 and 2 show the means and standard deviations of the first five frequencies, respectively. The tables show that for the eigenvalue problem considered here, the solutions computed using the proposed scheme are in excellent agreement with those of the MCS and that compared with the original computational scheme of the SPFEM, the proposed one is more accurate. In the case of σ = 0.2, the results of the original scheme are unacceptable. However, the proposed scheme still performs well. In the case of σ = 0.2, the means computed using the SGM are quite accurate, however the standard deviations are bad. For the SGM, it requires more calculations of deterministic eigenvalue problems, which lead the increase of the computational effort. Hence, the proposed scheme performs better than the SGM or the original computational scheme of the SPFEM, under the condition that the numbers of involved deterministic systems are the same.

Table 1
The means (rad/s) of the first five frequencies.

	σ = 0.1		σ = 0.2
Frequencies	MCS	OSPFEM	MSPFEM	SGM	MCS	OSPFEM	MSPFEM	SGM
1	3.061	3.062	3.061	3.062	3.013	3.019	3.012	3.017
2	9.166	9.167	9.166	9.167	9.018	9.039	9.019	9.033
3	15.217	15.219	15.217	15.219	14.977	15.006	14.974	14.996
4	21.179	21.182	21.179	21.181	20.847	20.886	20.844	20.873
5	27.019	27.020	27.017	27.019	26.596	26.644	26.594	26.628
emax	-	0.019%	0.006%	0.015%	-	0.231%	0.025%	0.162%

Table 2
The standard deviations (rad/s) of the first five frequencies.

	σ = 0.1		σ = 0.2
Frequencies	MCS	OSPFEM	MSPFEM	SGM	MCS	OSPFEM	MSPFEM	SGM
1	0.042	0.041	0.042	0.039	0.092	0.082	0.095	0.063
2	0.125	0.123	0.127	0.117	0.273	0.245	0.284	0.188
3	0.209	0.204	0.211	0.194	0.456	0.407	0.470	0.313
4	0.290	0.284	0.293	0.269	0.627	0.567	0.652	0.436
5	0.367	0.362	0.374	0.344	0.804	0.723	0.828	0.557
emax	-	2.280%	1.994%	7.138%	-	10.691%	3.889%	31.379%

Table 3 shows the means and standard deviations of the first five frequencies in the case of σ = 0.25. It can be seen from Table 3 that the means computed by the use of the original scheme, the proposed scheme and the SGM agree well with that computed using the MCS. However, the standard deviations computed using these methods are unacceptable, although the proposed scheme performs a little better. Hence, for the case of σ ≥ 0.2 the higher-order perturbation technique is required. However, for the case of low-level uncertainties, the classical SPFEM and the proposed scheme can give the acceptable estimates of the means and standard deviations with the small computational amount.

Table 3
The means (rad/s) and standard deviations (rad/s) of the first five frequencies for σ = 0.25.

	Means		Standard deviations
Frequencies	MCS	OSPFEM	MSPFEM	SGM	MCS	OSPFEM	MSPFEM	SGM
1	2.973	2.987	2.969	2.981	0.123	0.102	0.129	0.059
2	8.900	8.943	8.890	8.927	0.370	0.307	0.386	0.176
3	14.777	14.847	14.763	14.822	0.614	0.509	0.638	0.294
4	20.582	20.664	20.555	20.631	0.845	0.709	0.883	0.411
5	26.261	26.362	26.235	26.323	1.072	0.904	1.117	0.528
emax	-	0.483%	0.130%	0.309%	-	17.158%	4.737%	52.393%

The CPU times for all methods are compared in Table 4 which shows the computational cost of the MSPFEM is almost as much as the OSPFEM or the SGM. Combining Tables 1-3 with Table 4 shows that the proposed MSPFEM cost nearly as much CPU time as the OSPFEM to produce more accurate evalutations.

Table 4
CPU times (s) for example 4.2.

σ	MCS	OSPFEM	MSPFEM	SGM
0.1	3.438	0.075	0.005	0.006
0.2	2.924	0.079	0.004	0.005
0.25	3.117	0.958	0.006	0.004

Note that for the eigenvalue problem considered in this example, the form of the perturbation solution of the eigenvalue is quite complicated. To demonstrate this point, a brief introduction about the perturbation estimate of the eigenvalue is given. More information can be found in the work of (Wu and Zhong 2013). If the eigenproblem considered here is defined by

And

in which ε is the zero-mean random vector, xk(ε) and λk(ε) are the kth eigenvector and eigenvalue, respectively. Using the Taylor series expansion for the stochastic stiffness matrix k(ε), the kth stochastic eigenvector xk(ε) and eigenvalue λk(ε) gives

and

where K₀, K 1,i, K 2,ij, xk,0, xk,1,i, xk,2,ij, λk,0, λk,1,i and λk,2,ij are defined by Eq. (2). Substituting Eqs. (76)-(78) into Eq. (74) gives

and

where

and

In terms of Eqs. (79)-(84), λk,0, λk,1,i and λk,2,ii can be computed. Once λk,0, λk,1,i and λk,2,ii are obtained, the mean and covariance of the eigenvalue can be computed using the following equation

Comparing Eqs. (79)- (84) with Eqs. (6)-(8) can show that the original computational scheme of the SPFEM for the stochastic eigenvalue problems is far different from that for the stochastic linear static problem. However, it is very convenient to use the proposed scheme for this problem, i.e., we just need to compute the following deterministic eigenvalue problems:

And

in which K₀ denotes the random stiffness matrix K(ε) evaluated at the mean value ε = 0, and K(± bs) denotes the random stiffness matrix K(ε) evaluated at ε = (0, ..., ±σs, ..., 0).

4.3 The Nonlinear Truss Structure

The nonlinear truss structure represented in Fig. 4 is considered. The involved parameters are L = 10m and f = 1000N, and the sectional area of each bar is A = 0.04m². Figure 5 shows the constitutive relationship of each bar, in which E⁽⁺⁾ = 210x10⁹N/m² and E = 2E⁽⁺⁾(1+ ε i), i = 1,...,5, is considered to be the random parameter. The random variables ε i ∈ [-a,a] E are independent random variables with uniform probability density functions. For this example, the constitutive law is composed of two straight lines with different slopes, hence we need to judge the stress state of each bar which may be tensed or be compressed. Actually the governing equations of the nonlinear truss cannot be solved using directly the Newton-iterative method. The parametric quadratic programming method is used for the governing equations; (Zhang et al. 2002).

Figure 4
The truss.

Figure 5
The constitutive relationship.

For this nonlinear truss, it is difficult for the perturbation technique to develop the computational scheme of the original SPFEM, but it is convenient to address this problem using the proposed scheme. To test the accuracy of the proposed scheme, the MCS with 10,000 samples is performed, and these solutions are taken as the reference solutions. Compared with MCS, the proposed scheme involves only 11 calculations of deterministic nonlinear equations. The SGM is also employed to this problem. The Gauss-Hermite quadrature with two integration points are used to the SGM, and hence, the SGM also requires 11 calculations of deterministic nonlinear equations.

Tables 5~7 show the means and standard deviations of the horizontal displacements at nodes A and B computed using the MCS, the proposed approach and the SGM for the case of a = 0.2, a = 0.4.and a = 0.6, respectively. The three tables show that for this nonlinear truss, the results of the proposed scheme agree quite well with the reference solutions. Even in the case of a = 0.6, the proposed scheme still performs well. It can seen from Tables 5-7 that the means computed using the proposed scheme agree well with that computed using the SGM under the condition that the numbers of involved deterministic systems for both methods are the same. However, the standard deviations computed using the proposed scheme are better than that computed using the SGM in terms of their percentage errors.

Table 5
The means (m) and standard deviations (m) of the horizontal displacements at node A and B for a = 0.2.

	Means		Standard deviations
	MCS	MSPFEM	SGM	MCS	MSPFEM	SGM
Node A	1.7724E-6	1.7714E-6	1.77E-06	7.6078E-8	7.5895E-8	7.44E-08
Percentage errors	-	0.0551%	0.06%	-	0.2410%	2.17%
Node B	1.5385E-6	1.5380E-6	1.54E06	6.7859E-8	6.7519E-8	6.67E-08
Percentage errors	-	0.0373%	0.06%	-	0.5005%	1.76%

Table 6
The means (m) and standard deviations (m) of the horizontal displacements at node A and B for a = 0.4.

	Means		Standard deviations
	MCS	MSPFEM	SGM	MCS	MSPFEM	SGM
Node A	1.8007E-6	1.7998E-6	1.80E-06	1.6195E-7	1.6038E-7	1.48E-07
Percentage errors	-	-	0.0514%	0.11%	0.9705%	8.62%
Node B	1.5561E-6	1.5557E-6	1.56E-06	1.4364E-7	1.4282E-7	1.36E-07
Percentage errors	-	0.0297%	0.07%	-	0.5708%	5.58%

Table 7
The means (m) and standard deviations (m) of the horizontal displacements at node A and B for a = 0.6.

	Means		Standard deviations
	MCS	MSPFEM	SGM	MCS	MSPFEM	SGM
Node A	1.85E-06	1.85E-06	1.85E-06	2.59E-07	2.65E-07	2.19E-07
Percentage errors	-	0.24%	-0.10%	-	-2.19%	15.55%
Node B	1.59E-06	1.59E-06	1.59E-06	2.29E-07	2.36E-07	2.10E-07
Percentage errors	-	0.34%	0.13%	-	-3.35%	8.35%

Table 8 shows the CPU times for all methods. Combining Tables 5-7 with Table 8 shows that the porposed MSPFEM cost nearly as much CPU time as the SGM to produce more accurate evalutations.

Table 8
CPU times (s) for example 4.3.

A	MCS	MSPFEM	SGM
0.2	8.900	0.007	0.008
0.4	9.009	0.007	0.007
0.6	8.978	0.006	0.007

4.4 The Heat Conduction Problem

Consider the non-dimensional transient heat conduction in a square domain. The governing equation is

where c = 1 is the heat capacity, ρ = 1 is the mass density, T is the temperature, k(x,y) is the thermal conductivity, and W(x,y,t) is the volumetric heat supply. The boundary conditions and initial conditions are

and

respectively. The heat supply is defined by

The thermal conductivity is selected as the uncertain parameter which can be defined by

where ε(x,y) is a zero-mean homogeneous stochastic filed with a constant standard deviation σ = 1.5. The correlation function is defined by

in which λ = 3L is the correlation length. The midpoint method is adopted for the random field, and the 40x40 mesh of four-node linear rectangular elements is used for the space discretization. The heat conduction equations are integrated by the Crank-Nicholson method with the time step Δt = 0.01. The integral interval is t ∈ [0,2]. The number of random variables is 1600. The correlated random variables are transformed into a new random vector b with the diagonal covariance matrix, as shown in Eqs. (55)-(57), and assume the random vector b is uniform random vector. Only 50 highest variables are used in computations. The MCS, the OSPFEM and the proposed MSPFEM are employed for this problem. The MCS is performed with 10000 samples. The MCS solutions computed with 10000 samples are seen as the reference solutions to test the accuracy of the proposed method.

The means and standard deviations of the temperatures T(0,0,t) at t =0.5, 1, 1.5 and 2 computed using the MCS, the OSPFEM and the proposed approach are compared in Table 9. The CPU times for all methods are shown in Table 10. Figures 6(a) and 6(b) show the temperature contours of the means and standard deviations at t = 2, respectively. It can be seen from Tables 9 and 10 that the MSPFEM cost almost as much CPU time as the OSPFEM to produce more accurate evaluations.

Figure 6
Contours of temperature distribution at t = 2: (a) means and (b) standard deviations.

Table 9
Means and standard deviations of the temperatures T (0, 0, t) at different times.

	Means		Standard deviations
t	0.5	1	1.5	2	0.5	1	1.5	2
MCS	8.076	8.928	9.016	9.026	0.807	1.169	1.231	1.241
MSPFEM	8.075	0.926	9.013	9.023	0.807	1.165	1.226	1.235
Percentage errors	0.016%	0.028%	0.034%	0.035%	-0.017%	0.298%	0.401%	0.428%
OSPFEM	8.075	8.923	9.008	9.017	0.799	1.131	1.181	1.187
Percentage errors	0.017%	0.058%	0.084%	0.093%	0.959%	3.319%	4.255%	4.502%

Table 10
CPU times (s) for example 4.4.

MCS	MSPFEM	SGM
2877.16	15.29	15.92

5 CONCLUSIONS

In actual engineering, there is a considerable number of structures with low-level uncertainties that can be analyzed using the SPFEM in theory. However, the traditional computational scheme depends on the type of problem considered. Namely, each type of problem has its own computational scheme, e.g., the computational scheme of the static problems is far different from that of the eigenvalue problems. Furthermore, the derivatives of the system matrices with respect to the random variables are often required, and they may be very complicated for some problems. In this paper, an unified computational scheme of the SPFEM is developed for linear or nonlinear structures with correlated or uncorrelated random variables. The proposed computational scheme can provide second-order estimates of the mean and variance without using the derivatives of the system matrices; therefore, compared with the traditional computational scheme of the SPFEM, the proposed scheme is more convenient to implement for a complex stochastic system. In the case of one random variable with a symmetric probability density function, the proposed computational scheme can even provide a result with fifth-order accuracy. When the proposed scheme is used, it involves 2q + 1 (q is the number of random variables) analyses of deterministic systems. The numerical examples show that the proposed scheme can be used for complicated problems and provides very accurate results.

Acknowledgements

The authors are grateful for the support of the Natural Science Foundation of China (No. 11272076), the Fundamental Research Funds for the Central Universities (No. DUT13LK12) and the Key Laboratory Fund of Liaoning (No. L2013019).

References

• Boyaval S. E. B. (2012). A fast Monte-Carlo method with a reduced basis of control variates applied to uncertainty propagation and Bayesian estimation. Computer Methods in Applied Mechanics and Engineering, 241:190-205.
• Chang T., Liu M., Chang H. (2008). Finite element analysis of nonlinear shell structures with uncertain material property. Thin-Walled Structures, 46:1055-1065.
• Deb M. K., Babuska I. M., Oden J. T. (2001). Solution of stochastic partial differential equations using Galerkin finite element techniques. Computer Methods in Applied Mechanics and Engineering, 190:6359-6372.
• Falsone G., Ferro G. (2007). An exact solution for the static and dynamic analysis of FE discretized uncertain structures. Computer Methods in Applied Mechanics and Engineering, 196:2390-2400.
• Falsone G., Impollonia N. (2002). A new approach for the stochastic analysis of finite element modelled structures with uncertain parameters. Computer Methods in Applied Mechanics and Engineering, 191:5067-5085.
• Ganapathysubramanian B., Zabaras N. (2007). Sparse grid collocation schemes for stochastic natural convection problems. Journal of Computational Physics, 225:652-685.
• Holtz M. (2011). Sparse Grid Quadrature in High Dimensions with Applications in Finance and Insurance. Springer Berlin Heidelberg,
• Johnson E. A., Proppe C., Spencer Jr B. F., Bergman L. A., Székely G. S., Schuёller G. I. (2003). Parallel processing in computational stochastic dynamics. Probabilistic Engineering Mechanics, 18:37-60.
• Kamiński M. (2006). On generalized stochastic perturbation-based finite element method. Communications in Numerical Methods in Engineering, 22:23-31.
• Kamiński M. (2007). Generalized perturbation-based stochastic finite element method in elastostatics. Computers and Structures, 85:586-594.
• Kamiński M. (2010). Potential problems with random parameters by the generalized perturbation-based stochastic finite element method. Computers and Structures, 88:437-445.
• Kamiński M. (2011). On semi-analytical probabilistic finite element method for homogenization of the periodic fiber-reinforced composites. International Journal for Numerical Methods in Engineering, 86:1144-1162.
• Kamiński M. (2013). The Stochastic Perturbation Method for Computational Mechanics: Practical Applications in Science and Engineering. John Wiley & Sons, New Jersey.
• Kamiński M., Szafran J. (2009). Eigenvalue analysis for high telecommunication towers with lognormal stiffness by the response function method and SFEM. Computer Assited Mechanics and Engineering Sciences, 16:279-290.
• Kim H. M., Inoue J. (2004). A spectral stochastic element free Galerkin method for the problems with random material parameter. International Journal for Numerical Methods in Engineering, 61:1957-1975.
• Kleiber M., Hien T. D. (1992). The stochastic finite element method: basic perturbation technique and computer implementation. Wiley,
• Lei Z., Qiu C. (2000). Neumann dynamic stochastic finite element method of vibration for structures with stochastic parameters to random excitation. Computers and Structures, 77:651-657.
• Li J., Chen J. B. (2004). Probability density evolution method for dynamic response analysis of structures with uncertain parameters. Computational Mechanics, 34:400-409.
• Lin J. H., Zhao Y., Zhang Y. H. (2001). Accurate and highly efficient algorithms for structural stationary/non-stationary random responses. Computer Methods in Applied Mechanics and Engineering, 191:103-111.
• Liu W. K., Belytschko T., Mani A. (1986). Random field finite elements. International Journal for Numerical Methods in Engineering, 23:1831-1845.
• Papadrakakis M., Kotsopulos A. (1999). Parallel solution methods for stochastic finite element analysis using Monte Carlo simulation. Computer Methods in Applied Mechanics and Engineering, 168:305-320.
• Rahman S. (2008). A polynomial dimensional decomposition for stochastic computing. International Journal for Numerical Methods in Engineering, 76:2091-2116.
• Sakata S., Ashida F., Enya K. (2011). Perturbation-based stochastic microscopic stress analysis of a particle-reinforced composite material via stochastic homogenization analysis considering uncertainty in material properties. International Journal for Multiscale Computational Engineering, 9:395-408.
• Sakata S., Ashida F., Shimizu Y. (2011). Inverse stochastic homogenization analysis for a particle-reinforced composite material with the Monte Carlo simulation. International Journal for Multiscale Computational Engineering, 9:409-423.
• Stefanou G. (2009). The stochastic finite element method: past, present and future. Computer Methods in Applied Mechanics and Engineering, 198:1031-1051.
• Székely G. S., Schuёller G. I. (2001). Computational procedure for a fast calculation of eigenvectors and eigenvalues of structures with random properties. Computer Methods in Applied Mechanics and Engineering, 191:799-816.
• Van den Nieuwenhof B., Coyette J. (2003). Modal approaches for the stochastic finite element analysis of structures with material and geometric uncertainties. Computer Methods in Applied Mechanics and Engineering, 192:3705-3729.
• Wang Q., Yao J. (2010). Study on high order perturbation-based nonlinear stochastic finite element method for dynamic problems. Journal of Marine Science and Application, 9:386-392.
• Wu F., Zhong W. X. (2013). A hybrid approach for the time domain analysis of linear stochastic structures. Computer Methods in Applied Mechanics and Engineering, 265:71-82.
• Xia B. Z., Yu D. J. (2012). Modified sub-interval perturbation finite element method for 2D acoustic field prediction with large uncertain-but-bounded parameters. Journal of Sound and Vibration, 331:3774-3790.
• Xiu D., Karniadakis G. E. (2002). The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal On Scientific Computing, 24:619-644.
• Zhang H. W., Xu W. L., Di S. L., Thomson P. F. (2002). Quadratic programming method in numerical simulation of metal forming process. Computer Methods in Applied Mechanics and Engineering, 191:5555-5578.

Publication Dates

• Publication in this collection
  2015

History

• Received
  12 Dec 2014
• Reviewed
  30 Apr 2015
• Accepted
  06 Aug 2015