- Email: [email protected]
An efficient matrix tridiagonalization method for 3D finite element analysis of free vibration
Journal Pre-proof An efficient matrix tridiagonalization method for 3D finite element analysis of free vibration A. Malatip, N. Prasomsuk, C. Siriparu, N. Paoprasert, S. Otarawanna
PII: DOI: Reference:
S0378-4754(19)30391-X https://doi.org/10.1016/j.matcom.2019.12.017 MATCOM 4921
To appear in:
Mathematics and Computers in Simulation
Received date : 16 February 2018 Revised date : 29 September 2019 Accepted date : 27 December 2019 Please cite this article as: A. Malatip, N. Prasomsuk, C. Siriparu et al., An efficient matrix tridiagonalization method for 3D finite element analysis of free vibration, Mathematics and Computers in Simulation (2020), doi: https://doi.org/10.1016/j.matcom.2019.12.017. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). ⃝ Published by Elsevier B.V. All rights reserved.
Journal Pre-proof
An Efficient Matrix Tridiagonalization Method for 3D Finite Element Analysis of Free Vibration
A. Malatip1,*, N. Prasomsuk2, C. Siriparu2, N. Paoprasert1 and S. Otarawanna1
1
Computer-Aided Engineering Research Team (CAET), Engineering Design and Advanced Manufacturing
Research Group, National Metal and Materials Technology Center (MTEC), National Science and Technology Development Agency (NSTDA), 114 Thailand Science Park, Thanon Pahonyothin, Tambon Khlong Nueng, Amphoe Khlong Luang, Pathum Thani 12120, Thailand
Department of Mechanical Engineering, Faculty of Engineering, Chulalongkorn University,
f
2
roo
Bangkok 10330, Thailand
rep
* Corresponding author, Email: [email protected]
Abstract
lP
In order to avoid mechanical resonance, a vibrating structure needs to be designed such that its working frequency interval is sufficiently far from its natural frequencies. Natural
rna
frequencies of a mechanical system are obtained from free vibration analysis commonly done by the finite element method. One important step in the analysis is matrix tridiagonalization
Jo u
commonly performed by the block Lanczos method. However, the classical block Lanczos method suffers numerical instability due to the loss of matrix orthogonality. In this paper, we demonstrate the implementation of the block Lanczos method with an orthogonality fixing scheme for 3D free vibration problems. The solution accuracy and computational time of this method were compared with those of the classical block Lanczos method and the Householder method. The results show that the block Lanczos method with the orthogonality fixing scheme
Journal Pre-proof employed in this work can effectively avoid the numerical instability due to the loss of matrix orthogonality. Furthermore, the block Lanczos method with the orthogonality fixing scheme provides solution accuracy as good as the Householder method while using significantly less computational time.
Keywords: matrix tridiagonalization; eigenvalue problem; free vibration; block Lanczos; orthogonality
1. Introduction One of the most common problems for a vibrating structure is mechanical resonance.
f
When the resonance occurs, the structure oscillates with higher amplitude than it does at other
roo
frequencies. A vibrating system is prone to mechanical resonance when the frequency of its oscillation matches one of the system's natural frequencies. Mechanical resonance can cause
rep
unwanted violent motions and even structural failure in structures including buildings, bridges, machines and vehicles. Therefore, design engineers must ensure that the working frequency interval of the mechanical system is sufficiently far from its natural frequencies.
lP
Natural frequencies of a mechanical system can be determined from free vibration analysis [Ginberg 2001; Bottega 2006]. For 3D complex geometries, this is commonly
rna
performed by the finite element method. In the finite element derivation of 3D free vibration problems, the mass and stiffness matrices can be transformed into a single standard-form
Jo u
matrix. If the size of the standard-form matrix is relatively small, eigenvalues and eigenvectors of the matrix can be directly solved within a relatively short period of time. For 3D free vibration problems, the size of the standard-form matrix is relatively large. Therefore, the matrix needs to be transformed into the tridiagonalized form (containing non-zero elements only on the main diagonal, and the first diagonal below and above the main diagonal) by a tridiagonalization algorithm.
Journal Pre-proof For the tridiagonalization operation, there are three renowned methods which are the Householder method [Bunse-Gerstner and Gragg 1988], the Lanczos method [Lanczos 1950, Parlett and Reid 1981, Cullum and Willoughby 1985] and the block Lanczos method [Golub and Underwood 1977]. The Householder method has a disadvantage when dealing with a sparse or structured matrix because the method destroys the matrix’s sparsity or structure [Guo and Qiao 2003]. On the contrary, the Lanczos method and the block Lanczos method can retain the matrix’s sparsity or structure when used as they involves only matrix-vector multiplication [Luk and Qiao 2003]. However, cancellation and roundoff errors make the Lanczos method and the block Lanczos method unstable due to the loss of orthogonality [Paige 1976, Iguchi et al. 1992].
roo
f
The loss of orthogonality in the Lanczos method can be fixed by different schemes, e.g. selective orthogonalization (SO) [Parlett and Scott 1979], partial orthogonalization (PRO) [Simon 1984] and reorthogonalization with improved convergence-check (RIC) [Imai et al.
rep
1991]. However, for relatively-large matrices, the block Lanczos method [Montgomery 1995] is preferred due to its better performance compared to the Lanczos method. This is because the
lP
block Lanczos method exploits computer memory hierarchies [Yu et al. 2003, Qiao et al. 2005]. To stabilize the block Lanczos algorithm, a componentwise technique for detecting and fixing
rna
the loss of orthogonality has been proposed and demonstrated on general eigenvalue problems [Qiao et al. 2005].
In this work, it aims to investigate the effectiveness of the block Lanczos method with
Jo u
this orthogonality fixing scheme for 3D free vibration problems. In the following sections, the implementation of this method is demonstrated and compared with the classical block Lanczos method and the Householder method [Householder 1958] in terms of solution accuracy and computational time.
Journal Pre-proof 2. Finite Element Derivation 2.1 Finite Element Equations The differential equation describing the motion in the Euler-Lagrange system is d L L − = 0. dt ri ri
(1)
The Lagrangian for a system of particles can be defined by L = T − U , where T is the kinetic energy, U is the potential energy, t is time and ri are axes in a coordinate system. By finite element derivation, the partial differential equation (Eq. 1) is transformed into algebraic equation for 3D free vibration analysis without damping by M δ + Kδ = 0 ,
(2)
roo
f
where M is mass matrix, K is stiffness matrix, and δ is displacement vector. The oscillation is assumed to be simple harmonic motion which can be described by
rep
δ = sin t ,
(3)
where is amplitude of motion and is angular frequency.
lP
Angular frequency is related to frequency f by
=2πf.
(4)
δ = − 2 sin t .
(5)
rna
The second derivative of Eq. (3) is
where
Jo u
By substituting Eqs. (3) and (5) into Eq. (2), Eq. (2) can be written as
(K − M =
2
)
M =0,
N
T
N dV ,
(6)
(7)
V
where is density and
K = B T CB dV = B T CB V , V
where C is stress-strain relation matrix and
(8)
Journal Pre-proof B = LN ,
(9)
where L is strain-displacement relation matrix and N is interpolation function matrix.
2.2 Tetrahedral Element The four-node tetrahedral element type [Daryl 1992] was used in this study. The distribution of variables can be written as = N 11 + N 2 2 + N 3 3 + N 4 4 = NΦ .
(10)
For 3D free vibration analysis, the displacement matrix can be written as δ = Nδ ,
(11)
where
v1
u4
w1
N1
N1
(12) v4
f
w ,
rep
N1 N=
v
roo
δ T = u1
δ T = u
N4
N4
Jo u
rna
lP
where u, v, and w are displacement in x, y, and z directions.
w4 ,
, N 4
(13)
(14)
Journal Pre-proof 3. Numerical Methods Fig. 1 shows the procedure of 3D free vibration analysis in this work. It begins by creating geometry and finite element mesh of the analyzed problem in the commercial finite element software ANSYS®. The information of finite element model is contained in an input file subsequently transferred to the computer programs developed for comparing the classical block Lanczos method, the block Lanczos method with an orthogonality fixing scheme and the Householder method. The mass and stiffness matrices in Eq. (6) ( M and K , respectively) are computed and subsequently transformed into a single standard-form matrix ( A in Eq. (15)). After that, finite element equations are solved by the QL algorithm with implicit shifts to obtain
Jo u
rna
lP
rep
roo
f
natural frequencies and corresponding mode shapes.
Fig.1 Procedure of finite element analysis of free vibration in this study
Journal Pre-proof 3.1 Finite Element Equations in a Standard Form Eq. (6) can be transformed into a standard form of the eigenvalue problem.
(A − I ) Ψ = 0 ,
(15)
where = 2 and I is the identity matrix. The mass matrix in Eq. (6) can be written as M = LL T ,
(16)
where L is a lower triangular matrix which can be obtained by using the Cholesky's symmetric decomposition. Substituting Eq. (16) into Eq. (6) leads to
(K − LL ) Φ = 0 . T
(17)
f
Premultiplying Eq. (17) with L−1 and substituting
in to Eq. (17) leads to
(L
)
KL − T − I Ψ = 0 .
rep
−1
roo
Φ = L− T Ψ
(18)
(19)
Therefore, matrix A of the standard-form equation is (20)
lP
A = L−1KL − T .
rna
3.2 Matrix Tridiagonalization
A tridiagonal matrix is a matrix containing nonzero values only on the main diagonal,
Jo u
the upper subdiagonal and the lower subdiagonal. In finite element analysis of eigenvalue problems, it is common to transform matrices into the tridiagonal form. This leads to the reduction in computational effort required while the eigenvalues of the original matrices are still preserved. In this paper, three tridiagonalization methods, i.e. (1) the classical block Lanczos method, (2) the block Lanczos method with an orthogonality fixing scheme and (3) the Householder method, were compared used for comparison purposes.
Journal Pre-proof 3.2.1 The Block Lanczos Method The Lanczos algorithm was introduced by C. Lanczos [Lanczos 1950]. The algorithm was widely used as a method for reducing a symmetric matrix to the tridiagonal form [Hughes 2000]. The Lanczos algorithm has been further developed to be the so-called block Lanczos algorithm [Montgomery 1995]. The idea of the block Lanczos algorithm is to change matrix A into a block tridiagonal matrix containing small square matrices on the diagonal and on the
upper and lower subdiagonals, as shown in Eq. (21).
0
D2
0
0
D p −1
0
B p −1
0 0 0 . B Tp −1 D p
(21)
f
B 1T
roo
D1 B1 J=0 0 0
From Eq. (21), to compute D i where i = 1, 2,..., p and B j where j = 1, 2,..., p − 1 , the block
rep
Lanczos algorithm can be used as described below.
Set Q 0 = B 0 = 0 , and let Q 1 be an orthonormal column matrix. The iteration block
D i = Q iT AQ j ,
(22)
R i = AQ i − Q i D i − Q i −1B iT−1 .
(23)
rna
lP
Lanczos for i=1,2,...,p-1 is
Jo u
Decompose R i using the QR algorithm to obtain Q i +1 and B i , hence Q i +1B i = R i .
(24)
Finally, we obtain D p = Q Tp AQ p . The next step is to transform the block tridiagonal matrix to tridiagonal matrix T as shown in Eq. (25) by applying the normal Lanczos tridiagonalization.
Journal Pre-proof 1 1 T=0 0 0
1
0
2
0
0
n −1
0
n −1
0 0 0 n −1 n
(25)
To compute i and i in matrix T , set P0 = 0 , 0 = 0 , and let b be an arbitrary starting vector. Then, P1 =
b b
.
(26)
For i = 1, 2,..., n , from orthogonality of matrix P , the steps are as follows (27)
R i = JPi − i Pi − i −1 Pi −1 ,
(28)
roo
f
i = PiT JPi ,
1 = R i .
rep
If i = 0 , quit and end.
Pi +1 =
Ri
i
.
(30)
lP
Else,
(29)
3.2.2 The Block Lanczos Method with an Orthogonality Fixing Scheme
rna
The block Lanczos algorithm may lose orthogonality when computing matrix Q at each iteration. To detect the loss of orthogonality, the componentwise detection algorithm
Jo u
proposed by Qiao et al. [Qiao et al. 2005] was applied in this work to test free vibration real symmetric matrices. At each iteration, if we incorporate rounding error term F j into Eq. (24) where R j is defined in Eq. (23), the equation for j = 1, 2,..., p becomes Q j +1 B j + F j = AQ j − Q j D j − Q j −1 B Tj −1 .
(31)
Journal Pre-proof By moving the term F j to the right hand side, the above equation can be written for steps j and i as Q j +1 B j = AQ j − Q j D j − Q j −1 B Tj −1 − F j ,
(32)
Q i +1 B i = AQ i − Q i D i − Q i −1 B iT−1 − Fi .
(33)
Premultiplying the above two equations with Q iT and Q Tj , respectively, and let W j , i = Q Tj Q i , the equations become W i , j +1 B j = Q iT AQ j − W i , j D j − W i , j −1 B Tj −1 − Q iT F j ,
(34)
W j , i +1 B i = Q Tj AQ i − W j , i D i − W j , i −1 B iT−1 − Q Tj Fi .
(35)
Eq. (35) can be rearranged as
Because the transpose of Q Tj AQ i is Q iT AQ j
j
in Eq. (34), and W iT, j = W j , i , substituting
in Eq. (34) with the transpose of Q Tj AQ i in Eq. (36) to obtain
rep
Q iT AQ
(36)
roo
f
Q Tj AQ i = W j , i +1 B i + W j , i D i + W j , i −1 B iT−1 + Q Tj Fi .
Wi , j +1B j = B iT Wi +1, j + D iT Wi , j + B i −1 Wi −1, j − Wi , j D j − Wi , j −1B Tj −1 + G i , j ,
(37)
lP
where G i , j = FiT Q j − Q iT F j . Also, set W j −1, j −1 = W j , j = I , and W j , j +1B j = (block size ) B 1Φ ,
Φ N (0, 0 .6 ) ,
(38)
rna
where Φ is a normally distributed random vector with zero mean and variance equals to 0.6, and is a small round off unit.
Jo u
In the componentwise detection algorithm [Qiao et al. 2005], for i = 1, 2,..., j − 1 , generate random numbers for G i , j such that G i , j = (B i + B j ) Θ ,
Θ N (0 , 0 .3 ) .
(39)
Journal Pre-proof The above expression implies that each entry in G i , j is the product of B i + B j and a normally distributed random number vector with zero mean and variance equals to 0.3. Hence, by substitute G i , j in Eq. (37), W i , j +1 can be obtained. When the absolute value of a component in W i , j +1 exceeds , find the largest interval of indices li , u i , such that i li , u i , and for all k li , u i , W i , j +1 has at least one component that is larger than a tolerance. A tolerance can be any value between and
, but typically,
7 /8 it is chosen to be .
Experiments show that there is usually one interval during the iteration. The lower
f
bound (li ) is normally 1, so the iteration only needs to search for the upper bound (u i ) . To fix
W k , j +1 = Ω ,
roo
rounding error, after the columns of Q j +1 are orthogonalized against the columns of Q i , set Ω N (0, 1 .5 ) ,
k = 1, 2,..., u i .
(40)
rep
So that the entries of W k , j +1 are adjusted after detecting a loss of orthogonality.
lP
3.2.3 The Householder Method
The Householder method was introduced by A.S. Householder [Householder 1958]. By
rna
applying a sequence of Householder reflection matrices, any symmetric matrix can be converted into the tridiagonal form. Each sequence produces a complete row and column of zeroes apart from the elements within the tridiagonal and subdiagonal. The Householder
Jo u
reflection matrix can be written as H i = I − 2 u i u iT ,
(41)
where u i is a unit vector: ui =
xi − y i xi − y i
.
(42)
Journal Pre-proof Vector x i contains all the off-diagonal entries of each column of matrix. Vector y i contains subdiagonal values of that column. For example, x 1T = 0
a 21
y 1T = 0
a 31
a n1 ,
(43)
0
0 .
(44)
r1
where ri = x 1 = y 1 . The sequence is computed by A 2 = H 1 AH 1
0 a~
23
a~32
a~33
a~
a~
22
n2
0 a~2 n ~ a3 n . a~nn
(46)
f
r1 a~
roo
a11 r 1 =0 0
(45)
n3
Thus, by a single Householder transformation, matrix A is converted into a similar matrix A 2
rep
whose the first row and column are in tridiagonal form. By repeating the process on the lower right (n − 1) (n − 1) submatrix of A 2 , the next matrix A 3 whose the first two rows and
rna
matrix is achieved.
lP
columns are in tridiagonal form are constructed. The process is repeated until the tridiagonal
3.3 Solver for Finite Element Equations
Jo u
After tridiagonalization, the QL algorithm with implicit shifts [Press et al. 2002] is used to solve the finite element equations. The original QL algorithm decomposes matrix A into a product of an orthogonal matrix Q and a lower triangular matrix L . Therefore, matrix A can be expressed as
A n = Q nL n .
Then calculate the new matrix A by switching the order of Q n and L n in Eq. (47)
(47)
Journal Pre-proof A n +1 = L n Q n .
(48)
The sequence when n → , A n +1 converges to a diagonal matrix which the eigenvalues appear on the diagonal. Convergence can be accelerated by using the shifting technique. The acceleration begins by letting k n be any value, then A n − knI = Q nL n
(49)
A n +1 = L n Q n + k n I .
(50)
Similarly,
So, at each iteration, the eigenvalues are shifted based on the chosen value of k n . Again, when n → , A n +1 converges to a diagonal matrix. Hence, the eigenvalues of the original matrix A
roo
f
are entries on the diagonal of matrix A n +1 subtracted by the shifted values k n at all iterations.
rep
4. Results
Three test cases were used to evaluate the classical block Lanczos (BL) method, the block Lanczos method with an orthogonality fixing scheme (BLO) and the Householder
lP
method. Natural frequencies (eigenvalues of free vibration problems) computed by the finite
rna
element codes developed were compared with exact or experimental solutions.
4.1 Solution Accuracy
Jo u
4.1.1 Solid Cube
The first test case is the free vibration analysis of a solid cube fixed at one side (Fig.2) [Petyt 2010]. Its dimension is 0.254 0.254 0.254 m3. The material properties of the solid cube are Young's Modulus of 68 .95 10 9 N/m2, Poisson's ratio of 0.3 and density of 2,560
Journal Pre-proof kg/m3. Fig. 2 shows a tetrahedral finite element mesh of the cube consisting of 651 nodes and 3,105 elements.
0.254 m
f
Fixed bottomsurface (u, v, w = 0)
roo
Fig. 2 Problem statement and finite element model of the solid cube
rep
Table 1 compares the numerical results obtained from the three finite element codes developed with exact solutions. The exact solutions of the cube's natural frequencies for modes 1, 3, 4 and 5 are 2.212, 3.020, 5.239 and 5.915, respectively [Petyt (2010)]. For mode 1, the
lP
natural frequency values computed by the three finite element codes developed are 2 .2231
rna
corresponding with the exact solution of 2.212. For mode 2, the natural frequency values computed by the BLO and Householder methods are equally 2 .2280 . However, the BL method provides the value of 2 .2231 for
Jo u
modes 2 and 3, which is similar to the value for mode 1. The BL method gives the value of 2 .2244 for mode 4 which still does not match the value of 2 .2280 computed by the BLO
and Householder methods. It is not until modes 5-7 where the BL method gives the solution of 2 .2280 . From these computed values from the BL method, it can be concluded that the BL
method gives repeated solutions for mode 1 three times and a pseudo solution for mode 2.
Journal Pre-proof Therefore, the first three repeated solutions must be considered as the natural frequency for mode 1. The forth solution should be neglected and hence the repeated fifth to seventh solutions can be regarded as the natural frequency for mode 2. For mode 3 and higher modes, we can consider the solutions from the BL method as the way we do for modes 1 and 2. For other mesh sizes used, the BL method gives natural frequency values as shown in Table 2.
Table 1 Comparison of natural frequency values obtained from different methods. The highlighted areas indicate pseudo solutions. Natural frequency results (kHz) BL (651 nodes)
BLO Householder (651 nodes)
Exact Solution
1
2.223122766273
2.223122766357
2.212
2
2.223122766318
2.228006674919
-
3
2.223122766344
3.108044422828
4
2.224442868885
5.235033399656
5
2.228006674826
5.910830950668
6
2.228006674880
5.922559328955
7
2.228006674902
7.451773096101
8
2.237011501732
8.817662344362
9
3.108044422770
9.204074398796
10
3.108044422819
9.349554489457
11
3.108044422837
9.378075276070
12
3.108871789966
13
3.108871790000
14
3.108871790023
15
5.235033399620
16
5.235033399637
10.411883473055
17
5.235033399647
11.565741202671
18
5.235053022983
11.766066238042
19
5.910830950636
11.883370422864
21
roo 3.020
5.239 5.915
rep
lP
9.452224042220 9.829216402140
10.333548638576 10.345002378012
rna
Jo u
20
f
Mode no. given by computer program
5.910830950661
12.040328323601
5.910830968241
13.259677585933
Journal Pre-proof Table 2 Natural frequency values for modes 1-5 computed by the FEM program using the basic block Lanczos method. Pseudo solutions are highlighted. Mode no. given by computer program 1 2
3 4 5 6 7
Number of nodes 67 1 2 3 4 5
97 1
236
423 1
1
651 1
1000 1
2.437832655815 2.335489618980 2.266372843190 2.241599214872 2.223122766273 2.211758240963 2
2.476365069336 2.348254906872 2.266372843193 2.241599214897 2.223122766318 2.211758240968 2
3.829660914548 2.363350310018 2.276089857858 2.244733116264 2.223122766344 2.211758240982 2
3
5.394682797078 3.546944445138 2.276089857861 2.249288061804 2.224442868885 2.211758240999 2
3
6.564503205184 3.547257660608 3.274194421664 2.249288061841 2.228006674826 2.211758241008 4
5
2
3
5.333476178205 3.274194422769 3.174730055251 2.228006674880 2.215305954915 6.315265080811 3.276901551950 3.174730635354 2.228006674902 2.215305954922
8
3.276901636771 3.174730635377 2.237011501732 2.215305954931 4
9 10
3
5.274798763407 3.448616564223 3.108044422770 2.215305954935 5.274798763413 3.448628932741 3.108044422819 2.215305954967
5
11
3
6.084888196557 3.448628932752 3.108044422837 3.068655115535 4
6.084888203727 5.249929937648 3.108871789966 3.068655115548
13
5.250131976911 3.108871790000 3.068655115573
roo
f
12
14
5.250131976972 3.108871790023 3.068655115602
5
15 16
4
5.969031855977 5.235033399620 3.068655115608 5.969032182487 5.235033399637 3.068655115612
17
rep
18 19 20 21
22
lP
23
5.235053022983 5.225901970504 5
5.910830950636 5.225901970512 5.910830950661 5.225901970512 5.910830968241 5.225901970552 5
5.871401652430 5.871401652436 5.871401652439
rna
24
4
5.235033399647 5.225901970503
Fig. 3 graphically shows how the BL method gives repeated natural frequency values.
Jo u
For example, the frequencies for modes 1 to 10 (the first five BL's circles) are relatively closed to each other. These values can be combined to obtain a single frequency for mode 1. Similarly, the next values (the sixth to the eighth BL's circles) can also be combined to obtain a single frequency for the next mode. For this reason, the plotted solutions of the BL method in Fig. 3 are very different to the solutions from other methods.
Journal Pre-proof 35
BL BLO Householder
Natural frequency (kHz)
30 25 20 15 10 5 0 0
20
40
60
80
100
Mode
Fig. 3 Result comparison for the cubic box problem using the basic block Lanczos method,
f
the orthogonalized block Lanczos method, and the Householder method
roo
Figs. 4-7 show the natural frequency values obtained from the three finite element codes developed when different mesh sizes are used. Since the original study (Petyt 2010) did not
rep
provide an exact solution for mode 2, the comparison for this mode is omitted. The solutions from the BL method shown in Figs. 4-7 have been manipulated to get rid of repeated and pseudo solutions. Figs. 4-7 show that the solutions obtained from the developed FEM programs
Jo u
rna
lP
converge to the exact solutions with increasing node numbers.
Fig. 4 Comparison of mode frequency values for mode 1 (swaying mode)
Journal Pre-proof
rep
roo
f
Fig. 5 Comparison of mode frequency values for mode 3 (torsion mode)
Jo u
rna
lP
Fig. 6 Comparison of mode frequency values for mode 4 (longitudinal mode)
Fig. 7 Comparison of mode frequency values for mode 5 (swaying mode)
Journal Pre-proof 4.1.2 Cantilever Beam The second case is the free vibration analysis of a cantilever beam in six different boundary conditions: Clamped-Free (C-F), Clamped-Roller (C-R), Clamped-Clamped (C-C), Clamped-Pined (C-P), Pined-Pined (P-P), and Pined-Roller (P-R) as demonstrated in Figure 8. The beam’s dimension is 0.30 0.61 3.66 m3. Its material properties are Young's Modulus of 2.068 1011 N/m2, Poisson's ratio of 0.3 and density of 8,058 kg/m3. The figure also indicates a tetrahedral finite element mesh of the beam at 2,624 nodes and 11,906 elements. Fixed surface
Boundary conditions
rep
C-F
roo
3.66 m
C-R
0.61 m 0.30 m
C-P
P-P
P-R
lP
C-C
f
(u, v, w = 0 )
rna
Fig. 8 Problem statement and finite element model of the cantilever beam
To ensure the versatility and accuracy of the model, Table 3 compares the numerical
Jo u
results obtained in six boundary conditions: Clamped-Free (C-F), Clamped-Roller (C-R), Clamped-Clamped (C-C), Clamped-Pined (C-P), Pined-Pined (P-P), and Pined-Roller (P-R). For BL method, similar natural frequency values can be observed in every condition. For example, under C-F condition, BL method provides the value of 19 .4367 for modes 2 and 3, which is similar to the value for mode 1. Four out of six boundary conditions (C-F, C-R, C-P,
Journal Pre-proof and P-R) provide pseudo solutions. In contrast, BLO and Householder methods indicate neither the repeated natural frequency values nor the pseudo solutions, as shown in Table 4.
Table 3 Comparison of natural frequency values obtained in six boundary conditions under BL method. The highlighted areas indicate pseudo solutions. M ode no. given by computer program 1
Natural frequency results (Hz) C-F
C-R
C-C
C-P
P-P
1 19.43671718 1 80.72083427 1 119.66682261 1 93.74972459 1 72.78994754 1 61.50725303
2
19.43672065
80.72083543
119.66682312
93.74972546
72.78994851
3
19.43672129
80.72083578
119.66682314
93.74972584
72.78994913
4
2 37.33992738
80.72083600 37.33992927 2 148.95783167
119.66682602
5
119.66682632
93.74972596 163.44125741 2 142.65347230
6
37.35098933
7
38.07562885
148.95819055
207.97214783
164.65117372
8
38.07562989 3 118.12942941
148.96217077
207.97214791
164.65117637
148.96217123
207.97215775
142.66887198 164.65117725 3 181.84562758
10
118.12943007 3 208.31595516
207.97215800
183.71879152
11
118.12943028
208.31595552 3 312.84645353 3 214.06310960
12
148.07517590
208.31595565
312.84645366
181.84562807 214.06311004 4 205.30234303
13
148.07517632
214.30044020
312.84645392
214.06311736
14
148.07517645 4 174.53874000
214.30044049
312.84709457
16
214.30044060
96.11568211
142.65347293
100.26035988 100.26036026
roo
f
142.65347274
291.13249465 291.13249495
257.76373524
355.70183195
291.13258219
257.76373544
19
209.45628438
335.30230870
355.74634222
20
209.45628465
335.30230888
21
209.45628477
26
335.30230895 355.74634235 5 347.54621618 5 493.74155632
rna
25
lP
355.70183178
283.77160374
291.13258239 355.74634231 5 396.69288473 396.83551738 396.83551760
347.54621640
493.74155641
396.83554337
347.54621649
493.74155648
396.83554337
493.74238200 493.74238207
181.55950412
193.28062634 211.72282440
257.76373501
283.77160366
24
145.06003282
3
181.55950424 4 205.30234341 193.28062621
174.53874052 5 209.40336726
23
145.06003263
181.84562789
18
22
92.52623236
214.06311754 205.30234357 312.84709468 4 291.13241149 5 257.74138838
174.53874044 4 283.77160337 4 355.70183176
Jo u
17
rep
15
61.50725367
61.50726535 2 72.79002429 92.52623212
148.95783216 2 207.97214775 2 164.65117324
9
P-R
211.72282448
5
215.66229796
215.66229816
Journal Pre-proof Table 4 Comparison of natural frequency values obtained in six boundary conditions under BLO and Householder methods. M ode no. given by computer program
Natural frequency results (Hz) C-F
C-R
C-C
C-P
P-P
P-R
19.437
80.721
119.667
93.750
72.790
61.507
2
37.340
148.958
207.972
164.651
142.653
92.526
3
118.129
208.316
312.846
214.063
181.846
181.560
4
174.539
283.772
355.702
291.132
205.302
193.281
5
209.456
347.546
493.742
396.836
257.764
215.662
1
Under C-F boundary condition, the BL solutions appear totally different from the BL and BLO solutions (Fig. 9). This is because the BL method gives repeated solutions similar to
roo
f
Fig. 3. Figs. 10-12 show that the solutions obtained from the developed FEM programs converge to the exact solutions (Petyt 2010) with increasing node numbers.
rep
BL BLO Householder
6
lP
4
2
rna
Natural frequency (kHz)
8
0
Jo u
0
20
40
60
80
100
Mode
Fig. 9 Result comparison for the cantilever beam problem using the classical block Lanczos method, the orthogonalized block Lanczos method, and the Householder method
Journal Pre-proof
rep
roo
f
Fig. 10 Comparison of mode frequency values for mode 1 (first bending in x direction)
Jo u
rna
lP
Fig. 11 Comparison of mode frequency values for mode 2 (first bending in y direction)
Fig. 12 Comparison of mode frequency values for mode 3 (second bending in x direction)
Journal Pre-proof 4.1.3 Anvil The third test case is the free vibration analysis of an unconstrained anvil (Fig.13) (To 1982). Its dimension is 0.461 0.461 0.152 m3. The material properties of the anvil are Young's Modulus of 2.07 1011 N/m2, Poisson's ratio of 0.3 and density of 7,860 kg/m3. Fig. 13 shows a tetrahedral finite element mesh of the anvil consisting of 2,461 nodes and 11,857 elements.
m
roo
0.461
f
0.
46
1
m
0.152 m
rep
Fig. 13 Problem statement and finite element model of the anvil
In Fig. 14, the BL solutions are different from others because the BL method gives
lP
repeated solutions like in Figs. 3 and 9. Figs. 15-18 show that the solutions obtained from the
Jo u
numbers.
rna
developed FEM programs converge to the experimental solutions with increasing node
Journal Pre-proof
Natural frequency (kHz)
20
BL BLO Householder
15
10
5
0 0
20
40
60
80
100
Mode
Fig. 14 Result comparison for the anvil problem using the basic block Lanczos method, the
rna
lP
rep
roo
f
orthogonalized block Lanczos method, and the Householder method
Jo u
Fig. 15 Comparison of mode frequency values for mode 1 (twist mode)
Journal Pre-proof
rep
roo
f
Fig. 16 Comparison of mode frequency values for mode 2 (saddle mode)
Jo u
rna
lP
Fig. 17 Comparison of mode frequency values for mode 3 (umbrella mode)
Fig. 18 Comparison of mode frequency values for mode 4 (in-plane shear mode)
Journal Pre-proof 4.1.4 Turbine blade The final test case is the free vibration analysis of a turbine blade fixed at one end (Fig. 19) [Salama 1976]. Its dimension is shown in Fig.19. The material properties of the beam are Young's Modulus of 2.05 1011 N/m2, Poisson's ratio of 0.3 and density of 7,800 kg/m3. Fig.
o
3.17 cm 5
30
1.65 1 cm
19 shows a tetrahedral finite element mesh consisting of 3,000 nodes and 12,323 elements.
rep
7.62 cm
roo
f
5.08 cm
Fixed surface
(u, v, w = 0)
lP
Fig. 19 Problem statement and finite element model of the turbine blade
rna
In Fig. 20, the BL solutions appear totally different from the BL and BLO solutions because the BL method gives repeated solutions like in Figs. 3, 9 and 14. Figs. 21-23 show that
Jo u
the solutions obtained from the developed FEM programs converge to the experimental solutions with increasing node numbers.
Journal Pre-proof 180
BL BLO Householder
Natural frequency (kHz)
150 120 90 60 30 0 0
20
40
60
80
100
Mode
Fig. 20 Result comparison for the turbine blade problem using the classical block Lanczos
rna
lP
rep
roo
f
method, the orthogonalized block Lanczos method, and the Householder method
Jo u
Fig. 21 Comparison of mode frequency values for mode 1 (twist mode)
Journal Pre-proof
lP
rep
roo
f
Fig. 22 Comparison of mode frequency values for mode 2 (saddle mode)
rna
Fig. 23 Comparison of mode frequency values for mode 3 (umbrella mode)
4.2 Computational Time
Jo u
Concerning the four test problems, there were many numbers of finite element nodes used for each problem. Fig. 24 shows the time used in the matrix tridiagonalization step for each number of nodes regardless the test cases. For relatively-small problems (number of nodes < about 250), the BL and BLO methods used similar time but significantly more than the Householder method did. For larger problems (number of nodes > about 250), the Householder method used the longest time followed by the BLO and BL methods, respectively.
Journal Pre-proof 1
Time ratio (t/t max )
0.1 0.01
0.001 BL BLO
0.0001
Householder 0.00001 0.000001 0
500
1000
1500
2000
2500
3000
Number of nodes
roo
f
Fig. 24 Comparison of computation time used for different problem sizes
5. Discussion
The results described in the previous section show that there are two problems with the
rep
numerical solutions from the BL method. The first problem is the repeated natural frequency values in Table 1 and Table 3. The second problem is the pseudo modes, which does not match
lP
any exact nor experimental solutions, as performed under a single condition or different boundary conditions (Table 1 and Table 3). These problems attribute to the loss of matrix
rna
orthogonality during tridiagonalization subsequently causing rounding errors in later steps. On the other hand, the BLO and Householder methods do not suffer these problems (Table 1, Table 4). This is because the orthogonality fixing scheme in the BLO method is effective in detecting
Jo u
and fixing the loss of matrix orthogonality. For the Householder method, it can be deduced that there are no matrices losing orthogonality. In terms of computational time used, the results show that the BL method used less computational time than the Householder method in all cases except in case of relatively-small problems. This trend is in agreement with the literature [de Silva 2007] mentioning that the
Journal Pre-proof Householder method suits small problems better than the BL method. For the comparison of the BL method and the BLO method, the BLO method always used significantly more time than the BL method because some computer resource was used to perform the orthogonality fixing scheme. Nevertheless, the BLO method still used less computational time than the Householder method in all cases except in case of relatively-small problems. Considering the solution accuracy and computational time aspects, the BLO method is a good choice for the matrix tridiagonalization operation in 3D free vibration analysis. This is because the BLO method provides solution accuracy as good as the Householder method but uses significantly less computational time in most cases.
roo
f
6. Conclusions
In this work, we developed three computer programs for 3D finite element analysis of
rep
free vibration by using three matrix tridiagonalization methods: (1) the BL method, (2) the BLO method and (3) the Householder method. Numerical solutions obtained from the three developed programs were compared with exact or experimental solutions. Considering the
lP
analysis on several examples (e.g., solid cube, cantilever beam) and under different boundary conditions, the following conclusions can be drawn:
rna
1) The BL method experiences two problems, which are repeated natural frequency values and pseudo modes.
Jo u
2) The Householder method does not suffer the loss of matrix orthogonality and its accuracy is similar to that of the BLO method, which is satisfactory. 3) In terms of computational time, the BLO method always take significantly more hours than the BL method because some computer resource was occupied in performing the orthogonality fixing scheme.
Journal Pre-proof 4) The BLO method provides solution accuracy as good as the Householder method but uses significantly less computational time in most cases. Therefore, the BLO method is a good choice for the matrix tridiagonalization operation in 3D free vibration analysis.
Acknowledgements The authors would like to acknowledge the support of National Metal and Materials Technology Center (MTEC), the National Science and Technology Development Agency (NSTDA), Thailand for granting the Project No. P1100755.
f
References W.S. Bottega, Engineering Vibrations, Taylor and Francis Group, Florida, 2006.
[2]
A. Bunse-Gerstner, W. Gragg, Singular value decompositions of complex symmetric
roo
[1]
[3]
rep
matrices, J. Comput. Appl. Math. 21 (1988) 41–54. J.K. Cullum, R.A. Willoughby, (1985). Lanczos Algorithms for Large Symmetric Eigenvalue Computations. Vol. 1: Theory. Birkhauser, Boston, MA. Reissued by
lP
SIAM, Philadelphia, 2002.
L.L. Daryl, A First Course in Finite Element Method, PWS, USA, 1992.
[5]
J.H. Ginberg, Mechanical and Structural Vibrations: Theory and Applications, John
rna
[4]
Wiley and Sons, Inc., New York, 2001. G.H. Golub, R. Underwood, The Block Lanczos Method for Computing Eigenvalues, In
Jo u
[6]
Mathematical Software III, J.R. Rice (Ed.), Academic Press, New York, 1977. [7]
C. Guo and S. Qiao. A Stable Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS 03-08-SQ, Department of Computing and Software, McMaster University,Ontario, Canada, 2003.
Journal Pre-proof [8]
A.S. Householder, Unitary triangularization of a nonsymmetric matrix, JACM 5(4) (1958) 339-342.
[9]
T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey 07632, 2000.
[10] H. Iguchi, M. Natori, and H. Imai, Reorthogonalization in the Block Lanczos Algorithm, ISE-TR-92-100, 1992. [11] H. Imai, M. Natori, E. Kawamura, A new reorthogonalization in the Lanczos algorithm, J. lnfo. Proc. 14 (1991) 56-59. [12] C. Lanczos, An iteration method for solution of the eigenvalue problem of linear differential and integral operation, Journal of research of the national bureau of standards
roo
f
45(4) (1950) 255-282.
[13] F.T. Luk, S. Qiao, A Fast Singular Value Algorithm for Hankel Matrices, in: V. Olshevsky (Ed.), Fast Algorithms for Structured Matrices: Theory and Applications,
rep
Contemporary Mathematics, vol. 323, American Mathematical Society, 2003. [14] P.L. Montgomery, A block Lanczos algorithm for finding dependencies over GF(2),
lP
EUROCRYPT'95 921, Springer-Verlag. (1995) 106-120. [15] C. Paige, Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix,
rna
J. Inst. Math. Appl. 18 (1976) 341–349. [16] B.N. Parlett, J.K. Reid, Tracking the progress of the Lanczos algorithm for large symmetric eigenproblem, IMA J. Num. Anal. 1 (1981) 135-155.
Jo u
[17] B.N. Parlett, D. Scott, The Lanczos algorithm with selective orthogonalization, Math. Comp. 33 (1979) 217-238. [18] M. Petyt, Introduction to Finite Element Vibration Analysis, Cambridge University Press, New York, 2010.
Journal Pre-proof [19] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 90: the Art of Parallel Scientic Computing Volume 2 of Fortran Numerical Recipes, Cambridge University Press, New York, 2002. [20] S. Qiao, G. Liu, W. Xu, Block Lanczos tridiagonalization of complex symmetric matrices, Advanced signal processing algorithms, architectures, and implementations XV, Proc. of the SPIE, 5910 (2005) 285-295. [21] C.W.S. To, Application of the finite element method for the evaluation of velocity response of anvil, J. Sound Vib. 84(4) (1982) 529-548. [22] A.M. Salama, Finite Element Dynamic Analysis of Blade Packets and Bladed Disc Assemblies, PhD thesis, University of Southampton, UK., 1976.
roo
f
[23] C.W. de Silva, Computer Techniques in Vibration, CRC Press, Boca Raton, FL, 4-34, 2007.
[24] H.D. Simon, The Lanczos algorithm with partial reorthogonalization, Math. Comp. 42
rep
(1984) 115-142.
[25] S. Yu. Fialko, E.Z. Kriksunov, V.S. Karpilovskyy, A block Lanczos method with spectral
lP
transformations for natural vibrations and seismic analysis of large structures in SCAD software, in Proc. CMM-2003 - Computer Methods in Mechanics, Gliwice, Poland,
Jo u
rna
(2003) 129-130.
PII: DOI: Reference:
S0378-4754(19)30391-X https://doi.org/10.1016/j.matcom.2019.12.017 MATCOM 4921
To appear in:
Mathematics and Computers in Simulation
Received date : 16 February 2018 Revised date : 29 September 2019 Accepted date : 27 December 2019 Please cite this article as: A. Malatip, N. Prasomsuk, C. Siriparu et al., An efficient matrix tridiagonalization method for 3D finite element analysis of free vibration, Mathematics and Computers in Simulation (2020), doi: https://doi.org/10.1016/j.matcom.2019.12.017. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). ⃝ Published by Elsevier B.V. All rights reserved.
Journal Pre-proof
An Efficient Matrix Tridiagonalization Method for 3D Finite Element Analysis of Free Vibration
A. Malatip1,*, N. Prasomsuk2, C. Siriparu2, N. Paoprasert1 and S. Otarawanna1
1
Computer-Aided Engineering Research Team (CAET), Engineering Design and Advanced Manufacturing
Research Group, National Metal and Materials Technology Center (MTEC), National Science and Technology Development Agency (NSTDA), 114 Thailand Science Park, Thanon Pahonyothin, Tambon Khlong Nueng, Amphoe Khlong Luang, Pathum Thani 12120, Thailand
Department of Mechanical Engineering, Faculty of Engineering, Chulalongkorn University,
f
2
roo
Bangkok 10330, Thailand
rep
* Corresponding author, Email: [email protected]
Abstract
lP
In order to avoid mechanical resonance, a vibrating structure needs to be designed such that its working frequency interval is sufficiently far from its natural frequencies. Natural
rna
frequencies of a mechanical system are obtained from free vibration analysis commonly done by the finite element method. One important step in the analysis is matrix tridiagonalization
Jo u
commonly performed by the block Lanczos method. However, the classical block Lanczos method suffers numerical instability due to the loss of matrix orthogonality. In this paper, we demonstrate the implementation of the block Lanczos method with an orthogonality fixing scheme for 3D free vibration problems. The solution accuracy and computational time of this method were compared with those of the classical block Lanczos method and the Householder method. The results show that the block Lanczos method with the orthogonality fixing scheme
Journal Pre-proof employed in this work can effectively avoid the numerical instability due to the loss of matrix orthogonality. Furthermore, the block Lanczos method with the orthogonality fixing scheme provides solution accuracy as good as the Householder method while using significantly less computational time.
Keywords: matrix tridiagonalization; eigenvalue problem; free vibration; block Lanczos; orthogonality
1. Introduction One of the most common problems for a vibrating structure is mechanical resonance.
f
When the resonance occurs, the structure oscillates with higher amplitude than it does at other
roo
frequencies. A vibrating system is prone to mechanical resonance when the frequency of its oscillation matches one of the system's natural frequencies. Mechanical resonance can cause
rep
unwanted violent motions and even structural failure in structures including buildings, bridges, machines and vehicles. Therefore, design engineers must ensure that the working frequency interval of the mechanical system is sufficiently far from its natural frequencies.
lP
Natural frequencies of a mechanical system can be determined from free vibration analysis [Ginberg 2001; Bottega 2006]. For 3D complex geometries, this is commonly
rna
performed by the finite element method. In the finite element derivation of 3D free vibration problems, the mass and stiffness matrices can be transformed into a single standard-form
Jo u
matrix. If the size of the standard-form matrix is relatively small, eigenvalues and eigenvectors of the matrix can be directly solved within a relatively short period of time. For 3D free vibration problems, the size of the standard-form matrix is relatively large. Therefore, the matrix needs to be transformed into the tridiagonalized form (containing non-zero elements only on the main diagonal, and the first diagonal below and above the main diagonal) by a tridiagonalization algorithm.
Journal Pre-proof For the tridiagonalization operation, there are three renowned methods which are the Householder method [Bunse-Gerstner and Gragg 1988], the Lanczos method [Lanczos 1950, Parlett and Reid 1981, Cullum and Willoughby 1985] and the block Lanczos method [Golub and Underwood 1977]. The Householder method has a disadvantage when dealing with a sparse or structured matrix because the method destroys the matrix’s sparsity or structure [Guo and Qiao 2003]. On the contrary, the Lanczos method and the block Lanczos method can retain the matrix’s sparsity or structure when used as they involves only matrix-vector multiplication [Luk and Qiao 2003]. However, cancellation and roundoff errors make the Lanczos method and the block Lanczos method unstable due to the loss of orthogonality [Paige 1976, Iguchi et al. 1992].
roo
f
The loss of orthogonality in the Lanczos method can be fixed by different schemes, e.g. selective orthogonalization (SO) [Parlett and Scott 1979], partial orthogonalization (PRO) [Simon 1984] and reorthogonalization with improved convergence-check (RIC) [Imai et al.
rep
1991]. However, for relatively-large matrices, the block Lanczos method [Montgomery 1995] is preferred due to its better performance compared to the Lanczos method. This is because the
lP
block Lanczos method exploits computer memory hierarchies [Yu et al. 2003, Qiao et al. 2005]. To stabilize the block Lanczos algorithm, a componentwise technique for detecting and fixing
rna
the loss of orthogonality has been proposed and demonstrated on general eigenvalue problems [Qiao et al. 2005].
In this work, it aims to investigate the effectiveness of the block Lanczos method with
Jo u
this orthogonality fixing scheme for 3D free vibration problems. In the following sections, the implementation of this method is demonstrated and compared with the classical block Lanczos method and the Householder method [Householder 1958] in terms of solution accuracy and computational time.
Journal Pre-proof 2. Finite Element Derivation 2.1 Finite Element Equations The differential equation describing the motion in the Euler-Lagrange system is d L L − = 0. dt ri ri
(1)
The Lagrangian for a system of particles can be defined by L = T − U , where T is the kinetic energy, U is the potential energy, t is time and ri are axes in a coordinate system. By finite element derivation, the partial differential equation (Eq. 1) is transformed into algebraic equation for 3D free vibration analysis without damping by M δ + Kδ = 0 ,
(2)
roo
f
where M is mass matrix, K is stiffness matrix, and δ is displacement vector. The oscillation is assumed to be simple harmonic motion which can be described by
rep
δ = sin t ,
(3)
where is amplitude of motion and is angular frequency.
lP
Angular frequency is related to frequency f by
=2πf.
(4)
δ = − 2 sin t .
(5)
rna
The second derivative of Eq. (3) is
where
Jo u
By substituting Eqs. (3) and (5) into Eq. (2), Eq. (2) can be written as
(K − M =
2
)
M =0,
N
T
N dV ,
(6)
(7)
V
where is density and
K = B T CB dV = B T CB V , V
where C is stress-strain relation matrix and
(8)
Journal Pre-proof B = LN ,
(9)
where L is strain-displacement relation matrix and N is interpolation function matrix.
2.2 Tetrahedral Element The four-node tetrahedral element type [Daryl 1992] was used in this study. The distribution of variables can be written as = N 11 + N 2 2 + N 3 3 + N 4 4 = NΦ .
(10)
For 3D free vibration analysis, the displacement matrix can be written as δ = Nδ ,
(11)
where
v1
u4
w1
N1
N1
(12) v4
f
w ,
rep
N1 N=
v
roo
δ T = u1
δ T = u
N4
N4
Jo u
rna
lP
where u, v, and w are displacement in x, y, and z directions.
w4 ,
, N 4
(13)
(14)
Journal Pre-proof 3. Numerical Methods Fig. 1 shows the procedure of 3D free vibration analysis in this work. It begins by creating geometry and finite element mesh of the analyzed problem in the commercial finite element software ANSYS®. The information of finite element model is contained in an input file subsequently transferred to the computer programs developed for comparing the classical block Lanczos method, the block Lanczos method with an orthogonality fixing scheme and the Householder method. The mass and stiffness matrices in Eq. (6) ( M and K , respectively) are computed and subsequently transformed into a single standard-form matrix ( A in Eq. (15)). After that, finite element equations are solved by the QL algorithm with implicit shifts to obtain
Jo u
rna
lP
rep
roo
f
natural frequencies and corresponding mode shapes.
Fig.1 Procedure of finite element analysis of free vibration in this study
Journal Pre-proof 3.1 Finite Element Equations in a Standard Form Eq. (6) can be transformed into a standard form of the eigenvalue problem.
(A − I ) Ψ = 0 ,
(15)
where = 2 and I is the identity matrix. The mass matrix in Eq. (6) can be written as M = LL T ,
(16)
where L is a lower triangular matrix which can be obtained by using the Cholesky's symmetric decomposition. Substituting Eq. (16) into Eq. (6) leads to
(K − LL ) Φ = 0 . T
(17)
f
Premultiplying Eq. (17) with L−1 and substituting
in to Eq. (17) leads to
(L
)
KL − T − I Ψ = 0 .
rep
−1
roo
Φ = L− T Ψ
(18)
(19)
Therefore, matrix A of the standard-form equation is (20)
lP
A = L−1KL − T .
rna
3.2 Matrix Tridiagonalization
A tridiagonal matrix is a matrix containing nonzero values only on the main diagonal,
Jo u
the upper subdiagonal and the lower subdiagonal. In finite element analysis of eigenvalue problems, it is common to transform matrices into the tridiagonal form. This leads to the reduction in computational effort required while the eigenvalues of the original matrices are still preserved. In this paper, three tridiagonalization methods, i.e. (1) the classical block Lanczos method, (2) the block Lanczos method with an orthogonality fixing scheme and (3) the Householder method, were compared used for comparison purposes.
Journal Pre-proof 3.2.1 The Block Lanczos Method The Lanczos algorithm was introduced by C. Lanczos [Lanczos 1950]. The algorithm was widely used as a method for reducing a symmetric matrix to the tridiagonal form [Hughes 2000]. The Lanczos algorithm has been further developed to be the so-called block Lanczos algorithm [Montgomery 1995]. The idea of the block Lanczos algorithm is to change matrix A into a block tridiagonal matrix containing small square matrices on the diagonal and on the
upper and lower subdiagonals, as shown in Eq. (21).
0
D2
0
0
D p −1
0
B p −1
0 0 0 . B Tp −1 D p
(21)
f
B 1T
roo
D1 B1 J=0 0 0
From Eq. (21), to compute D i where i = 1, 2,..., p and B j where j = 1, 2,..., p − 1 , the block
rep
Lanczos algorithm can be used as described below.
Set Q 0 = B 0 = 0 , and let Q 1 be an orthonormal column matrix. The iteration block
D i = Q iT AQ j ,
(22)
R i = AQ i − Q i D i − Q i −1B iT−1 .
(23)
rna
lP
Lanczos for i=1,2,...,p-1 is
Jo u
Decompose R i using the QR algorithm to obtain Q i +1 and B i , hence Q i +1B i = R i .
(24)
Finally, we obtain D p = Q Tp AQ p . The next step is to transform the block tridiagonal matrix to tridiagonal matrix T as shown in Eq. (25) by applying the normal Lanczos tridiagonalization.
Journal Pre-proof 1 1 T=0 0 0
1
0
2
0
0
n −1
0
n −1
0 0 0 n −1 n
(25)
To compute i and i in matrix T , set P0 = 0 , 0 = 0 , and let b be an arbitrary starting vector. Then, P1 =
b b
.
(26)
For i = 1, 2,..., n , from orthogonality of matrix P , the steps are as follows (27)
R i = JPi − i Pi − i −1 Pi −1 ,
(28)
roo
f
i = PiT JPi ,
1 = R i .
rep
If i = 0 , quit and end.
Pi +1 =
Ri
i
.
(30)
lP
Else,
(29)
3.2.2 The Block Lanczos Method with an Orthogonality Fixing Scheme
rna
The block Lanczos algorithm may lose orthogonality when computing matrix Q at each iteration. To detect the loss of orthogonality, the componentwise detection algorithm
Jo u
proposed by Qiao et al. [Qiao et al. 2005] was applied in this work to test free vibration real symmetric matrices. At each iteration, if we incorporate rounding error term F j into Eq. (24) where R j is defined in Eq. (23), the equation for j = 1, 2,..., p becomes Q j +1 B j + F j = AQ j − Q j D j − Q j −1 B Tj −1 .
(31)
Journal Pre-proof By moving the term F j to the right hand side, the above equation can be written for steps j and i as Q j +1 B j = AQ j − Q j D j − Q j −1 B Tj −1 − F j ,
(32)
Q i +1 B i = AQ i − Q i D i − Q i −1 B iT−1 − Fi .
(33)
Premultiplying the above two equations with Q iT and Q Tj , respectively, and let W j , i = Q Tj Q i , the equations become W i , j +1 B j = Q iT AQ j − W i , j D j − W i , j −1 B Tj −1 − Q iT F j ,
(34)
W j , i +1 B i = Q Tj AQ i − W j , i D i − W j , i −1 B iT−1 − Q Tj Fi .
(35)
Eq. (35) can be rearranged as
Because the transpose of Q Tj AQ i is Q iT AQ j
j
in Eq. (34), and W iT, j = W j , i , substituting
in Eq. (34) with the transpose of Q Tj AQ i in Eq. (36) to obtain
rep
Q iT AQ
(36)
roo
f
Q Tj AQ i = W j , i +1 B i + W j , i D i + W j , i −1 B iT−1 + Q Tj Fi .
Wi , j +1B j = B iT Wi +1, j + D iT Wi , j + B i −1 Wi −1, j − Wi , j D j − Wi , j −1B Tj −1 + G i , j ,
(37)
lP
where G i , j = FiT Q j − Q iT F j . Also, set W j −1, j −1 = W j , j = I , and W j , j +1B j = (block size ) B 1Φ ,
Φ N (0, 0 .6 ) ,
(38)
rna
where Φ is a normally distributed random vector with zero mean and variance equals to 0.6, and is a small round off unit.
Jo u
In the componentwise detection algorithm [Qiao et al. 2005], for i = 1, 2,..., j − 1 , generate random numbers for G i , j such that G i , j = (B i + B j ) Θ ,
Θ N (0 , 0 .3 ) .
(39)
Journal Pre-proof The above expression implies that each entry in G i , j is the product of B i + B j and a normally distributed random number vector with zero mean and variance equals to 0.3. Hence, by substitute G i , j in Eq. (37), W i , j +1 can be obtained. When the absolute value of a component in W i , j +1 exceeds , find the largest interval of indices li , u i , such that i li , u i , and for all k li , u i , W i , j +1 has at least one component that is larger than a tolerance. A tolerance can be any value between and
, but typically,
7 /8 it is chosen to be .
Experiments show that there is usually one interval during the iteration. The lower
f
bound (li ) is normally 1, so the iteration only needs to search for the upper bound (u i ) . To fix
W k , j +1 = Ω ,
roo
rounding error, after the columns of Q j +1 are orthogonalized against the columns of Q i , set Ω N (0, 1 .5 ) ,
k = 1, 2,..., u i .
(40)
rep
So that the entries of W k , j +1 are adjusted after detecting a loss of orthogonality.
lP
3.2.3 The Householder Method
The Householder method was introduced by A.S. Householder [Householder 1958]. By
rna
applying a sequence of Householder reflection matrices, any symmetric matrix can be converted into the tridiagonal form. Each sequence produces a complete row and column of zeroes apart from the elements within the tridiagonal and subdiagonal. The Householder
Jo u
reflection matrix can be written as H i = I − 2 u i u iT ,
(41)
where u i is a unit vector: ui =
xi − y i xi − y i
.
(42)
Journal Pre-proof Vector x i contains all the off-diagonal entries of each column of matrix. Vector y i contains subdiagonal values of that column. For example, x 1T = 0
a 21
y 1T = 0
a 31
a n1 ,
(43)
0
0 .
(44)
r1
where ri = x 1 = y 1 . The sequence is computed by A 2 = H 1 AH 1
0 a~
23
a~32
a~33
a~
a~
22
n2
0 a~2 n ~ a3 n . a~nn
(46)
f
r1 a~
roo
a11 r 1 =0 0
(45)
n3
Thus, by a single Householder transformation, matrix A is converted into a similar matrix A 2
rep
whose the first row and column are in tridiagonal form. By repeating the process on the lower right (n − 1) (n − 1) submatrix of A 2 , the next matrix A 3 whose the first two rows and
rna
matrix is achieved.
lP
columns are in tridiagonal form are constructed. The process is repeated until the tridiagonal
3.3 Solver for Finite Element Equations
Jo u
After tridiagonalization, the QL algorithm with implicit shifts [Press et al. 2002] is used to solve the finite element equations. The original QL algorithm decomposes matrix A into a product of an orthogonal matrix Q and a lower triangular matrix L . Therefore, matrix A can be expressed as
A n = Q nL n .
Then calculate the new matrix A by switching the order of Q n and L n in Eq. (47)
(47)
Journal Pre-proof A n +1 = L n Q n .
(48)
The sequence when n → , A n +1 converges to a diagonal matrix which the eigenvalues appear on the diagonal. Convergence can be accelerated by using the shifting technique. The acceleration begins by letting k n be any value, then A n − knI = Q nL n
(49)
A n +1 = L n Q n + k n I .
(50)
Similarly,
So, at each iteration, the eigenvalues are shifted based on the chosen value of k n . Again, when n → , A n +1 converges to a diagonal matrix. Hence, the eigenvalues of the original matrix A
roo
f
are entries on the diagonal of matrix A n +1 subtracted by the shifted values k n at all iterations.
rep
4. Results
Three test cases were used to evaluate the classical block Lanczos (BL) method, the block Lanczos method with an orthogonality fixing scheme (BLO) and the Householder
lP
method. Natural frequencies (eigenvalues of free vibration problems) computed by the finite
rna
element codes developed were compared with exact or experimental solutions.
4.1 Solution Accuracy
Jo u
4.1.1 Solid Cube
The first test case is the free vibration analysis of a solid cube fixed at one side (Fig.2) [Petyt 2010]. Its dimension is 0.254 0.254 0.254 m3. The material properties of the solid cube are Young's Modulus of 68 .95 10 9 N/m2, Poisson's ratio of 0.3 and density of 2,560
Journal Pre-proof kg/m3. Fig. 2 shows a tetrahedral finite element mesh of the cube consisting of 651 nodes and 3,105 elements.
0.254 m
f
Fixed bottomsurface (u, v, w = 0)
roo
Fig. 2 Problem statement and finite element model of the solid cube
rep
Table 1 compares the numerical results obtained from the three finite element codes developed with exact solutions. The exact solutions of the cube's natural frequencies for modes 1, 3, 4 and 5 are 2.212, 3.020, 5.239 and 5.915, respectively [Petyt (2010)]. For mode 1, the
lP
natural frequency values computed by the three finite element codes developed are 2 .2231
rna
corresponding with the exact solution of 2.212. For mode 2, the natural frequency values computed by the BLO and Householder methods are equally 2 .2280 . However, the BL method provides the value of 2 .2231 for
Jo u
modes 2 and 3, which is similar to the value for mode 1. The BL method gives the value of 2 .2244 for mode 4 which still does not match the value of 2 .2280 computed by the BLO
and Householder methods. It is not until modes 5-7 where the BL method gives the solution of 2 .2280 . From these computed values from the BL method, it can be concluded that the BL
method gives repeated solutions for mode 1 three times and a pseudo solution for mode 2.
Journal Pre-proof Therefore, the first three repeated solutions must be considered as the natural frequency for mode 1. The forth solution should be neglected and hence the repeated fifth to seventh solutions can be regarded as the natural frequency for mode 2. For mode 3 and higher modes, we can consider the solutions from the BL method as the way we do for modes 1 and 2. For other mesh sizes used, the BL method gives natural frequency values as shown in Table 2.
Table 1 Comparison of natural frequency values obtained from different methods. The highlighted areas indicate pseudo solutions. Natural frequency results (kHz) BL (651 nodes)
BLO Householder (651 nodes)
Exact Solution
1
2.223122766273
2.223122766357
2.212
2
2.223122766318
2.228006674919
-
3
2.223122766344
3.108044422828
4
2.224442868885
5.235033399656
5
2.228006674826
5.910830950668
6
2.228006674880
5.922559328955
7
2.228006674902
7.451773096101
8
2.237011501732
8.817662344362
9
3.108044422770
9.204074398796
10
3.108044422819
9.349554489457
11
3.108044422837
9.378075276070
12
3.108871789966
13
3.108871790000
14
3.108871790023
15
5.235033399620
16
5.235033399637
10.411883473055
17
5.235033399647
11.565741202671
18
5.235053022983
11.766066238042
19
5.910830950636
11.883370422864
21
roo 3.020
5.239 5.915
rep
lP
9.452224042220 9.829216402140
10.333548638576 10.345002378012
rna
Jo u
20
f
Mode no. given by computer program
5.910830950661
12.040328323601
5.910830968241
13.259677585933
Journal Pre-proof Table 2 Natural frequency values for modes 1-5 computed by the FEM program using the basic block Lanczos method. Pseudo solutions are highlighted. Mode no. given by computer program 1 2
3 4 5 6 7
Number of nodes 67 1 2 3 4 5
97 1
236
423 1
1
651 1
1000 1
2.437832655815 2.335489618980 2.266372843190 2.241599214872 2.223122766273 2.211758240963 2
2.476365069336 2.348254906872 2.266372843193 2.241599214897 2.223122766318 2.211758240968 2
3.829660914548 2.363350310018 2.276089857858 2.244733116264 2.223122766344 2.211758240982 2
3
5.394682797078 3.546944445138 2.276089857861 2.249288061804 2.224442868885 2.211758240999 2
3
6.564503205184 3.547257660608 3.274194421664 2.249288061841 2.228006674826 2.211758241008 4
5
2
3
5.333476178205 3.274194422769 3.174730055251 2.228006674880 2.215305954915 6.315265080811 3.276901551950 3.174730635354 2.228006674902 2.215305954922
8
3.276901636771 3.174730635377 2.237011501732 2.215305954931 4
9 10
3
5.274798763407 3.448616564223 3.108044422770 2.215305954935 5.274798763413 3.448628932741 3.108044422819 2.215305954967
5
11
3
6.084888196557 3.448628932752 3.108044422837 3.068655115535 4
6.084888203727 5.249929937648 3.108871789966 3.068655115548
13
5.250131976911 3.108871790000 3.068655115573
roo
f
12
14
5.250131976972 3.108871790023 3.068655115602
5
15 16
4
5.969031855977 5.235033399620 3.068655115608 5.969032182487 5.235033399637 3.068655115612
17
rep
18 19 20 21
22
lP
23
5.235053022983 5.225901970504 5
5.910830950636 5.225901970512 5.910830950661 5.225901970512 5.910830968241 5.225901970552 5
5.871401652430 5.871401652436 5.871401652439
rna
24
4
5.235033399647 5.225901970503
Fig. 3 graphically shows how the BL method gives repeated natural frequency values.
Jo u
For example, the frequencies for modes 1 to 10 (the first five BL's circles) are relatively closed to each other. These values can be combined to obtain a single frequency for mode 1. Similarly, the next values (the sixth to the eighth BL's circles) can also be combined to obtain a single frequency for the next mode. For this reason, the plotted solutions of the BL method in Fig. 3 are very different to the solutions from other methods.
Journal Pre-proof 35
BL BLO Householder
Natural frequency (kHz)
30 25 20 15 10 5 0 0
20
40
60
80
100
Mode
Fig. 3 Result comparison for the cubic box problem using the basic block Lanczos method,
f
the orthogonalized block Lanczos method, and the Householder method
roo
Figs. 4-7 show the natural frequency values obtained from the three finite element codes developed when different mesh sizes are used. Since the original study (Petyt 2010) did not
rep
provide an exact solution for mode 2, the comparison for this mode is omitted. The solutions from the BL method shown in Figs. 4-7 have been manipulated to get rid of repeated and pseudo solutions. Figs. 4-7 show that the solutions obtained from the developed FEM programs
Jo u
rna
lP
converge to the exact solutions with increasing node numbers.
Fig. 4 Comparison of mode frequency values for mode 1 (swaying mode)
Journal Pre-proof
rep
roo
f
Fig. 5 Comparison of mode frequency values for mode 3 (torsion mode)
Jo u
rna
lP
Fig. 6 Comparison of mode frequency values for mode 4 (longitudinal mode)
Fig. 7 Comparison of mode frequency values for mode 5 (swaying mode)
Journal Pre-proof 4.1.2 Cantilever Beam The second case is the free vibration analysis of a cantilever beam in six different boundary conditions: Clamped-Free (C-F), Clamped-Roller (C-R), Clamped-Clamped (C-C), Clamped-Pined (C-P), Pined-Pined (P-P), and Pined-Roller (P-R) as demonstrated in Figure 8. The beam’s dimension is 0.30 0.61 3.66 m3. Its material properties are Young's Modulus of 2.068 1011 N/m2, Poisson's ratio of 0.3 and density of 8,058 kg/m3. The figure also indicates a tetrahedral finite element mesh of the beam at 2,624 nodes and 11,906 elements. Fixed surface
Boundary conditions
rep
C-F
roo
3.66 m
C-R
0.61 m 0.30 m
C-P
P-P
P-R
lP
C-C
f
(u, v, w = 0 )
rna
Fig. 8 Problem statement and finite element model of the cantilever beam
To ensure the versatility and accuracy of the model, Table 3 compares the numerical
Jo u
results obtained in six boundary conditions: Clamped-Free (C-F), Clamped-Roller (C-R), Clamped-Clamped (C-C), Clamped-Pined (C-P), Pined-Pined (P-P), and Pined-Roller (P-R). For BL method, similar natural frequency values can be observed in every condition. For example, under C-F condition, BL method provides the value of 19 .4367 for modes 2 and 3, which is similar to the value for mode 1. Four out of six boundary conditions (C-F, C-R, C-P,
Journal Pre-proof and P-R) provide pseudo solutions. In contrast, BLO and Householder methods indicate neither the repeated natural frequency values nor the pseudo solutions, as shown in Table 4.
Table 3 Comparison of natural frequency values obtained in six boundary conditions under BL method. The highlighted areas indicate pseudo solutions. M ode no. given by computer program 1
Natural frequency results (Hz) C-F
C-R
C-C
C-P
P-P
1 19.43671718 1 80.72083427 1 119.66682261 1 93.74972459 1 72.78994754 1 61.50725303
2
19.43672065
80.72083543
119.66682312
93.74972546
72.78994851
3
19.43672129
80.72083578
119.66682314
93.74972584
72.78994913
4
2 37.33992738
80.72083600 37.33992927 2 148.95783167
119.66682602
5
119.66682632
93.74972596 163.44125741 2 142.65347230
6
37.35098933
7
38.07562885
148.95819055
207.97214783
164.65117372
8
38.07562989 3 118.12942941
148.96217077
207.97214791
164.65117637
148.96217123
207.97215775
142.66887198 164.65117725 3 181.84562758
10
118.12943007 3 208.31595516
207.97215800
183.71879152
11
118.12943028
208.31595552 3 312.84645353 3 214.06310960
12
148.07517590
208.31595565
312.84645366
181.84562807 214.06311004 4 205.30234303
13
148.07517632
214.30044020
312.84645392
214.06311736
14
148.07517645 4 174.53874000
214.30044049
312.84709457
16
214.30044060
96.11568211
142.65347293
100.26035988 100.26036026
roo
f
142.65347274
291.13249465 291.13249495
257.76373524
355.70183195
291.13258219
257.76373544
19
209.45628438
335.30230870
355.74634222
20
209.45628465
335.30230888
21
209.45628477
26
335.30230895 355.74634235 5 347.54621618 5 493.74155632
rna
25
lP
355.70183178
283.77160374
291.13258239 355.74634231 5 396.69288473 396.83551738 396.83551760
347.54621640
493.74155641
396.83554337
347.54621649
493.74155648
396.83554337
493.74238200 493.74238207
181.55950412
193.28062634 211.72282440
257.76373501
283.77160366
24
145.06003282
3
181.55950424 4 205.30234341 193.28062621
174.53874052 5 209.40336726
23
145.06003263
181.84562789
18
22
92.52623236
214.06311754 205.30234357 312.84709468 4 291.13241149 5 257.74138838
174.53874044 4 283.77160337 4 355.70183176
Jo u
17
rep
15
61.50725367
61.50726535 2 72.79002429 92.52623212
148.95783216 2 207.97214775 2 164.65117324
9
P-R
211.72282448
5
215.66229796
215.66229816
Journal Pre-proof Table 4 Comparison of natural frequency values obtained in six boundary conditions under BLO and Householder methods. M ode no. given by computer program
Natural frequency results (Hz) C-F
C-R
C-C
C-P
P-P
P-R
19.437
80.721
119.667
93.750
72.790
61.507
2
37.340
148.958
207.972
164.651
142.653
92.526
3
118.129
208.316
312.846
214.063
181.846
181.560
4
174.539
283.772
355.702
291.132
205.302
193.281
5
209.456
347.546
493.742
396.836
257.764
215.662
1
Under C-F boundary condition, the BL solutions appear totally different from the BL and BLO solutions (Fig. 9). This is because the BL method gives repeated solutions similar to
roo
f
Fig. 3. Figs. 10-12 show that the solutions obtained from the developed FEM programs converge to the exact solutions (Petyt 2010) with increasing node numbers.
rep
BL BLO Householder
6
lP
4
2
rna
Natural frequency (kHz)
8
0
Jo u
0
20
40
60
80
100
Mode
Fig. 9 Result comparison for the cantilever beam problem using the classical block Lanczos method, the orthogonalized block Lanczos method, and the Householder method
Journal Pre-proof
rep
roo
f
Fig. 10 Comparison of mode frequency values for mode 1 (first bending in x direction)
Jo u
rna
lP
Fig. 11 Comparison of mode frequency values for mode 2 (first bending in y direction)
Fig. 12 Comparison of mode frequency values for mode 3 (second bending in x direction)
Journal Pre-proof 4.1.3 Anvil The third test case is the free vibration analysis of an unconstrained anvil (Fig.13) (To 1982). Its dimension is 0.461 0.461 0.152 m3. The material properties of the anvil are Young's Modulus of 2.07 1011 N/m2, Poisson's ratio of 0.3 and density of 7,860 kg/m3. Fig. 13 shows a tetrahedral finite element mesh of the anvil consisting of 2,461 nodes and 11,857 elements.
m
roo
0.461
f
0.
46
1
m
0.152 m
rep
Fig. 13 Problem statement and finite element model of the anvil
In Fig. 14, the BL solutions are different from others because the BL method gives
lP
repeated solutions like in Figs. 3 and 9. Figs. 15-18 show that the solutions obtained from the
Jo u
numbers.
rna
developed FEM programs converge to the experimental solutions with increasing node
Journal Pre-proof
Natural frequency (kHz)
20
BL BLO Householder
15
10
5
0 0
20
40
60
80
100
Mode
Fig. 14 Result comparison for the anvil problem using the basic block Lanczos method, the
rna
lP
rep
roo
f
orthogonalized block Lanczos method, and the Householder method
Jo u
Fig. 15 Comparison of mode frequency values for mode 1 (twist mode)
Journal Pre-proof
rep
roo
f
Fig. 16 Comparison of mode frequency values for mode 2 (saddle mode)
Jo u
rna
lP
Fig. 17 Comparison of mode frequency values for mode 3 (umbrella mode)
Fig. 18 Comparison of mode frequency values for mode 4 (in-plane shear mode)
Journal Pre-proof 4.1.4 Turbine blade The final test case is the free vibration analysis of a turbine blade fixed at one end (Fig. 19) [Salama 1976]. Its dimension is shown in Fig.19. The material properties of the beam are Young's Modulus of 2.05 1011 N/m2, Poisson's ratio of 0.3 and density of 7,800 kg/m3. Fig.
o
3.17 cm 5
30
1.65 1 cm
19 shows a tetrahedral finite element mesh consisting of 3,000 nodes and 12,323 elements.
rep
7.62 cm
roo
f
5.08 cm
Fixed surface
(u, v, w = 0)
lP
Fig. 19 Problem statement and finite element model of the turbine blade
rna
In Fig. 20, the BL solutions appear totally different from the BL and BLO solutions because the BL method gives repeated solutions like in Figs. 3, 9 and 14. Figs. 21-23 show that
Jo u
the solutions obtained from the developed FEM programs converge to the experimental solutions with increasing node numbers.
Journal Pre-proof 180
BL BLO Householder
Natural frequency (kHz)
150 120 90 60 30 0 0
20
40
60
80
100
Mode
Fig. 20 Result comparison for the turbine blade problem using the classical block Lanczos
rna
lP
rep
roo
f
method, the orthogonalized block Lanczos method, and the Householder method
Jo u
Fig. 21 Comparison of mode frequency values for mode 1 (twist mode)
Journal Pre-proof
lP
rep
roo
f
Fig. 22 Comparison of mode frequency values for mode 2 (saddle mode)
rna
Fig. 23 Comparison of mode frequency values for mode 3 (umbrella mode)
4.2 Computational Time
Jo u
Concerning the four test problems, there were many numbers of finite element nodes used for each problem. Fig. 24 shows the time used in the matrix tridiagonalization step for each number of nodes regardless the test cases. For relatively-small problems (number of nodes < about 250), the BL and BLO methods used similar time but significantly more than the Householder method did. For larger problems (number of nodes > about 250), the Householder method used the longest time followed by the BLO and BL methods, respectively.
Journal Pre-proof 1
Time ratio (t/t max )
0.1 0.01
0.001 BL BLO
0.0001
Householder 0.00001 0.000001 0
500
1000
1500
2000
2500
3000
Number of nodes
roo
f
Fig. 24 Comparison of computation time used for different problem sizes
5. Discussion
The results described in the previous section show that there are two problems with the
rep
numerical solutions from the BL method. The first problem is the repeated natural frequency values in Table 1 and Table 3. The second problem is the pseudo modes, which does not match
lP
any exact nor experimental solutions, as performed under a single condition or different boundary conditions (Table 1 and Table 3). These problems attribute to the loss of matrix
rna
orthogonality during tridiagonalization subsequently causing rounding errors in later steps. On the other hand, the BLO and Householder methods do not suffer these problems (Table 1, Table 4). This is because the orthogonality fixing scheme in the BLO method is effective in detecting
Jo u
and fixing the loss of matrix orthogonality. For the Householder method, it can be deduced that there are no matrices losing orthogonality. In terms of computational time used, the results show that the BL method used less computational time than the Householder method in all cases except in case of relatively-small problems. This trend is in agreement with the literature [de Silva 2007] mentioning that the
Journal Pre-proof Householder method suits small problems better than the BL method. For the comparison of the BL method and the BLO method, the BLO method always used significantly more time than the BL method because some computer resource was used to perform the orthogonality fixing scheme. Nevertheless, the BLO method still used less computational time than the Householder method in all cases except in case of relatively-small problems. Considering the solution accuracy and computational time aspects, the BLO method is a good choice for the matrix tridiagonalization operation in 3D free vibration analysis. This is because the BLO method provides solution accuracy as good as the Householder method but uses significantly less computational time in most cases.
roo
f
6. Conclusions
In this work, we developed three computer programs for 3D finite element analysis of
rep
free vibration by using three matrix tridiagonalization methods: (1) the BL method, (2) the BLO method and (3) the Householder method. Numerical solutions obtained from the three developed programs were compared with exact or experimental solutions. Considering the
lP
analysis on several examples (e.g., solid cube, cantilever beam) and under different boundary conditions, the following conclusions can be drawn:
rna
1) The BL method experiences two problems, which are repeated natural frequency values and pseudo modes.
Jo u
2) The Householder method does not suffer the loss of matrix orthogonality and its accuracy is similar to that of the BLO method, which is satisfactory. 3) In terms of computational time, the BLO method always take significantly more hours than the BL method because some computer resource was occupied in performing the orthogonality fixing scheme.
Journal Pre-proof 4) The BLO method provides solution accuracy as good as the Householder method but uses significantly less computational time in most cases. Therefore, the BLO method is a good choice for the matrix tridiagonalization operation in 3D free vibration analysis.
Acknowledgements The authors would like to acknowledge the support of National Metal and Materials Technology Center (MTEC), the National Science and Technology Development Agency (NSTDA), Thailand for granting the Project No. P1100755.
f
References W.S. Bottega, Engineering Vibrations, Taylor and Francis Group, Florida, 2006.
[2]
A. Bunse-Gerstner, W. Gragg, Singular value decompositions of complex symmetric
roo
[1]
[3]
rep
matrices, J. Comput. Appl. Math. 21 (1988) 41–54. J.K. Cullum, R.A. Willoughby, (1985). Lanczos Algorithms for Large Symmetric Eigenvalue Computations. Vol. 1: Theory. Birkhauser, Boston, MA. Reissued by
lP
SIAM, Philadelphia, 2002.
L.L. Daryl, A First Course in Finite Element Method, PWS, USA, 1992.
[5]
J.H. Ginberg, Mechanical and Structural Vibrations: Theory and Applications, John
rna
[4]
Wiley and Sons, Inc., New York, 2001. G.H. Golub, R. Underwood, The Block Lanczos Method for Computing Eigenvalues, In
Jo u
[6]
Mathematical Software III, J.R. Rice (Ed.), Academic Press, New York, 1977. [7]
C. Guo and S. Qiao. A Stable Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS 03-08-SQ, Department of Computing and Software, McMaster University,Ontario, Canada, 2003.
Journal Pre-proof [8]
A.S. Householder, Unitary triangularization of a nonsymmetric matrix, JACM 5(4) (1958) 339-342.
[9]
T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey 07632, 2000.
[10] H. Iguchi, M. Natori, and H. Imai, Reorthogonalization in the Block Lanczos Algorithm, ISE-TR-92-100, 1992. [11] H. Imai, M. Natori, E. Kawamura, A new reorthogonalization in the Lanczos algorithm, J. lnfo. Proc. 14 (1991) 56-59. [12] C. Lanczos, An iteration method for solution of the eigenvalue problem of linear differential and integral operation, Journal of research of the national bureau of standards
roo
f
45(4) (1950) 255-282.
[13] F.T. Luk, S. Qiao, A Fast Singular Value Algorithm for Hankel Matrices, in: V. Olshevsky (Ed.), Fast Algorithms for Structured Matrices: Theory and Applications,
rep
Contemporary Mathematics, vol. 323, American Mathematical Society, 2003. [14] P.L. Montgomery, A block Lanczos algorithm for finding dependencies over GF(2),
lP
EUROCRYPT'95 921, Springer-Verlag. (1995) 106-120. [15] C. Paige, Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix,
rna
J. Inst. Math. Appl. 18 (1976) 341–349. [16] B.N. Parlett, J.K. Reid, Tracking the progress of the Lanczos algorithm for large symmetric eigenproblem, IMA J. Num. Anal. 1 (1981) 135-155.
Jo u
[17] B.N. Parlett, D. Scott, The Lanczos algorithm with selective orthogonalization, Math. Comp. 33 (1979) 217-238. [18] M. Petyt, Introduction to Finite Element Vibration Analysis, Cambridge University Press, New York, 2010.
Journal Pre-proof [19] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 90: the Art of Parallel Scientic Computing Volume 2 of Fortran Numerical Recipes, Cambridge University Press, New York, 2002. [20] S. Qiao, G. Liu, W. Xu, Block Lanczos tridiagonalization of complex symmetric matrices, Advanced signal processing algorithms, architectures, and implementations XV, Proc. of the SPIE, 5910 (2005) 285-295. [21] C.W.S. To, Application of the finite element method for the evaluation of velocity response of anvil, J. Sound Vib. 84(4) (1982) 529-548. [22] A.M. Salama, Finite Element Dynamic Analysis of Blade Packets and Bladed Disc Assemblies, PhD thesis, University of Southampton, UK., 1976.
roo
f
[23] C.W. de Silva, Computer Techniques in Vibration, CRC Press, Boca Raton, FL, 4-34, 2007.
[24] H.D. Simon, The Lanczos algorithm with partial reorthogonalization, Math. Comp. 42
rep
(1984) 115-142.
[25] S. Yu. Fialko, E.Z. Kriksunov, V.S. Karpilovskyy, A block Lanczos method with spectral
lP
transformations for natural vibrations and seismic analysis of large structures in SCAD software, in Proc. CMM-2003 - Computer Methods in Mechanics, Gliwice, Poland,
Jo u
rna
(2003) 129-130.