- Email: [email protected]
An efficient method for evaluating the natural frequencies of structures with uncertain-but-bounded parameters
Computers and Structures 87 (2009) 582–590
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
An efficient method for evaluating the natural frequencies of structures with uncertain-but-bounded parameters Su Huan Chen *, Liang Ma, Guang Wei Meng, Rui Guo Department of Mechanics, Nanling Campus, Jilin University, Changchun 130025, PR China
a r t i c l e
i n f o
Article history: Received 30 August 2008 Accepted 24 February 2009 Available online 2 April 2009 Keywords: Efficient method Natural frequencies Uncertain-but-bounded parameters Epsilon-algorithm
a b s t r a c t Using interval theory and the second-order Taylor series, the eigenvalue problems of structures with multi-parameter can be transformed into those with single parameter. The epsilon-algorithm is used to accelerate the convergence of the Neumann series to obtain the bounds of eigenvalues of structures with single interval parameter, thus increasing the computing accuracy and reducing the computational effort. Finally, the effect of uncertain parameters on natural frequencies is evaluated. Two engineering examples show that the proposed method can give better results than those obtained by the first-order approximation, even if the uncertainties of parameters are fairly large. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction The natural frequency analysis for structures with deterministic parameters has been extensively developed. However, in engineering situations, the structural parameters are often uncertain, such as the inaccuracy of the measurement, errors in the manufacturing and assembly process, invalidity of some components and uncertainty in boundary conditions etc. The uncertainties of structural parameters may lead to large and unexpected excursion of responses that may lead to drastic reduction in accuracy and precision of the operation. Therefore, uncertainty plays an important role in the modern engineering structural analysis. Over the past decades, a number of methods have been developed that include uncertain model properties in the finite element (FE) analysis and aim at the quantification of the uncertainty on the analysis result. The probabilistic concept is by far the most popular method for numerical uncertainty modeling. Its popularity has led to a number of probabilistic FE procedures [1–3]. However, probabilistic modeling is not the only way to describe the uncertainty, and uncertainty is not equal to randomness. Indeed, the probabilistic approaches are not able to deliver reliable results at the required precision without sufficient experimental data to validate regarding the joint probability densities of the random variables or functions involved. Therefore, one may recognize that uncertainties in parameters can be modeled on the basis of alterative, non-probabilistic conceptual frameworks. One such approach, based on a set theoretic formulation, is an unknown-but-bounded model (convex models). Such set models of uncertainty have been * Corresponding author. Tel.: +86 431 8509 5259; fax: +86 431 8509 5288. E-mail address: [email protected] (S.H. Chen). 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.02.009
applied by Deif in linear programming and system theory [4]. Recently, such set models of uncertainties in parameters have drawn interest both from the system control robustness analysis field and from the structural failure measures field. For example, the convex model was introduced by Ben-Haim and Elishakoff [5], discussed later by Lindberg [6], for the study of dynamic response and failure of structures under pulsed parametric loading; the convex model has been applied in determining the upper and lower bounds of static response for structures by Liu et al. [7]. The convex model has also been applied to impulsive response, buckling analysis and optimal design of structures with uncertain parameters [8–11]. Since the mid-1960s, a new method called the interval analysis has appeared. Moore [12], Alefeld and Herzberger [13] have done the pioneering work. Mathematically, linear interval equations, nonlinear interval equation, and interval eigenvalue problems have been partially resolved. But because of the complexity of the algorithm, it is difficult to apply these results to practical engineering problems. Recently, the interval finite element (IFE) method was presented by Chen et al. [14] which makes the method easier to deal with the interval eigenvalues for closed-loop systems of structures with uncertain parameters. It should be pointed out that the previous interval analysis methods for computing the upper and lower bounds of response of structures with interval parameters are based on the first-order perturbation or the first-order Taylor series. Although the procedure is easy to implement on the computer and incorporates the finite element code, the applications of the previous methods are limited to the case when the interval uncertainties of parameters are small. However, if the parameter uncertainties are fairly large or the combination of a large number uncertain parameters, the
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
accuracy of the computational results will become unacceptable. Thus, it is highly desirable to present a more accurate method for computing the upper and lower bounds of responses of structures with fairly large uncertainties of interval parameters. To this end, this paper presents an efficient method to estimate the natural frequencies of structures for the case with fairly large uncertainties of parameters. As we known, the eigenvalues k is defined as x2 in mode analysis of FEM problems, where x denotes the natural frequencies of structures. The idea of the proposed method is that the eigenvalues are considered as functions of the structural parameters; using interval theory and the second-order Taylor series expansion, we transform the eigenvalue problems of structures with multi-parameter into those with single parameter; The epsilon-algorithm is used to accelerate the convergence of Neumann series to compute the eigenvalues of structures with single interval parameter, thus the interval eigenvalues of the structures with multi-parameter are obtained and the effect of uncertain parameters on natural frequencies is evaluated finally. This paper is organized as follows. A brief review of mathematical background on the interval analysis is given in the Appendix. In the Section 2, the epsilon-algorithm for eigensolution reanalysis of structures with fairly large changes of parameters is discussed [15]. The Section 3 presents the definition of the interval eigenvalue problems of the interval parameter structures. In the Section 4, the proposed method for evaluating the natural frequencies of structures with fairly large uncertainties of parameters is developed. In Section 5, two engineering examples are given to illustrate the application of the proposed method. The results obtained by the proposed method are compared with those obtained by the exact solutions and the first-order approximation. The conclusions are drawn in Section 6.
where K0 and M0 are the stiffness and mass matrices of the finite element assemblage. This equation will be referred to as the initial problem in the following discussions. The perturbation method studies the changes of eigenvalues of the system subjected to changes in its design parameters. Therefore, if the initial system is represented by Eq. (5), the problem becomes that to determine u and k when K0 and M0 are perturbed to the form K0 þ DK and M0 þ DM, respectively.
Ku ¼ kMu
ð6Þ
where
K ¼ K0 þ DK
ð7Þ
M ¼ M0 þ DM
ð8Þ
in Eqs. (6)–(8), if we introduce the following notations
f 0 ¼ k0 M0 u0
ð9Þ
f ¼ kMu
ð10Þ
Df ¼ kMu k0 M0 u0
ð11Þ
then we obtain
ðK0 þ DKÞu ¼ f 0 þ Df:
u ¼ ðK0 þ DKÞ1 ðf 0 þ DfÞ 1 1 ¼ ðI þ K1 0 DKÞ K0 ðk0 M0 u0 þ k0 DMu0 Þ
¼ ðI þ BÞ1 K1 0 ðk0 M0 u0 þ k0 DMu0 Þ
2.1. The epsilon-algorithm
then obtaining a series
ðjÞ kþ1
e
ð1Þ
ðjþ1Þ k1
¼e
þ
h
ðjþ1Þ k
e
ðjÞ k
e
i1
ð2Þ j; k ¼ 0; 1; 2; . . .
ð3Þ
The iterative formulae (1)–(3) are similar to the scalar case except that it requires the inverse of a vector. The definition of the inverse of the vector is given by Roberts [19]
u1 ¼
u u ¼ Pn H 2 ðu uÞ i¼1 jui j
ð4Þ
ð13Þ
In Eq. (13), B ¼ K1 0 DK; k and u were approximated by k0 and u0 , respectively. By using the Neumann series expansion, we have
~ ðI B þ B2 ÞK1 u 0 ðk0 M0 u0 þ k0 DMu0 Þ
eðjÞ 1 ¼ 0 ðjÞ e0 ¼ sj
ð12Þ
It follows that
2. Eigensolution reanalysis with the epsilon-algorithm
The epsilon-algorithm was presented [16–18] to accelerate the convergence of a sequence. Given a vector sequence fs0 ; s1 ; s2 ; . . .g, we construct the iterative form to obtain the vector sequence in the epsilon-algorithm as follows:
583
ð14Þ
~ 0 ¼ K1 u 0 ðk0 M0 u0 þ k0 DMu0 Þ 1 1 ~ 1 ¼ K1 ~ u 0 DKK0 ðk0 M0 u0 þ k0 DMu0 Þ ¼ K0 DKu0
~ ~ 2 ¼ K1 u 0 DKu1 .. . e ~ s ¼ K1 u 0 DK u s1
ð15Þ
s ¼ 3; 4; . . .
Assume the solution of the Eq. (6) has the following form
~1 þ þ u ~s þ ~¼u ~0 þ u u
ð16Þ
~ 0 ; s1 ¼ We define a vector sequence fs0 ; s1 ; . . . ; ss ; . . .g where s0 ¼ u ~ 1 ; s2 ¼ u ~0 þ u ~1 þ u ~ 2 , and in general ~0 þ u u
si ¼
i X
~j u
i ¼ 0; 1; 2; . . . ; s; . . .
ð17Þ
j¼0
where the u denotes the conjugate of u and uH the conjugate and transpose of u. It should be noted that for the undamped structure the modal vectors are real. Thus, u in Eq. (4) should be replaced by u, and uH should be replaced by uT . In this case, u1 in Eq. (4) is just a normalized version of u itself. The vector epsilon-algorithm table can be constructed by Eqs. (1)–(3).
To accelerate the convergence of the sequence (17), using the vector epsilon-algorithm (Eqs. (1)–(3)), the sequence, fs0 ; s1 ; . . . ; ss ; . . .g, can be obtained. Refs. [16–18] have pointed out that in ðjÞ the epsilon iterative formulae, when i is odd, the ei are meaningðjÞ less, and that when the i ¼ 2k ðk ¼ 1; 2; . . .Þ, the vector e2k is the Shanks transformation. The solution is the last even row in the epsilon-algorithm,
2.2. The epsilon-algorithm for accelerating the convergence of Neumann series
u ¼ e2k
ðjÞ
ð18Þ
Compute the eigenvalues using the Rayleigh quotients The eigenproblem for the structure is as follows
K0 u0 ¼ k0 M0 u0
ð5Þ
ki ¼
uTi Kui uTi Mui
ð19Þ
584
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
3. Generalized interval eigenvalue problems 3.1. The definition of generalized interval eigenvalue problems The eigenvalue problem is described by equation
KðbÞu ¼ kMðbÞu
ð20Þ
under the parameter constraint
b6b6b
ð21Þ
Using the interval matrix expression and the interval extension, Eq. (20) can be expressed as I
I
Kðb ÞuI ¼ kI Mðb ÞuI
ð22Þ
Eq. (22) is called interval eigenvalue problem, under constraint (21). I The basic problem is, for given interval parameter vector b , to I find the interval eigenvalue ki , which is the smallest interval and encloses all possible eigenvalue ki , satisfying KðbÞu ¼ kMðbÞu. In other words, we seek a hull, i.e.
n
C ¼ ki : ki 2 R; KðbÞu ¼ ki MðbÞu; ki –0; b 2 bI
o
ð23Þ
to the set
kIi ¼ ki ; ki
ð24Þ
where I
b2b
ki ¼ max ki ðhKðbÞ; MðbÞiÞ ¼ max ki ðbÞ;
ð25Þ I
b2b
ð26Þ
where ki ðhKðbÞ; MðbÞiÞ means solving the eigenvalues for the generalized eigenvalue problems with KðbÞ and MðbÞ described by Eq. (20). 3.2. The interval eigenvalues of structures Using Eq. (A1), interval parameter vector can be written as
h i I c c b ¼ b; b ¼ b Db; b þ Db
ð27Þ
where
Db ¼ ðDb1 ; Db2 ; . . . ; Dbm ÞT
ð28Þ
Eqs. (27) and (28) mean that the uncertain parameters are constrained into a m-dimension rectangle with widths 2Db1 ; 2Db2 ; . . . ; 2Dbm . For a m-dimension rectangle, we can define its boundary points P bd which are placed on the boundary surfaces and the extreme points P ex which are placed at the intersections of the boundary surfaces. For example, we consider a 2-dimension rectangle shown in Fig. 1. The boundary points P bd are placed on P1ex P2ex ; P1ex P 3ex ; P3ex P4ex and P2ex P4ex , and the extreme points are the intersections P 1ex ; P 2ex ; P 3ex and P 4ex . It means that the extreme points Psex can be determined by the parameter vector
¼ Pex ðb1 ; b2 ; . . . ; bm Þ : bk ¼
k ¼ 1; 2; . . . ; m s ¼ 1; 2; . . . ; 2
on the surface’s boundary line and the extrema of any boundary line will occur on the line’s extreme points. According to the definitions of the extreme points Pex and the boundary points P bd , the extrema of Eqs. (25) and (26) will occur on the extreme points Psex ; s ¼ 1; 2; . . . ; 2m , thus we obtain the upper and lower bounds of the interval eigenvalues.
ki ¼ min ki ðPsex Þ;
i ¼ 1; 2; . . . ; n s ¼ 1; 2; . . . 2m
ð30Þ
max ki ðPsex Þ;
i ¼ 1; 2; . . . ; n s ¼ 1; 2; . . . 2m
ð31Þ
ki ¼
ki ¼ min ki ðhKðbÞ; MðbÞiÞ ¼ min ki ðbÞ;
Psex
Fig. 1. 2-Dimension rectangle.
c bk m
þ Dbk or
c bk
Dbk
From Eqs. (30) and (31), it can be seen that in order to obtain the exact solutions of the upper and lower bounds of eigenvalues, we have to solve 2m eigenvalue problems given by Eq. (20) on the extreme points P sex ðs ¼ 1; 2; . . . ; 2m Þ of the m-dimension rectangle, then take min ki ðPsex Þ as the lower bounds of eigenvalues ki ði ¼ 1; 2; . . . ; nÞ, and take max ki ðPsex Þ as the upper bounds of eigenvalues ki ði ¼ 1; 2; . . . ; nÞ. This is a tedious job which would not be acceptable if m is too large. Therefore, we need to seek an efficient method for computing the upper and lower bounds of eigenvalues of structures with interval parameters. 4. Methods for evaluating the natural frequencies of engineering structures with uncertain parameters As we known, the eigenvalues k is defined as x2 in mode analysis of FEM problems, where x denotes the natural frequencies of structures. If the bounds of eigenvalues are obtained, the natural frequencies of structures are evaluated correspondingly. 4.1. The method based on the first-order approximation Using the first-order Taylor series and the first-order perturbation for eigenvalue analysis, Ref. [21] developed a method for computing the bounds of eigenvalues as follows
c I c k b ¼ k b Dk; k b þ Dk
ð32Þ
ð29Þ m
where m is the number of the interval parameters, 2 denotes the number of the extreme points. Eqs. (25) and (26) call for finding the minimum and maximum I of the function ki ðbÞ in the interval set b . According to optimal theory [20], the extrema of Eqs. (25) and (26) will occur on the boundary surfaces of the m-dimension rectangle described by Eqs. (27) and (28). Similarly, the extrema of any boundary surface will occur
where
Dk ¼
m
X c
c T @KðbÞ
u b
c Db i u b
@b i b ¼b i¼1 i i
c c T @MðbÞ
c
k b u b Db i u b
@bi bi ¼bc
ð33Þ
i
c
I
in which bi and Dbi denote the mean (or midpoint) value of bi and I the uncertainty (or the maximum width) in bi , respectively.
585
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
4.2. The proposed method for evaluating the natural frequencies of structures with fairly large uncertainties of the parameters Using the second-order Taylor series to expand ki ðb1 ; b2 ; . . . ; bm Þ c c c c about the point b ¼ b1 ; b2 ; . . . ; bm , we have
c c c 1 c c T c ki ðbÞ ¼ ki b þ GTki b b b þ b b Hki b b b 2
ð34Þ
c c where b and Hki b are the gradient vector and the Hessian c matrix of eigenvalues at the point b , respectively. The gradient vector and the Hessian matrix can be written as follows GTki
GTki ðbÞ ¼
@ki @ki @ki @b1 @b2 @bm
2
@ 2 ki @b1 @b1
@ 2 ki @b1 @b2
6 6 @ 2 ki 6 Hki ðbÞ ¼ 6 @b2 @b1 6 4
@ 2 ki @b2 @b2
@ 2 ki @bm @b1
@ 2 ki @bm @b2
ð35Þ @ 2 ki @b1 @bm
3
7 7 7 7 7 5
@ 2 ki @b2 @bm
ð36Þ
@ 2 ki @bm @bm
c @ki
c c c
ki b1 ; b2 ; bm ¼ ki b þ b1 b1
@b1 b1 ¼bc 1
2 1 @ ki
c 2 þ b1 b1
2 @b1 @b1
c
2
I I I I c c c I c c c I ki b1 ; b2 ; b3 ki b1 ; b2 ; b3 þ ki b1 ; b2 ; b3 þ ki b1 ; b2 ; b3 c c c 2ki b1 ; b2 ; b3
ð42Þ
As we have discussed in the Eqs. (30) and (31), exact solutions I I I for ki b1 ; b2 ; b3 can be obtained by solving 23 finite element prob-
I I I ki b1 ; b2 ; b3 " " # ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ¼ min ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ " ## ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ð43Þ Similarly, we have that for 1-dimension case
1 @ ki
c ðb2 b2 Þ þ
2 @b2 @b2
b2 ¼bc 2
2
c 2 b2 b2
b2 ¼bc2
ð38Þ Substituting Eqs. (37) and (38) into Eq. (34), we have
c c c c ki ðbÞ ¼ ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b c 1 c T c þ b b HLki b b b 2
ð41Þ
From Eq. (41), it can be seen that the computational problems of the interval eigenvalues of structures with multi-parameter have been transformed into those with single parameter. The equation can be written as follows while m = 3
ð37Þ
b1 ¼b1
c @ki
c c ki b1 ;b2 ; .. .bm ¼ ki ðb Þ þ @b
I I c c c I c ki b ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c I þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b
lems, that is
The computation of the Hessian matrix is a tedious job if m is c large. For the sake of simplicity, we consider b2 ¼ b2 ; b3 ¼ c c b3 ; . . . ; bm ¼ bm in Eq. (34), so ki ðbÞ becomes a function of b1 . We have
Similarly, we obtain
Because of the complexities of solving mixed partial derivatives and the widening of computed intervals caused by the multiplications and accumulations of different interval parameters, we neglect the last part of the right side of Eq. (40) and have
ð39Þ
h h i I c c c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð44Þ
h h i c I c c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð45Þ
h h i c c I c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð46Þ
Substituting Eqs. (44)–(46) into Eq. (42), we have
h i h i 3 c c c c c c c c min ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 þ min ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 7 6 h 7 6 i c c c c c c c 7 6 ; k 2k ; þ min k b ; b ; b b ; b ; b b ; b ; b 3 3 i i i 7 6 1 2 1 2 1 2 3 I I I h i h i 7 ki b1 ; b2 ; b3 6 7 6 c c c c c c c c 6 max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 þ max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 7 7 6 5 4 h i c c c c c c c þ max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 2ki b1 ; b2 ; b3 2
extension of the Using the natural interval function, yields I I c c c I c ki b ¼ ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c I þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b T 1 I c c I c þ b b HLki b b b ð40Þ 2
c where HLki b denotes the Hessian matrix in which the diagonal elements are zero.
ð47Þ
From Eq. (47), it can be seen that, using the proposed method, the problem for solving 23 FE equations has been transformed into that for solving 2 3 FE equations. In order to further reduce the computation effort, we suggest that based on the Neumann series expansions, the accurate approximations of solutions for the upper and lower bounds of eigenvalues can be obtained by using the reanalysis method with the epsilon-algorithm given by the Section 2 for structures with a single interval parameter.
586
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
In general, the exact solution for the interval eigenvalues can be obtained by solving 2m FE problems, and using proposed method the approximate one can be obtained by 2m reanalysis problems, thus reducing the computational effort. Compared with the first-order approximation, the proposed method replaces the first-order perturbation with the reanalysis method based on epsilon-algorithm which considers many higher order items to make the result more accurate. Although we neglect the mixed second-order items of Taylor series in Eq. (40), both the first-order and the second-order items of each parameter are included. However, only the first-order items of Taylor series are considered in the first-order approximation. The complete method presented in this paper is summarized below: 1. Transform the interval eigenvalues of structures with multiparameter into those of structures with single parameter using Eq. (41); 2. Compute the upper and lower bounds of eigenvalues of structures with a single interval parameter using the reanalysis method given by Eqs. (15) and (1)–(3); 3. Using Eq. (41), the interval eigenvalues of structures with multi-parameter are obtained so that the natural frequencies are evaluated finally.
5. Numerical examples 5.1. Example 1 We consider a 40-story structure as shown in Fig. 2. The finite element model of the given structure consists of 202 nodes, 357 beam elements and 1212 degrees of freedom. The parameters of the structure are as follows. The length of the longitudinal beam is 5.0 m and that of cross beam is 3.0 m; The Young’s modulus of the used material is E ¼ 2:1 1011 Pa; the mass density q ¼ 7:8 103 kg=m3 ; the height and width of the cross-section of the beams are H ¼ 0:1 m, and W ¼ 0:2 m; Assume that E; H and W are the interval parameters given by
EI ¼ Ec þ aEc e HI ¼ Hc þ bHc e W I ¼ W c þ cW c e e ¼ ½1; 1 The results are listed by frequencies. Case 1 (one parameter). Consider the case with one parameter with uncertainty (a ¼ 5%). The upper and lower bounds of the frequencies of the structure obtained by the firstorder approximation are listed in Table 1. The corresponding results obtained by the proposed method are listed in Table 2. For comparison, the exact solutions of the bounds of frequencies obtained by Eq. (43) are also listed in Tables 1 and 2. From results in Tables 1 and 2, it can be seen that the max error of the upper bound of the mode 2 for the first-order approximation is 1.3%, while the corresponding error is 0.0% for the proposed method. Case 2 (two parameters). Consider the case with two parameters with uncertainties (a ¼ b ¼ 5%Þ. The results are listed in Tables 3 and 4, and show that the max error obtained by the first-order approximation of the lower bound of the mode 5 reaches 4.4%, while the corresponding error obtained by the proposed method is 0.0%. Case 3 (three parameters). Consider the case with three parameters with uncertainties (a ¼ b ¼ c ¼ 5%Þ. The results are listed in Tables 5 and 6, and show that the max error obtained by the first-order approximation for the lower bound of the mode 5 reaches 4.4%, while the corresponding error obtained by the proposed method is 0.0%. The max error of the lower bound of the mode 3 is 0.6%.If the uncertainties are a ¼ b ¼ c ¼ 10%,
Fig. 2. 40-Story frame structure.
Table 1 Bounds of frequencies using the first-order approximation (a ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.29D03, [2.11D02, [4.98D02, [6.00D02, [1.19D01,
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
1.2 0.9 1.2 1.0 0.8
0.9 1.3 0.9 0.9 0.8
3.53D03] 2.27D02] 5.35D02] 6.43D02] 1.27D01]
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
587
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590 Table 2 Bounds of frequencies using the proposed method (a ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
Table 3 Bounds of frequencies using the first-order approximation (a ¼ b ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.03D03, [1.95D02, [4.68D02, [5.53D02, [1.09D01,
[3.16D03, [2.03D02, [4.84D02, [5.76D02, [1.14D01,
4.1 3.9 3.3 4.0 4.4
2.2 2.1 2.2 2.2 2.3
3.75D03] 2.41D02] 5.62D02] 6.84D02] 1.35D01]
3.67D03] 2.36D02] 5.50D02] 6.69D02] 1.32D01]
Table 4 Bounds of frequencies using the proposed method (a ¼ b ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.15D03, [2.02D02, [4.82D02, [5.74D02, [1.14D01,
[3.16D03, [2.03D02, [4.84D02, [5.76D02, [1.14D01,
0.3 0.5 0.4 0.3 0.0
0.3 0.4 0.2 0.1 0.0
3.66D03] 2.35D02] 5.49D02] 6.68D02] 1.32D01]
3.67D03] 2.36D02] 5.50D02] 6.69D02] 1.32D01]
Table 5 Bounds of frequencies using the first-order approximation (a ¼ b ¼ c ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.03D03, [1.95D02, [4.59D02, [5.53D02, [1.09D01,
[3.16D03, [2.03D02, [4.79D02, [5.76D02, [1.14D01,
4.1 3.9 4.2 4.0 4.4
2.5 2.1 2.3 2.2 2.3
3.76D03] 2.41D02] 5.69D02] 6.84D02] 1.35D01]
3.67D03] 2.36D02] 5.56D02] 6.69D02] 1.32D01]
Table 6 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.15D03, [2.02D02, [4.76D02, [5.74D02, [1.14D01,
[3.16D03, [2.03D02, [4.79D02, [5.76D02, [1.14D01,
0.3 0.5 0.6 0.3 0.0
0.3 0.4 0.4 0.1 0.0
3.66D03] 2.35D02] 5.54D02] 6.68D02] 1.32D01]
3.67D03] 2.36D02] 5.56D02] 6.69D02] 1.32D01]
Table 7 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 10%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[2.88D03, [1.85D02, [4.32D02, [5.24D02, [1.04D01,
[2.91D03, [1.87D02, [4.41D02, [5.31D02, [1.05D01,
1.0 1.1 2.0 1.3 1.0
0.8 0.8 1.2 0.8 0.7
3.91D03] 2.51D02] 5.89D02] 7.12D02] 1.41D01]
3.94D03] 2.53D02] 5.96D02] 7.18D02] 1.42D01]
588
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
The results are listed by frequencies.
from Table 7 it can be seen that the max error of the lower bound of the mode 3 obtained by the proposed method is only 2.0%.
Case 1 (one parameter). Consider the one parameter case with uncertainty (a ¼ 3%Þ. The upper and lower bounds of the frequencies of the structure obtained by the first-order approximation are listed in Table 8. The corresponding results obtained by the proposed method are listed in Table 9. For comparison, the exact solutions of the bounds of frequencies obtained by Eq. (43) are also listed in Tables 8 and 9. From results in Tables 8 and 9, it can be seen that the max error of the upper bound of the mode 2 for the first-order approximation is 10.5%, while the corresponding error is 0.0% for the proposed method. Case 2 (two parameters). Consider the case with two parameters with uncertainties (a ¼ b ¼ 3%Þ. The
5.2. Example 2 We consider a chassis structure as shown in Fig. 3. The finite element model of the given structure consists of 20,953 nodes, 19,880 thin shell elements and 125,718 degrees of freedom. The parameters of the structure are as follows. The length of the chassis is 3.0 m, the front width of the chassis 0.9 m, the back width 1.2 m; The Young’s modulus of the material used is E ¼ 2:1 1011 Pa; the mass density q ¼ 7:8 103 kg=m3 ; the thickness of the thin shell in the longitudinal beam is t1 ¼ 0:01 m, and that of the cross beam is t2 ¼ 0:01 m. Assume that E; t 1 and t 2 are the interval parameters given by
EI ¼ Ec þ aEc e t I1 ¼ t c1 þ bt c1 e tI2 ¼ t c2 þ ct c2 e e ¼ ½1; 1
Fig. 3. Chassis structure.
Table 8 Bounds of frequencies using the first-order approximation (a ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.26D+00, [1.88D+01, [3.38D+01, [3.39D+01, [4.20D+01,
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
1.3 10.5 0.6 6.4 5.8
1.3 8.8 0.6 5.6 5.2
8.73D+00] 2.36D+01] 3.52D+01] 3.94D+01] 4.83D+01]
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
Table 9 Bounds of frequencies using the proposed method (a ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
Table 10 Bounds of frequencies using the first-order approximation (a ¼ b ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.03D+00, [1.75D+01, [3.37D+01, [3.36D+01, [4.11D+01,
[8.25D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
2.7 16.3 0.6 7.2 7.4
2.2 12.8 0.6 6.1 6.3
8.94D+00] 2.46D+01] 3.53D+01] 3.97D+01] 4.90D+01]
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
589
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590 Table 11 Bounds of frequencies using the proposed method (a ¼ b ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.24D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
[8.25D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
0.1 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
Table 12 Bounds of frequencies using the first-order approximation (a ¼ b ¼ c ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[7.92D+00, [1.39D+01, [3.32D+01, [3.31D+01, [3.99D+01,
[8.20D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
3.4 32.2 1.2 8.3 9.5
2.7 20.7 0.8 7.2 8.0
9.04D+00] 2.68D+01] 3.58D+01] 4.02D+01] 5.00D+01]
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
Table 13 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.19D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
[8.20D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
0.1 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
Table 14 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 10%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[7.46D+00, [1.84D+01, [3.13D+01, [3.43D+01, [4.14D+01,
[7.51D+00, [1.87D+01, [3.14D+01, [3.43D+01, [4.15D+01,
0.7 1.6 0.3 0.0 0.2
0.4 0.4 0.5 0.0 0.2
9.49D+00] 2.40D+01] 3.75D+01] 3.90D+01] 4.87D+01]
results are listed in Tables 10 and 11, and show that the max error obtained by the first-order approximation of the lower bound of the mode 2 reaches 16.3%, while the corresponding error obtained by the proposed method is 0.0%. Case 3 (three parameters). Consider the case with three parameters with uncertainties (a ¼ b ¼ c ¼ 3%Þ. The results are listed in Tables 12 and 13, and show that the max error obtained by the first-order approximation for the lower bound of mode 2 reaches 32.2%, while the corresponding error obtained by the proposed method is 0.0%. The max error of the lower bound of mode 1 is 0.1%.If the uncertainties of the parameters are a ¼ b ¼ c ¼ 10%, from Table 14 it can be seen that the max error of the lower bound of the mode 2 obtained by the proposed method is only 1.6%.
9.53D+00] 2.41D+01] 3.77D+01] 3.90D+01] 4.88D+01]
The above numerical results indicate that the first-order approximation is valid for the case when the uncertainties of parameters are small. The proposed method can give good results even if the uncertainties of parameters are fairly large. The errors of the results obtained by the proposed method are attributed to neglecting of the mixed second-order items of Taylor series in Eq. (40). So the corresponding errors will go up as the uncertainties become very large, but they are much smaller than those obtained by the first-order approximation.
6. Conclusions An efficient method for estimating the natural frequencies of structures with fairly large uncertainties of parameters is developed. Based on interval theory and the second-order Taylor series, the solution for the eigenvalues have been reduced from 2m finite element analysis problems to 2m eigenvalue reanalysis ones. The epsilon-algorithm was used to accelerate the convergence of the Neumann series to compute the bounds of eigenvalues of structures with single interval parameter, thus increasing the computational accuracy and reducing the computing effort for the interval
590
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
eigenvalues of structures with multi-parameter. Two numerical examples were given to illustrate the application and validity of the proposed method. The three cases, one parameter, two parameters and three parameters were discussed. The comparison of numerical results obtained by the exact solution, the first-order approximation and the proposed method was given. The results indicated that the first-order approximation is valid only for the case when the uncertainties of parameters are small and the proposed method can give good results even if the uncertainties of parameters are fairly large. Acknowledgements This work is supported by Science and Technology Development Foundation of Jilin Province (20070541), 985 Engineering of Jilin Univ. and Innovation fund for 985 Engineering of Jilin Univ. No. 20080104. Appendix A. The definitions of the interval and interval operations Assuming that R is the set of all real numbers, IðRÞ; IðRn Þand IðR Þ denote the sets of all closed real interval numbers, n-dimensional real interval vectors and n n real interval matrices, respecx 2 IðRÞ can be usually written in the following tively. X I ¼ ½x; form: nn
X I ¼ ½X c DX; X c þ DX
ðA1Þ
c
I
in which X and DX denote the mean (or midpoint) value of X and the uncertainty (or the maximum width) in X I , respectively. It follows that
x þ x 2 x x DX ¼ 2 Xc ¼
ðA2Þ ðA3Þ
In terms of the interval addition, Eq. (A1) can be put into the more useful form
X I ¼ X c þ DX I
ðA4Þ
I
where DX ¼ ½DX; DX . A n-dimensional real interval vector XI 2 IðRn Þ can be expressed as
XI ¼ ðX I1 ; X I2 ; . . . ; X In ÞT
ðA5Þ
The mean value and uncertainty of XI are
T Xc ¼ X c1 ; X c2 ; . . . ; X cn
ðA6Þ
DX ¼ ðDX 1 ; DX 2 ; . . . ; DX n ÞT
ðA7Þ
Similar expression 2 IðRnn Þ AI ¼ A; A
exists
for
AI ¼ Ac þ DAI I
an
nn
interval
matrix
ðA8Þ c
where DA ¼ ½DA; DA; A and DA denote the mean matrix of AI and the uncertain (or the maximum width) matrix of AI , respectively. It follows that
ij þ aij AþA a or acij ¼ 2 2 ij aij AA a DA ¼ or Daij ¼ 2 2 Ac ¼
ðA9Þ ðA10Þ
where Ac ¼ acij and DA ¼ Daij . Appendix B. Natural extension of the real rational function An interval function is an interval-valued function of one or more interval arguments. Assuming that F XI ¼ F X I1 ; X I2 ; . . . ; X In is the interval value function of interval variable XI ¼ X I1 ; X I2 ; . . . ; X In Þ, if X Ii # Y Ii ; i ¼ 1; 2; . . . ; n, one has
F X I1 ; X I2 ; . . . ; X In # F Y I1 ; Y I2 ; . . . ; Y In
ðA11Þ
We say that the interval value function FðXI Þ ¼ F X I1 ; X I2 ; . . . ; X In of the interval XI ¼ X I1 ; X I2 ; . . . ; X In is inclusion monotonic, if f is the real function of nreal variables x1 ; x2 ; . . . ; xn and the interval value function F of n interval variables X I1 ; X I2 ; . . . ; X In satisfy
Fðx1 ; x2 ; . . . ; xn Þ ¼ f ðx1 ; x2 ; . . . ; xn Þ;
xi 2 X Ii ; i ¼ 1; 2; . . . ; n
ðA12Þ
F is known as the interval extension of f. Real rational functions of n real variables may have natural extensions. Given a rational expression in real variables, we can replace the real variables by corresponding interval variables and replace the real arithmetic operations by the corresponding interval arithmetic operations to obtain a rational interval function, which is called natural extension of the real rational function. The extensions of the real rational function are inclusion monotonic and they can be calculated through finite-interval arithmetic operations. References [1] Nagpal VK et al. Probabilistic structural analysis to quantify uncertainties associated with turbopump blades. AIAA J 1989;27:809–13. [2] Hien TD, Kleiber M. Finite element analysis based on stochastic Hamilton variation principle. Comput Struct 1990;37:893–902. [3] Bellomo N, Casciati F. Nonlinear stochastic mechanics. Berlin: Springer; 1992. [4] Deif A. Sensitivity analysis in linear systems. Berlin: Springer; 1986. [5] Ben-Haim Y, Elishakoff I. Convex models of uncertainty in applied mechanics. New York: Elsevier; 1990. [6] Lindberg HE. Dynamic response and buckling failure measures for structures with bounded and random imperfections. Trans ASME J Appl Mech 1991;58:1092–4. [7] Liu ZS, Chen SH, Han WZ. Solving the extremum of static response for structural systems with uncertain unknown-but-bound parameters. Comput Struct 1994;50:557–61. [8] Tzan SR, Pantelides CP. Convex models for impulsive response of structures. J Eng Mech 1996;122:521–9. [9] Pantelides CP. Buckling and postbuckling of stiffened elements with uncertainty. Thin-Wall Struct 1996;26:1–17. [10] Ganzerli S, Pantelides CP. Optimum structural design via convex model superposition. Comput Struct 2000;74:639–47. [11] Pantelides CP, Booth BC. Computer-aided design of optimal structures with uncertainty. Comput Struct 2000;74:293–307. [12] Moore RE. Interval analysis. Englewood Cliffs, New Jersey: Prentice-Hall; 1966. [13] Alefeld G, Herzberger J. Introduction to interval computations. New York: Academic Press; 1983. [14] Chen SH, Zhang XM, Chen YD. Interval eigenvalues of closed-loop systems of uncertain structures. Comput Struct 2006;84:243–53. [15] Chen SH, Wu XM, Yang ZJ. Eigensolution reanalysis of modified structures with epsilon-algorithm. Int J Numer Method Eng 2006;66:2115–30. [16] Graves-Morris PR, Roberts DE, Salam A. The epsilon-algorithm and related topics. J Comput Appl Math 2000;122:51–80. [17] Brezinski C. Convergence acceleration of some sequence by epsilon-algorithm. Numer Math 1978;29:173–7. [18] Sidi A. Extension and completion of Wynn’s theory on convergence of columns of the epsilon table. J Approx Theory 1996;86:21–40. [19] Roberts DE. On the convergence of rows of vector Pade’ approximation. J Comput Appl Math 1996;70:95–109. [20] Rao SS. Optimization theory and applications. 2nd ed. Wiley Eastern; 1984. [21] Yang XW, Chen SH, Lian HD. Interval finite element method based on the element for eigenvalue analysis of structures with interval parameters. Struct Eng Mech 2001;12:669–84.
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
An efficient method for evaluating the natural frequencies of structures with uncertain-but-bounded parameters Su Huan Chen *, Liang Ma, Guang Wei Meng, Rui Guo Department of Mechanics, Nanling Campus, Jilin University, Changchun 130025, PR China
a r t i c l e
i n f o
Article history: Received 30 August 2008 Accepted 24 February 2009 Available online 2 April 2009 Keywords: Efficient method Natural frequencies Uncertain-but-bounded parameters Epsilon-algorithm
a b s t r a c t Using interval theory and the second-order Taylor series, the eigenvalue problems of structures with multi-parameter can be transformed into those with single parameter. The epsilon-algorithm is used to accelerate the convergence of the Neumann series to obtain the bounds of eigenvalues of structures with single interval parameter, thus increasing the computing accuracy and reducing the computational effort. Finally, the effect of uncertain parameters on natural frequencies is evaluated. Two engineering examples show that the proposed method can give better results than those obtained by the first-order approximation, even if the uncertainties of parameters are fairly large. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction The natural frequency analysis for structures with deterministic parameters has been extensively developed. However, in engineering situations, the structural parameters are often uncertain, such as the inaccuracy of the measurement, errors in the manufacturing and assembly process, invalidity of some components and uncertainty in boundary conditions etc. The uncertainties of structural parameters may lead to large and unexpected excursion of responses that may lead to drastic reduction in accuracy and precision of the operation. Therefore, uncertainty plays an important role in the modern engineering structural analysis. Over the past decades, a number of methods have been developed that include uncertain model properties in the finite element (FE) analysis and aim at the quantification of the uncertainty on the analysis result. The probabilistic concept is by far the most popular method for numerical uncertainty modeling. Its popularity has led to a number of probabilistic FE procedures [1–3]. However, probabilistic modeling is not the only way to describe the uncertainty, and uncertainty is not equal to randomness. Indeed, the probabilistic approaches are not able to deliver reliable results at the required precision without sufficient experimental data to validate regarding the joint probability densities of the random variables or functions involved. Therefore, one may recognize that uncertainties in parameters can be modeled on the basis of alterative, non-probabilistic conceptual frameworks. One such approach, based on a set theoretic formulation, is an unknown-but-bounded model (convex models). Such set models of uncertainty have been * Corresponding author. Tel.: +86 431 8509 5259; fax: +86 431 8509 5288. E-mail address: [email protected] (S.H. Chen). 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.02.009
applied by Deif in linear programming and system theory [4]. Recently, such set models of uncertainties in parameters have drawn interest both from the system control robustness analysis field and from the structural failure measures field. For example, the convex model was introduced by Ben-Haim and Elishakoff [5], discussed later by Lindberg [6], for the study of dynamic response and failure of structures under pulsed parametric loading; the convex model has been applied in determining the upper and lower bounds of static response for structures by Liu et al. [7]. The convex model has also been applied to impulsive response, buckling analysis and optimal design of structures with uncertain parameters [8–11]. Since the mid-1960s, a new method called the interval analysis has appeared. Moore [12], Alefeld and Herzberger [13] have done the pioneering work. Mathematically, linear interval equations, nonlinear interval equation, and interval eigenvalue problems have been partially resolved. But because of the complexity of the algorithm, it is difficult to apply these results to practical engineering problems. Recently, the interval finite element (IFE) method was presented by Chen et al. [14] which makes the method easier to deal with the interval eigenvalues for closed-loop systems of structures with uncertain parameters. It should be pointed out that the previous interval analysis methods for computing the upper and lower bounds of response of structures with interval parameters are based on the first-order perturbation or the first-order Taylor series. Although the procedure is easy to implement on the computer and incorporates the finite element code, the applications of the previous methods are limited to the case when the interval uncertainties of parameters are small. However, if the parameter uncertainties are fairly large or the combination of a large number uncertain parameters, the
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
accuracy of the computational results will become unacceptable. Thus, it is highly desirable to present a more accurate method for computing the upper and lower bounds of responses of structures with fairly large uncertainties of interval parameters. To this end, this paper presents an efficient method to estimate the natural frequencies of structures for the case with fairly large uncertainties of parameters. As we known, the eigenvalues k is defined as x2 in mode analysis of FEM problems, where x denotes the natural frequencies of structures. The idea of the proposed method is that the eigenvalues are considered as functions of the structural parameters; using interval theory and the second-order Taylor series expansion, we transform the eigenvalue problems of structures with multi-parameter into those with single parameter; The epsilon-algorithm is used to accelerate the convergence of Neumann series to compute the eigenvalues of structures with single interval parameter, thus the interval eigenvalues of the structures with multi-parameter are obtained and the effect of uncertain parameters on natural frequencies is evaluated finally. This paper is organized as follows. A brief review of mathematical background on the interval analysis is given in the Appendix. In the Section 2, the epsilon-algorithm for eigensolution reanalysis of structures with fairly large changes of parameters is discussed [15]. The Section 3 presents the definition of the interval eigenvalue problems of the interval parameter structures. In the Section 4, the proposed method for evaluating the natural frequencies of structures with fairly large uncertainties of parameters is developed. In Section 5, two engineering examples are given to illustrate the application of the proposed method. The results obtained by the proposed method are compared with those obtained by the exact solutions and the first-order approximation. The conclusions are drawn in Section 6.
where K0 and M0 are the stiffness and mass matrices of the finite element assemblage. This equation will be referred to as the initial problem in the following discussions. The perturbation method studies the changes of eigenvalues of the system subjected to changes in its design parameters. Therefore, if the initial system is represented by Eq. (5), the problem becomes that to determine u and k when K0 and M0 are perturbed to the form K0 þ DK and M0 þ DM, respectively.
Ku ¼ kMu
ð6Þ
where
K ¼ K0 þ DK
ð7Þ
M ¼ M0 þ DM
ð8Þ
in Eqs. (6)–(8), if we introduce the following notations
f 0 ¼ k0 M0 u0
ð9Þ
f ¼ kMu
ð10Þ
Df ¼ kMu k0 M0 u0
ð11Þ
then we obtain
ðK0 þ DKÞu ¼ f 0 þ Df:
u ¼ ðK0 þ DKÞ1 ðf 0 þ DfÞ 1 1 ¼ ðI þ K1 0 DKÞ K0 ðk0 M0 u0 þ k0 DMu0 Þ
¼ ðI þ BÞ1 K1 0 ðk0 M0 u0 þ k0 DMu0 Þ
2.1. The epsilon-algorithm
then obtaining a series
ðjÞ kþ1
e
ð1Þ
ðjþ1Þ k1
¼e
þ
h
ðjþ1Þ k
e
ðjÞ k
e
i1
ð2Þ j; k ¼ 0; 1; 2; . . .
ð3Þ
The iterative formulae (1)–(3) are similar to the scalar case except that it requires the inverse of a vector. The definition of the inverse of the vector is given by Roberts [19]
u1 ¼
u u ¼ Pn H 2 ðu uÞ i¼1 jui j
ð4Þ
ð13Þ
In Eq. (13), B ¼ K1 0 DK; k and u were approximated by k0 and u0 , respectively. By using the Neumann series expansion, we have
~ ðI B þ B2 ÞK1 u 0 ðk0 M0 u0 þ k0 DMu0 Þ
eðjÞ 1 ¼ 0 ðjÞ e0 ¼ sj
ð12Þ
It follows that
2. Eigensolution reanalysis with the epsilon-algorithm
The epsilon-algorithm was presented [16–18] to accelerate the convergence of a sequence. Given a vector sequence fs0 ; s1 ; s2 ; . . .g, we construct the iterative form to obtain the vector sequence in the epsilon-algorithm as follows:
583
ð14Þ
~ 0 ¼ K1 u 0 ðk0 M0 u0 þ k0 DMu0 Þ 1 1 ~ 1 ¼ K1 ~ u 0 DKK0 ðk0 M0 u0 þ k0 DMu0 Þ ¼ K0 DKu0
~ ~ 2 ¼ K1 u 0 DKu1 .. . e ~ s ¼ K1 u 0 DK u s1
ð15Þ
s ¼ 3; 4; . . .
Assume the solution of the Eq. (6) has the following form
~1 þ þ u ~s þ ~¼u ~0 þ u u
ð16Þ
~ 0 ; s1 ¼ We define a vector sequence fs0 ; s1 ; . . . ; ss ; . . .g where s0 ¼ u ~ 1 ; s2 ¼ u ~0 þ u ~1 þ u ~ 2 , and in general ~0 þ u u
si ¼
i X
~j u
i ¼ 0; 1; 2; . . . ; s; . . .
ð17Þ
j¼0
where the u denotes the conjugate of u and uH the conjugate and transpose of u. It should be noted that for the undamped structure the modal vectors are real. Thus, u in Eq. (4) should be replaced by u, and uH should be replaced by uT . In this case, u1 in Eq. (4) is just a normalized version of u itself. The vector epsilon-algorithm table can be constructed by Eqs. (1)–(3).
To accelerate the convergence of the sequence (17), using the vector epsilon-algorithm (Eqs. (1)–(3)), the sequence, fs0 ; s1 ; . . . ; ss ; . . .g, can be obtained. Refs. [16–18] have pointed out that in ðjÞ the epsilon iterative formulae, when i is odd, the ei are meaningðjÞ less, and that when the i ¼ 2k ðk ¼ 1; 2; . . .Þ, the vector e2k is the Shanks transformation. The solution is the last even row in the epsilon-algorithm,
2.2. The epsilon-algorithm for accelerating the convergence of Neumann series
u ¼ e2k
ðjÞ
ð18Þ
Compute the eigenvalues using the Rayleigh quotients The eigenproblem for the structure is as follows
K0 u0 ¼ k0 M0 u0
ð5Þ
ki ¼
uTi Kui uTi Mui
ð19Þ
584
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
3. Generalized interval eigenvalue problems 3.1. The definition of generalized interval eigenvalue problems The eigenvalue problem is described by equation
KðbÞu ¼ kMðbÞu
ð20Þ
under the parameter constraint
b6b6b
ð21Þ
Using the interval matrix expression and the interval extension, Eq. (20) can be expressed as I
I
Kðb ÞuI ¼ kI Mðb ÞuI
ð22Þ
Eq. (22) is called interval eigenvalue problem, under constraint (21). I The basic problem is, for given interval parameter vector b , to I find the interval eigenvalue ki , which is the smallest interval and encloses all possible eigenvalue ki , satisfying KðbÞu ¼ kMðbÞu. In other words, we seek a hull, i.e.
n
C ¼ ki : ki 2 R; KðbÞu ¼ ki MðbÞu; ki –0; b 2 bI
o
ð23Þ
to the set
kIi ¼ ki ; ki
ð24Þ
where I
b2b
ki ¼ max ki ðhKðbÞ; MðbÞiÞ ¼ max ki ðbÞ;
ð25Þ I
b2b
ð26Þ
where ki ðhKðbÞ; MðbÞiÞ means solving the eigenvalues for the generalized eigenvalue problems with KðbÞ and MðbÞ described by Eq. (20). 3.2. The interval eigenvalues of structures Using Eq. (A1), interval parameter vector can be written as
h i I c c b ¼ b; b ¼ b Db; b þ Db
ð27Þ
where
Db ¼ ðDb1 ; Db2 ; . . . ; Dbm ÞT
ð28Þ
Eqs. (27) and (28) mean that the uncertain parameters are constrained into a m-dimension rectangle with widths 2Db1 ; 2Db2 ; . . . ; 2Dbm . For a m-dimension rectangle, we can define its boundary points P bd which are placed on the boundary surfaces and the extreme points P ex which are placed at the intersections of the boundary surfaces. For example, we consider a 2-dimension rectangle shown in Fig. 1. The boundary points P bd are placed on P1ex P2ex ; P1ex P 3ex ; P3ex P4ex and P2ex P4ex , and the extreme points are the intersections P 1ex ; P 2ex ; P 3ex and P 4ex . It means that the extreme points Psex can be determined by the parameter vector
¼ Pex ðb1 ; b2 ; . . . ; bm Þ : bk ¼
k ¼ 1; 2; . . . ; m s ¼ 1; 2; . . . ; 2
on the surface’s boundary line and the extrema of any boundary line will occur on the line’s extreme points. According to the definitions of the extreme points Pex and the boundary points P bd , the extrema of Eqs. (25) and (26) will occur on the extreme points Psex ; s ¼ 1; 2; . . . ; 2m , thus we obtain the upper and lower bounds of the interval eigenvalues.
ki ¼ min ki ðPsex Þ;
i ¼ 1; 2; . . . ; n s ¼ 1; 2; . . . 2m
ð30Þ
max ki ðPsex Þ;
i ¼ 1; 2; . . . ; n s ¼ 1; 2; . . . 2m
ð31Þ
ki ¼
ki ¼ min ki ðhKðbÞ; MðbÞiÞ ¼ min ki ðbÞ;
Psex
Fig. 1. 2-Dimension rectangle.
c bk m
þ Dbk or
c bk
Dbk
From Eqs. (30) and (31), it can be seen that in order to obtain the exact solutions of the upper and lower bounds of eigenvalues, we have to solve 2m eigenvalue problems given by Eq. (20) on the extreme points P sex ðs ¼ 1; 2; . . . ; 2m Þ of the m-dimension rectangle, then take min ki ðPsex Þ as the lower bounds of eigenvalues ki ði ¼ 1; 2; . . . ; nÞ, and take max ki ðPsex Þ as the upper bounds of eigenvalues ki ði ¼ 1; 2; . . . ; nÞ. This is a tedious job which would not be acceptable if m is too large. Therefore, we need to seek an efficient method for computing the upper and lower bounds of eigenvalues of structures with interval parameters. 4. Methods for evaluating the natural frequencies of engineering structures with uncertain parameters As we known, the eigenvalues k is defined as x2 in mode analysis of FEM problems, where x denotes the natural frequencies of structures. If the bounds of eigenvalues are obtained, the natural frequencies of structures are evaluated correspondingly. 4.1. The method based on the first-order approximation Using the first-order Taylor series and the first-order perturbation for eigenvalue analysis, Ref. [21] developed a method for computing the bounds of eigenvalues as follows
c I c k b ¼ k b Dk; k b þ Dk
ð32Þ
ð29Þ m
where m is the number of the interval parameters, 2 denotes the number of the extreme points. Eqs. (25) and (26) call for finding the minimum and maximum I of the function ki ðbÞ in the interval set b . According to optimal theory [20], the extrema of Eqs. (25) and (26) will occur on the boundary surfaces of the m-dimension rectangle described by Eqs. (27) and (28). Similarly, the extrema of any boundary surface will occur
where
Dk ¼
m
X c
c T @KðbÞ
u b
c Db i u b
@b i b ¼b i¼1 i i
c c T @MðbÞ
c
k b u b Db i u b
@bi bi ¼bc
ð33Þ
i
c
I
in which bi and Dbi denote the mean (or midpoint) value of bi and I the uncertainty (or the maximum width) in bi , respectively.
585
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
4.2. The proposed method for evaluating the natural frequencies of structures with fairly large uncertainties of the parameters Using the second-order Taylor series to expand ki ðb1 ; b2 ; . . . ; bm Þ c c c c about the point b ¼ b1 ; b2 ; . . . ; bm , we have
c c c 1 c c T c ki ðbÞ ¼ ki b þ GTki b b b þ b b Hki b b b 2
ð34Þ
c c where b and Hki b are the gradient vector and the Hessian c matrix of eigenvalues at the point b , respectively. The gradient vector and the Hessian matrix can be written as follows GTki
GTki ðbÞ ¼
@ki @ki @ki @b1 @b2 @bm
2
@ 2 ki @b1 @b1
@ 2 ki @b1 @b2
6 6 @ 2 ki 6 Hki ðbÞ ¼ 6 @b2 @b1 6 4
@ 2 ki @b2 @b2
@ 2 ki @bm @b1
@ 2 ki @bm @b2
ð35Þ @ 2 ki @b1 @bm
3
7 7 7 7 7 5
@ 2 ki @b2 @bm
ð36Þ
@ 2 ki @bm @bm
c @ki
c c c
ki b1 ; b2 ; bm ¼ ki b þ b1 b1
@b1 b1 ¼bc 1
2 1 @ ki
c 2 þ b1 b1
2 @b1 @b1
c
2
I I I I c c c I c c c I ki b1 ; b2 ; b3 ki b1 ; b2 ; b3 þ ki b1 ; b2 ; b3 þ ki b1 ; b2 ; b3 c c c 2ki b1 ; b2 ; b3
ð42Þ
As we have discussed in the Eqs. (30) and (31), exact solutions I I I for ki b1 ; b2 ; b3 can be obtained by solving 23 finite element prob-
I I I ki b1 ; b2 ; b3 " " # ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ¼ min ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ " ## ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ð43Þ Similarly, we have that for 1-dimension case
1 @ ki
c ðb2 b2 Þ þ
2 @b2 @b2
b2 ¼bc 2
2
c 2 b2 b2
b2 ¼bc2
ð38Þ Substituting Eqs. (37) and (38) into Eq. (34), we have
c c c c ki ðbÞ ¼ ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b c 1 c T c þ b b HLki b b b 2
ð41Þ
From Eq. (41), it can be seen that the computational problems of the interval eigenvalues of structures with multi-parameter have been transformed into those with single parameter. The equation can be written as follows while m = 3
ð37Þ
b1 ¼b1
c @ki
c c ki b1 ;b2 ; .. .bm ¼ ki ðb Þ þ @b
I I c c c I c ki b ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c I þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b
lems, that is
The computation of the Hessian matrix is a tedious job if m is c large. For the sake of simplicity, we consider b2 ¼ b2 ; b3 ¼ c c b3 ; . . . ; bm ¼ bm in Eq. (34), so ki ðbÞ becomes a function of b1 . We have
Similarly, we obtain
Because of the complexities of solving mixed partial derivatives and the widening of computed intervals caused by the multiplications and accumulations of different interval parameters, we neglect the last part of the right side of Eq. (40) and have
ð39Þ
h h i I c c c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð44Þ
h h i c I c c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð45Þ
h h i c c I c c c c ki ðb1 ; b2 ; b3 Þ ¼ min ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ ; h ii c c c c max ki ðb1 ; b2 ; b3 Þ; ki ðb1 ; b2 ; b3 Þ
ð46Þ
Substituting Eqs. (44)–(46) into Eq. (42), we have
h i h i 3 c c c c c c c c min ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 þ min ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 7 6 h 7 6 i c c c c c c c 7 6 ; k 2k ; þ min k b ; b ; b b ; b ; b b ; b ; b 3 3 i i i 7 6 1 2 1 2 1 2 3 I I I h i h i 7 ki b1 ; b2 ; b3 6 7 6 c c c c c c c c 6 max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 þ max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 7 7 6 5 4 h i c c c c c c c þ max ki b1 ; b2 ; b3 ; ki b1 ; b2 ; b3 2ki b1 ; b2 ; b3 2
extension of the Using the natural interval function, yields I I c c c I c ki b ¼ ki b1 ; b2 ; . . . ; bm þ ki b1 ; b2 ; . . . ; bm þ c c c I þ ki b1 ; b2 ; . . . ; bm ðm 1Þki b T 1 I c c I c þ b b HLki b b b ð40Þ 2
c where HLki b denotes the Hessian matrix in which the diagonal elements are zero.
ð47Þ
From Eq. (47), it can be seen that, using the proposed method, the problem for solving 23 FE equations has been transformed into that for solving 2 3 FE equations. In order to further reduce the computation effort, we suggest that based on the Neumann series expansions, the accurate approximations of solutions for the upper and lower bounds of eigenvalues can be obtained by using the reanalysis method with the epsilon-algorithm given by the Section 2 for structures with a single interval parameter.
586
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
In general, the exact solution for the interval eigenvalues can be obtained by solving 2m FE problems, and using proposed method the approximate one can be obtained by 2m reanalysis problems, thus reducing the computational effort. Compared with the first-order approximation, the proposed method replaces the first-order perturbation with the reanalysis method based on epsilon-algorithm which considers many higher order items to make the result more accurate. Although we neglect the mixed second-order items of Taylor series in Eq. (40), both the first-order and the second-order items of each parameter are included. However, only the first-order items of Taylor series are considered in the first-order approximation. The complete method presented in this paper is summarized below: 1. Transform the interval eigenvalues of structures with multiparameter into those of structures with single parameter using Eq. (41); 2. Compute the upper and lower bounds of eigenvalues of structures with a single interval parameter using the reanalysis method given by Eqs. (15) and (1)–(3); 3. Using Eq. (41), the interval eigenvalues of structures with multi-parameter are obtained so that the natural frequencies are evaluated finally.
5. Numerical examples 5.1. Example 1 We consider a 40-story structure as shown in Fig. 2. The finite element model of the given structure consists of 202 nodes, 357 beam elements and 1212 degrees of freedom. The parameters of the structure are as follows. The length of the longitudinal beam is 5.0 m and that of cross beam is 3.0 m; The Young’s modulus of the used material is E ¼ 2:1 1011 Pa; the mass density q ¼ 7:8 103 kg=m3 ; the height and width of the cross-section of the beams are H ¼ 0:1 m, and W ¼ 0:2 m; Assume that E; H and W are the interval parameters given by
EI ¼ Ec þ aEc e HI ¼ Hc þ bHc e W I ¼ W c þ cW c e e ¼ ½1; 1 The results are listed by frequencies. Case 1 (one parameter). Consider the case with one parameter with uncertainty (a ¼ 5%). The upper and lower bounds of the frequencies of the structure obtained by the firstorder approximation are listed in Table 1. The corresponding results obtained by the proposed method are listed in Table 2. For comparison, the exact solutions of the bounds of frequencies obtained by Eq. (43) are also listed in Tables 1 and 2. From results in Tables 1 and 2, it can be seen that the max error of the upper bound of the mode 2 for the first-order approximation is 1.3%, while the corresponding error is 0.0% for the proposed method. Case 2 (two parameters). Consider the case with two parameters with uncertainties (a ¼ b ¼ 5%Þ. The results are listed in Tables 3 and 4, and show that the max error obtained by the first-order approximation of the lower bound of the mode 5 reaches 4.4%, while the corresponding error obtained by the proposed method is 0.0%. Case 3 (three parameters). Consider the case with three parameters with uncertainties (a ¼ b ¼ c ¼ 5%Þ. The results are listed in Tables 5 and 6, and show that the max error obtained by the first-order approximation for the lower bound of the mode 5 reaches 4.4%, while the corresponding error obtained by the proposed method is 0.0%. The max error of the lower bound of the mode 3 is 0.6%.If the uncertainties are a ¼ b ¼ c ¼ 10%,
Fig. 2. 40-Story frame structure.
Table 1 Bounds of frequencies using the first-order approximation (a ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.29D03, [2.11D02, [4.98D02, [6.00D02, [1.19D01,
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
1.2 0.9 1.2 1.0 0.8
0.9 1.3 0.9 0.9 0.8
3.53D03] 2.27D02] 5.35D02] 6.43D02] 1.27D01]
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
587
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590 Table 2 Bounds of frequencies using the proposed method (a ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
[3.33D03, [2.13D02, [5.04D02, [6.06D02, [1.20D01,
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
3.50D03] 2.24D02] 5.30D02] 6.37D02] 1.26D01]
Table 3 Bounds of frequencies using the first-order approximation (a ¼ b ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.03D03, [1.95D02, [4.68D02, [5.53D02, [1.09D01,
[3.16D03, [2.03D02, [4.84D02, [5.76D02, [1.14D01,
4.1 3.9 3.3 4.0 4.4
2.2 2.1 2.2 2.2 2.3
3.75D03] 2.41D02] 5.62D02] 6.84D02] 1.35D01]
3.67D03] 2.36D02] 5.50D02] 6.69D02] 1.32D01]
Table 4 Bounds of frequencies using the proposed method (a ¼ b ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.15D03, [2.02D02, [4.82D02, [5.74D02, [1.14D01,
[3.16D03, [2.03D02, [4.84D02, [5.76D02, [1.14D01,
0.3 0.5 0.4 0.3 0.0
0.3 0.4 0.2 0.1 0.0
3.66D03] 2.35D02] 5.49D02] 6.68D02] 1.32D01]
3.67D03] 2.36D02] 5.50D02] 6.69D02] 1.32D01]
Table 5 Bounds of frequencies using the first-order approximation (a ¼ b ¼ c ¼ 5%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.03D03, [1.95D02, [4.59D02, [5.53D02, [1.09D01,
[3.16D03, [2.03D02, [4.79D02, [5.76D02, [1.14D01,
4.1 3.9 4.2 4.0 4.4
2.5 2.1 2.3 2.2 2.3
3.76D03] 2.41D02] 5.69D02] 6.84D02] 1.35D01]
3.67D03] 2.36D02] 5.56D02] 6.69D02] 1.32D01]
Table 6 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 5%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[3.15D03, [2.02D02, [4.76D02, [5.74D02, [1.14D01,
[3.16D03, [2.03D02, [4.79D02, [5.76D02, [1.14D01,
0.3 0.5 0.6 0.3 0.0
0.3 0.4 0.4 0.1 0.0
3.66D03] 2.35D02] 5.54D02] 6.68D02] 1.32D01]
3.67D03] 2.36D02] 5.56D02] 6.69D02] 1.32D01]
Table 7 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 10%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
3.41D03 2.19D02 5.17D02 6.22D02 1.23D01
[2.88D03, [1.85D02, [4.32D02, [5.24D02, [1.04D01,
[2.91D03, [1.87D02, [4.41D02, [5.31D02, [1.05D01,
1.0 1.1 2.0 1.3 1.0
0.8 0.8 1.2 0.8 0.7
3.91D03] 2.51D02] 5.89D02] 7.12D02] 1.41D01]
3.94D03] 2.53D02] 5.96D02] 7.18D02] 1.42D01]
588
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
The results are listed by frequencies.
from Table 7 it can be seen that the max error of the lower bound of the mode 3 obtained by the proposed method is only 2.0%.
Case 1 (one parameter). Consider the one parameter case with uncertainty (a ¼ 3%Þ. The upper and lower bounds of the frequencies of the structure obtained by the first-order approximation are listed in Table 8. The corresponding results obtained by the proposed method are listed in Table 9. For comparison, the exact solutions of the bounds of frequencies obtained by Eq. (43) are also listed in Tables 8 and 9. From results in Tables 8 and 9, it can be seen that the max error of the upper bound of the mode 2 for the first-order approximation is 10.5%, while the corresponding error is 0.0% for the proposed method. Case 2 (two parameters). Consider the case with two parameters with uncertainties (a ¼ b ¼ 3%Þ. The
5.2. Example 2 We consider a chassis structure as shown in Fig. 3. The finite element model of the given structure consists of 20,953 nodes, 19,880 thin shell elements and 125,718 degrees of freedom. The parameters of the structure are as follows. The length of the chassis is 3.0 m, the front width of the chassis 0.9 m, the back width 1.2 m; The Young’s modulus of the material used is E ¼ 2:1 1011 Pa; the mass density q ¼ 7:8 103 kg=m3 ; the thickness of the thin shell in the longitudinal beam is t1 ¼ 0:01 m, and that of the cross beam is t2 ¼ 0:01 m. Assume that E; t 1 and t 2 are the interval parameters given by
EI ¼ Ec þ aEc e t I1 ¼ t c1 þ bt c1 e tI2 ¼ t c2 þ ct c2 e e ¼ ½1; 1
Fig. 3. Chassis structure.
Table 8 Bounds of frequencies using the first-order approximation (a ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.26D+00, [1.88D+01, [3.38D+01, [3.39D+01, [4.20D+01,
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
1.3 10.5 0.6 6.4 5.8
1.3 8.8 0.6 5.6 5.2
8.73D+00] 2.36D+01] 3.52D+01] 3.94D+01] 4.83D+01]
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
Table 9 Bounds of frequencies using the proposed method (a ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
[8.37D+00, [2.10D+01, [3.40D+01, [3.62D+01, [4.46D+01,
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
8.62D+00] 2.17D+01] 3.50D+01] 3.73D+01] 4.59D+01]
Table 10 Bounds of frequencies using the first-order approximation (a ¼ b ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.03D+00, [1.75D+01, [3.37D+01, [3.36D+01, [4.11D+01,
[8.25D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
2.7 16.3 0.6 7.2 7.4
2.2 12.8 0.6 6.1 6.3
8.94D+00] 2.46D+01] 3.53D+01] 3.97D+01] 4.90D+01]
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
589
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590 Table 11 Bounds of frequencies using the proposed method (a ¼ b ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.24D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
[8.25D+00, [2.09D+01, [3.39D+01, [3.62D+01, [4.44D+01,
0.1 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
8.75D+00] 2.18D+01] 3.51D+01] 3.74D+01] 4.61D+01]
Table 12 Bounds of frequencies using the first-order approximation (a ¼ b ¼ c ¼ 3%). Mode
Initial frequency (Hz)
The first-order approximate (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[7.92D+00, [1.39D+01, [3.32D+01, [3.31D+01, [3.99D+01,
[8.20D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
3.4 32.2 1.2 8.3 9.5
2.7 20.7 0.8 7.2 8.0
9.04D+00] 2.68D+01] 3.58D+01] 4.02D+01] 5.00D+01]
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
Table 13 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 3%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[8.19D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
[8.20D+00, [2.05D+01, [3.36D+01, [3.61D+01, [4.41D+01,
0.1 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
8.80D+00] 2.22D+01] 3.55D+01] 3.75D+01] 4.63D+01]
Table 14 Bounds of frequencies using the proposed method (a ¼ b ¼ c ¼ 10%). Mode
Initial frequency (Hz)
Approximate solution obtained by the proposed method (Hz)
Exact solution (Hz)
Error of lower bound (%)
Error of upper bound (%)
1 2 3 4 5
8.50D+00 2.13D+01 3.45D+01 3.68D+01 4.53D+01
[7.46D+00, [1.84D+01, [3.13D+01, [3.43D+01, [4.14D+01,
[7.51D+00, [1.87D+01, [3.14D+01, [3.43D+01, [4.15D+01,
0.7 1.6 0.3 0.0 0.2
0.4 0.4 0.5 0.0 0.2
9.49D+00] 2.40D+01] 3.75D+01] 3.90D+01] 4.87D+01]
results are listed in Tables 10 and 11, and show that the max error obtained by the first-order approximation of the lower bound of the mode 2 reaches 16.3%, while the corresponding error obtained by the proposed method is 0.0%. Case 3 (three parameters). Consider the case with three parameters with uncertainties (a ¼ b ¼ c ¼ 3%Þ. The results are listed in Tables 12 and 13, and show that the max error obtained by the first-order approximation for the lower bound of mode 2 reaches 32.2%, while the corresponding error obtained by the proposed method is 0.0%. The max error of the lower bound of mode 1 is 0.1%.If the uncertainties of the parameters are a ¼ b ¼ c ¼ 10%, from Table 14 it can be seen that the max error of the lower bound of the mode 2 obtained by the proposed method is only 1.6%.
9.53D+00] 2.41D+01] 3.77D+01] 3.90D+01] 4.88D+01]
The above numerical results indicate that the first-order approximation is valid for the case when the uncertainties of parameters are small. The proposed method can give good results even if the uncertainties of parameters are fairly large. The errors of the results obtained by the proposed method are attributed to neglecting of the mixed second-order items of Taylor series in Eq. (40). So the corresponding errors will go up as the uncertainties become very large, but they are much smaller than those obtained by the first-order approximation.
6. Conclusions An efficient method for estimating the natural frequencies of structures with fairly large uncertainties of parameters is developed. Based on interval theory and the second-order Taylor series, the solution for the eigenvalues have been reduced from 2m finite element analysis problems to 2m eigenvalue reanalysis ones. The epsilon-algorithm was used to accelerate the convergence of the Neumann series to compute the bounds of eigenvalues of structures with single interval parameter, thus increasing the computational accuracy and reducing the computing effort for the interval
590
S.H. Chen et al. / Computers and Structures 87 (2009) 582–590
eigenvalues of structures with multi-parameter. Two numerical examples were given to illustrate the application and validity of the proposed method. The three cases, one parameter, two parameters and three parameters were discussed. The comparison of numerical results obtained by the exact solution, the first-order approximation and the proposed method was given. The results indicated that the first-order approximation is valid only for the case when the uncertainties of parameters are small and the proposed method can give good results even if the uncertainties of parameters are fairly large. Acknowledgements This work is supported by Science and Technology Development Foundation of Jilin Province (20070541), 985 Engineering of Jilin Univ. and Innovation fund for 985 Engineering of Jilin Univ. No. 20080104. Appendix A. The definitions of the interval and interval operations Assuming that R is the set of all real numbers, IðRÞ; IðRn Þand IðR Þ denote the sets of all closed real interval numbers, n-dimensional real interval vectors and n n real interval matrices, respecx 2 IðRÞ can be usually written in the following tively. X I ¼ ½x; form: nn
X I ¼ ½X c DX; X c þ DX
ðA1Þ
c
I
in which X and DX denote the mean (or midpoint) value of X and the uncertainty (or the maximum width) in X I , respectively. It follows that
x þ x 2 x x DX ¼ 2 Xc ¼
ðA2Þ ðA3Þ
In terms of the interval addition, Eq. (A1) can be put into the more useful form
X I ¼ X c þ DX I
ðA4Þ
I
where DX ¼ ½DX; DX . A n-dimensional real interval vector XI 2 IðRn Þ can be expressed as
XI ¼ ðX I1 ; X I2 ; . . . ; X In ÞT
ðA5Þ
The mean value and uncertainty of XI are
T Xc ¼ X c1 ; X c2 ; . . . ; X cn
ðA6Þ
DX ¼ ðDX 1 ; DX 2 ; . . . ; DX n ÞT
ðA7Þ
Similar expression 2 IðRnn Þ AI ¼ A; A
exists
for
AI ¼ Ac þ DAI I
an
nn
interval
matrix
ðA8Þ c
where DA ¼ ½DA; DA; A and DA denote the mean matrix of AI and the uncertain (or the maximum width) matrix of AI , respectively. It follows that
ij þ aij AþA a or acij ¼ 2 2 ij aij AA a DA ¼ or Daij ¼ 2 2 Ac ¼
ðA9Þ ðA10Þ
where Ac ¼ acij and DA ¼ Daij . Appendix B. Natural extension of the real rational function An interval function is an interval-valued function of one or more interval arguments. Assuming that F XI ¼ F X I1 ; X I2 ; . . . ; X In is the interval value function of interval variable XI ¼ X I1 ; X I2 ; . . . ; X In Þ, if X Ii # Y Ii ; i ¼ 1; 2; . . . ; n, one has
F X I1 ; X I2 ; . . . ; X In # F Y I1 ; Y I2 ; . . . ; Y In
ðA11Þ
We say that the interval value function FðXI Þ ¼ F X I1 ; X I2 ; . . . ; X In of the interval XI ¼ X I1 ; X I2 ; . . . ; X In is inclusion monotonic, if f is the real function of nreal variables x1 ; x2 ; . . . ; xn and the interval value function F of n interval variables X I1 ; X I2 ; . . . ; X In satisfy
Fðx1 ; x2 ; . . . ; xn Þ ¼ f ðx1 ; x2 ; . . . ; xn Þ;
xi 2 X Ii ; i ¼ 1; 2; . . . ; n
ðA12Þ
F is known as the interval extension of f. Real rational functions of n real variables may have natural extensions. Given a rational expression in real variables, we can replace the real variables by corresponding interval variables and replace the real arithmetic operations by the corresponding interval arithmetic operations to obtain a rational interval function, which is called natural extension of the real rational function. The extensions of the real rational function are inclusion monotonic and they can be calculated through finite-interval arithmetic operations. References [1] Nagpal VK et al. Probabilistic structural analysis to quantify uncertainties associated with turbopump blades. AIAA J 1989;27:809–13. [2] Hien TD, Kleiber M. Finite element analysis based on stochastic Hamilton variation principle. Comput Struct 1990;37:893–902. [3] Bellomo N, Casciati F. Nonlinear stochastic mechanics. Berlin: Springer; 1992. [4] Deif A. Sensitivity analysis in linear systems. Berlin: Springer; 1986. [5] Ben-Haim Y, Elishakoff I. Convex models of uncertainty in applied mechanics. New York: Elsevier; 1990. [6] Lindberg HE. Dynamic response and buckling failure measures for structures with bounded and random imperfections. Trans ASME J Appl Mech 1991;58:1092–4. [7] Liu ZS, Chen SH, Han WZ. Solving the extremum of static response for structural systems with uncertain unknown-but-bound parameters. Comput Struct 1994;50:557–61. [8] Tzan SR, Pantelides CP. Convex models for impulsive response of structures. J Eng Mech 1996;122:521–9. [9] Pantelides CP. Buckling and postbuckling of stiffened elements with uncertainty. Thin-Wall Struct 1996;26:1–17. [10] Ganzerli S, Pantelides CP. Optimum structural design via convex model superposition. Comput Struct 2000;74:639–47. [11] Pantelides CP, Booth BC. Computer-aided design of optimal structures with uncertainty. Comput Struct 2000;74:293–307. [12] Moore RE. Interval analysis. Englewood Cliffs, New Jersey: Prentice-Hall; 1966. [13] Alefeld G, Herzberger J. Introduction to interval computations. New York: Academic Press; 1983. [14] Chen SH, Zhang XM, Chen YD. Interval eigenvalues of closed-loop systems of uncertain structures. Comput Struct 2006;84:243–53. [15] Chen SH, Wu XM, Yang ZJ. Eigensolution reanalysis of modified structures with epsilon-algorithm. Int J Numer Method Eng 2006;66:2115–30. [16] Graves-Morris PR, Roberts DE, Salam A. The epsilon-algorithm and related topics. J Comput Appl Math 2000;122:51–80. [17] Brezinski C. Convergence acceleration of some sequence by epsilon-algorithm. Numer Math 1978;29:173–7. [18] Sidi A. Extension and completion of Wynn’s theory on convergence of columns of the epsilon table. J Approx Theory 1996;86:21–40. [19] Roberts DE. On the convergence of rows of vector Pade’ approximation. J Comput Appl Math 1996;70:95–109. [20] Rao SS. Optimization theory and applications. 2nd ed. Wiley Eastern; 1984. [21] Yang XW, Chen SH, Lian HD. Interval finite element method based on the element for eigenvalue analysis of structures with interval parameters. Struct Eng Mech 2001;12:669–84.