Comparision and study of numerical methods for dynamic response evaluation of sdof

Nội dung text: Comparision and study of numerical methods for dynamic response evaluation of sdof

1. Tạp chí Khoa học và Công nghệ, Số 43A, 2020 COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC RESPONSE EVALUATION OF SDOF THAI PHUONG TRUC Department of Civil Engineering, Industrial University of Ho Chi Minh City, Vietnam; thaiphuongtruc@iuh.edu.vn Abstract. Written for senior-year undergraduates and first-year graduate students with solid backgrounds in differential and integral calculus, this paper is oriented toward engineers and applied mathematicians. Consequently, this paper should be useful to senior-year undergraduates the finite element method [1]. The scaled direct approach is adopted for this purpose and each step in the finite element solution process is given in full detail. For this reason, all students must be exposed to (and indeed should master). This paper provides the general framework for the development of nearly all (nonstructural) finite element models. The finite element method of analysis is a very powerful, modern computational tool. Applications range from deformation and stress analysis of automotive, aircraft, building, and bridge structures to field analysis of beat flux, fluid flow, magnetic flux, seepage, and other flow problems. This paper presents study and comparison of numerical methods which are used for evaluation of dynamic response. A Single Degree of Freedom (SDF)-linear problem is solved by means of Newmark’s Average acceleration method [2], Linear acceleration method [2], Central Difference method [6,7] with the help of MATLAB. The advantages, disadvantages, relative precision and applicability of these numerical methods are discussed throughout the analysis. Keywords. Finite element method, central difference method, Newmark’s constant average acceleration method, Newmark’s linear acceleration method. 1 INTRODUCTION The basic idea behind the finite element method is to divide the structure, body, or region being analyzed into a large number of finite elements, or simply elements. These elements may be one, two, or three dimensional. In 1941, Alexander Hrennikoff (was born in Russia, graduated from the Institute of Railway Engineering in Moscow) he developed the lattice analogy which models membrane and plate bending of structures as a lattice framework [1]. In the early 1960s, engineers used the method for approximate solution of problems in stress analysis, fluid flow, heat transfer. A book by Argyris in 1955 on energy theorems and matrix methods laid a foundation for further developments in finite element studies. In 1956, Turner et al derived stiffness matrices for truss beam and other elements. Today, the developments in mainframe computers and availability of powerful microcomputers has brought this method within reach of students and engineers working in industries. © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
2. COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC 89 RESPONSE EVALUATION OF SDOF Figure 1: FEM mesh created by an analyst prior to finding a solution to a magnetic problem using FEM software. 2 PROBLEM SOLVING METHODOLOGY Basic steps of finite element analysis [9] Step 1 - Discretization of real continuum or structure Step 2 - Identify primary unknown quantity Step 3 - Interpolation functions and the derivation of Interpolation functions Step 4 - Derivation of Element equation Step 5 - Derive overall Stiffness Equation The static equilibrium equation: Kq  R (2.1) Dynamic equilibrium equations Mqt   Cqt   Kqt  Rt  (2.2) where: K , M  , C are overall stiffness matrix, mass matrix, damping matrix; q t , nodal force vector q t , q t  - vector displacement node, velocity vector, acceleration vector (dynamics problems) R - nodal force vector (static problem) R t  - Nodal applied force vector (dynamic problem) K' TKTT (2.3) e  e  e  e M' TMTT (2.4) e  e  e  e T which  e is transfer matrix: L . 0 L  . T  e 0 . L © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
3. 90 COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC RESPONSE EVALUATION OF SDOF and L is the cosine of a matrix: ' ' ' cos x,,, x cos x y cos x z lxx' m xy ' n xz ' ' ' ' L  lyx' m yy ' n yz ' cosyx ,,, cosyy cosyz lzx' m zy ' n zz ' ' ' ' cos z,,, x cos z y cos z z where l ' , m , n ' - is the direction cosine of x axis; xx xy' xz l , m , n - is the direction cosine of y axis yx' yy' yz' l ' , m , n ' - is the direction cosine of z axis zx zy' zz Global stiffness matrix can be expressed as K DT E D dV (2.5)  e   eee 0    V where: D - matrix deformation - displacement is derived from the relationship  e deformation-displacement  D q  B q e e  e   e  e (2.6) D  B  e   e (2.7)  - The differential operator matrix is determined from elastic theory B U B q  e - The function matrix in the displacement function expression e e  e , This function represents the variation of the displacement within the element. E  0 e - is the elastic compliance matrix, and is symmetric. The general form of the  E  strain–stress relation e 0 e  e Mass matrix of elements in the local coordinate system M BT B dV (2.8)  e   e  e V - The density of a material is defined as its mass per unit volume C - is damping matrix of the system is combined linearly from the stiffness matrix K and mass matrix M  can be determined by the formula: C  M  K with is the mass proportional Rayleigh damping coefficient  is the stiffness proportional Rayleigh damping coefficient Classical Rayleigh damping results in different damping ratios for different response frequencies, according to the equation: © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
4. COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC 91 RESPONSE EVALUATION OF SDOF 2  (2.9) 2  is the damping ratio (a value of 1 corresponds to critical damping)  is the response frequency in rad/s. It can be seen from this that the mass proportional term gives a damping ratio inversely proportional to response frequency and the stiffness proportional term gives a damping ratio linearly proportional to response frequency. Nodal force vector: R' TRT (2.10) e  e e Dynamic load changes with time f t : Rt' TT Rft (2.11) e  e e Nodal force vector of elements in the local coordinate system RR VS R R o  e e e  e (2.12) where: T R B pdV V e  e V  V T R B p dS SS e  e  S T o Ro D E  dV  e  e o  e e V V , S - volume and area of the element; T pV  p Vx p Vy p Vz  - force distributed by volume T pS  p Sx p Sy  - force distributed by area 0 0 0 0 0 0 0 T     - Initial forced distortion in the element e x y z xy zy xz  Step 6 – Solve for the unknowns Step 7 – Interpretation of results A differential equation is used for calculating structure by Finite Element Method. It is tough to solve this equation by analytical method 2.1. The central difference method [6] Step 1 - Initial calculations R  C q  K q q  0 0 0 (2.13) 0 M  t 2 q  q t  q q  (2.14) 1 0 02 0 M  C K (2.15) 2 t 2 t M  C a  (2.16) t 2 2 t © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
5. 92 COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC RESPONSE EVALUATION OF SDOF 2M  b  K  (2.17) t 2 Step 2 - Calculate the values with each time step Rˆ R   a q   b q  (2.18) i i i 1 i R q i (2.19) i 1 K q q q  i 1 i 1 (2.20) i 2 t q 2 q q q  i 1 i i 1 (2.21) i t 2 Step 3: Calculate for next time steps Example 1: Determine the displacement of the single degree of freedom system subjected to dynamic loads for the figure 2 Figure 2: A load changes with time.  Solution: The mass of the object: m 0,2533 , stiffness k 10 , Damping Coefficient c 0,1592 . the initial parameters: At the beginning of the survey t 0 ; initial displacement and velocity q0 0 ; q0 0 ; Step 1: R c q k q q 0 0 0 0 0 m t 2 q q t q q 0 1 0 02 0 m c k 26,13 t 2 2 t m c a 24,53 t 2 2 t © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
6. COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC 93 RESPONSE EVALUATION OF SDOF 2m b k 40,66 t 2 Step 2: Calculated for each time step: Riiiii R a. q 1 b . q R 24,53. q i 1 40,66. q i R R q i i i 1 k 26,13 Step 3: Table 1: Calculate for next time steps (numerical results show in table results). ti Ri qi 1 qi Ri qi 1 0,0 0,0000 0,0000 0,0000 0,0000 0,0000 0,1 5,0000 0,0000 0,0000 5,0000 0,1914 0,2 8,6602 0,0000 0,1914 16,4419 0,6293 0,3 10,0000 0,1914 0,6293 30,8934 1,1825 0,4 8,6602 0,6293 1,1825 41,3001 1,5808 0,5 5,0000 1,1825 1,5808 40,2649 1,5412 0,6 0,0000 1,5808 1,5412 23,8809 0,9141 0,7 0,0000 1,5412 0,9141 -0,6456 -0,0247 0,8 0,0000 0,9141 -0,0247 -23,4309 -0,8968 0,9 0,0000 -0,0247 -0,8968 -35,8598 -1,3726 1,0 0,0000 -0,8968 -1,3726 -33,8058 -1,2940 2.2. Newmark’s constant average acceleration method [9] Step 1: R c q k q q 0 0 0 0 0 0 m The time step is the incremental change in time t 0,1sec 2 4 kk c m 114,5 t t 2 4 a m 2. c 10,45; b 2 m 0,5066 t Step 2: Calculated for each time step: Riiiii R a. q b . q R 10.45 q i 0.5066 q i Ri R i qi ki 114,5 2 q q2 q 20 q 2 q i t ii ii 4 qi q iii t. q 2. q 400 q iii 0,1 q 2. q t 2 qi 1 q i q i ; qi 1 q i q i ; qi 1 q i q i © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
7. 94 COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC RESPONSE EVALUATION OF SDOF Step 3: Table 2: Calculate for next time steps (numerical results show in table results). ti Ri qi Ri Ri qi qi qi qi qi 0,0 0,0000 0,0000 5,0000 5,0000 0,0437 0,8733 17,4666 0,0000 0,0000 0,1 5,0000 17,4666 3,6603 21,6356 0,1890 2,0323 5,7137 0,8733 0,0437 0,2 8,6602 23,1803 1,3398 43,4485 0,3794 1,7776 -10,8087 2,9057 0,2326 0,3 10,0000 12,3724 -1,3397 53,8708 0,4705 0,0428 -23,8893 4,6833 0,6121 0,4 8,6602 -11,5169 -3,6602 39,8948 0,3484 -2,483 -26,6442 4,7261 1,0825 0,5 5,0000 -38,1611 -5,0000 -0,9009 -0,0079 -4,641 -16,5122 2,2422 1,4309 0,6 0,0000 -54,6733 0,0000 -52,7740 -0,4609 -4,418 20,9716 -2,3995 1,4231 0,7 0,0000 -33,7017 0,0000 -88,3275 -0,7714 -1,791 31,5787 -6,8133 0,9622 0,8 0,0000 -2,1229 0,0000 -91,0486 -0,7952 1,3159 30,5646 -8,6095 0,1908 0,9 0,0000 28,4417 0,0000 -61,8123 -0,5398 3,7907 18,9297 -7,2936 -0,6044 1,0 0,0000 47,3714 0,0000 -3,5029 -1,1442 2.3. Newmark’s linear acceleration method [9] Step 1: R c q k q q 0 0 0 0 0 0 m t 0,1sec 3 6 kk c m 166,8 t t 2 6 t a mc 3. 15,68 ; bm 3 c 0,7679 t 2 Step 2: Riiiii R a. q b . q R 15,68 q i 0,7679 q i Ri R i qi ki 166,8 3 t q q3 q q 30 q 3 q 0,05 q i t iii2 ii i 6 qi q iii t. q 3. q 600 q iii 0,1 q 3. q t 2 qi 1 q i q i ; qi 1 q i q i ; qi 1 q i q i © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
8. COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC 95 RESPONSE EVALUATION OF SDOF Step 3: Table 3: Calculate for next time steps (numerical results show in table results). ti Ri qi Ri Ri qi qi qi qi qi 0,0 0,0000 0,0000 5,0000 5,0000 0,0300 0,8995 17,9903 0,0000 0,0000 0,1 5,0000 17,9903 3,6603 31,5749 0,1893 2,0824 5,6666 0,8995 0,0300 0,2 8,6602 23,6569 1,3398 66,2479 0,3973 1,7897 -11,5191 2,9819 0,2193 0,3 10,0000 12,1378 -1,3397 82,7784 0,4964 -0,0296 -24,8677 4,7716 0,6166 0,4 8,6602 -12,7299 -3,6602 60,8987 0,3652 -2,6336 -27,2127 4,7420 1,1130 0,5 5,0000 -39,9426 -5,0000 -2,6205 -0,0157 -4,7994 -16,1033 2,1084 1,4782 0,6 0,0000 -56,0459 0,0000 -85,2198 -0,5110 -4,4558 22,9749 -2,6911 1,4625 0,7 0,0000 -33,0710 0,0000 -137,426 -0,8241 -1,6292 33,5584 -7,1469 0,9514 0,8 0,0000 0,4874 0,0000 -137,196 -0,8227 1,6218 31,4613 -8,7761 0,1273 0,9 0,0000 31,9487 0,0000 -87,6156 -0,5254 4,1031 18,1644 -7,1543 -0,6954 1,0 0,0000 50,1130 0,0000 -3,0512 -1,2208 3 COMPUTATIONAL RESULTS Table 4: Compare results of displacement values between methods with theoretical results. Central Difference Newmark’s constant average Newmark’s linear Analytical results t i Method acceleration method acceleration method based on theory 0,0 0,0000 0,0000 0,0000 0.0000 0,1 0,1914 0,0437 0,0300 0,0328 0,2 0,6293 0,2326 0,2193 0,2332 0,3 1,1825 0,6121 0,6166 0,6487 0,4 1,5808 1,0825 1,1130 1,1605 0,5 1,5412 1,4309 1,4782 1,5241 0,6 0,9141 1,4231 1,4625 1,4814 0,7 -0,0247 0,9622 0,9514 0,9245 0,8 -0,8968 0,1908 0,1273 0,0593 0,9 -1,3726 -0,6044 -0,6954 -0,7751 1,0 -1,2940 -1,1442 -1,2208 -1,2718 © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
9. 96 COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC RESPONSE EVALUATION OF SDOF 20.000 15.000 10.000 5.000 0.000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 -5.000 -10.000 -15.000 -20.000 Central Difference Method Newmark’s constant average acceleration method Newmark’s linear acceleration method Analytical results based on theory Figure 3: Compare results of displacement values between 3 methods. 4 CONCLUSION This study deals with three methods of calculating dynamic responses of a single-degree of freedom oscillator, i.e., central difference method (CDM) and Newmark’s method (NBM), finite element method (FEM) using recorded ground acceleration for 60seconds. The study shows that by comparing all the methods for the SDOF linear problem with the time step of 0.1 sec the Newmark’s average acceleration method is the most accurate method from the other three methods which are being compared as it gives almost similar results to the theoretical one. Therefore, the study concludes that to get optimum accuracy we can use average acceleration method when the time step is 0.1 sec. REFERENCES [1] Hughes, T. J. R., 1987, The Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ. [2] Newmark, N. M., 1954, "A Method of Computation for Structural Dynamics," Journal of the Engineering Mechanics Division, ASCE, pp. 67-94. [3] A. Hrennikoff, Solution of Problems of Elasticity by the Frame-Work Method, (1941). ASME J. Appl. Mech. 8, pp 619–pp 715. [4] Oden, J. T., 1969, "A General Theory of Finite Elements II. Applications," International Journal for Numerical Methods in Engineering, Vol. 1, pp. 247-259. [5] Abramowitz, M. and Stegun, I. A. (Eds.). "Differences." §25.1 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. New York: Dover, pp. 877-878, 1972. © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh
10. COMPARISION AND STUDY OF NUMERICAL METHODS FOR DYNAMIC 97 RESPONSE EVALUATION OF SDOF [6] Jeffreys, H. and Jeffreys, B. S. "Central Differences Formula." §9.084 in Methods of Mathematical Physics, 3rd ed. Cambridge, England: Cambridge University Press, pp. 284-286, 1988. [7] Sheppard, W. F. "Central-Difference Formula." Proc. London Math. Soc. 31, 449-488, 1899. [8] Whittaker, E. T. and Robinson, G. "Central-Difference Formula." Ch. 3 in The Calculus of Observations: A Treatise on Numerical Mathematics, 4th ed. New York: Dover, pp. 35-52, 1967. [9] Hashamdar, H., Ibrahim, Z., & Jameel, M. (2011). Finite element analysis of nonlinear structures with Newmark method, 6(6), 1395–1403. NGHIÊN CỨU VÀ SO SÁNH CÁC PHƯƠNG PHÁP SỐ TRONG PHÂN TÍCH BÀI TOÁN ĐỘNG LỰC HỌC HỆ MỘT BẬC TỰ DO THAI PHUONG TRUC Khoa Kỹ thuật Xây dựng, Đại học Công nghiệp thành phố Hồ Chí Minh; thaiphuongtruc@iuh.edu.vn Tóm tắt. Ngày nay với việc phát triển ngày càng mạnh mẽ của khoa học máy tính. Việc ứng dụng phương pháp phần tử hữu hạn để giải các bài toán kỹ thuật ngày càng trở nên phổ biến. Tuy nhiên, việc giải và phân tích số chỉ cho kết quả gần đúng với lời giải giải tích. Hiện có nhiều phương pháp số cho các kỹ sư lựa chọn, mỗi phương pháp có những ưu và nhược điểm riêng và mức độ chính xác cũng khác nhau. Trong bài báo này, tác giả so sánh kết quả của ba phương pháp: Phương pháp sai phân trung tâm, phương pháp gia tốc trung bình không đổi, phương pháp gia tốc tuyến tính với lời giải lý thuyết chính xác. Thông qua kết quả thu được, các kỹ sư sẽ có được cái nhìn tổng quát hơn về các phương pháp và lựa chọn cho mình một công cụ tính toán hợp lý để giải quyết các vấn đề tương tự. Từ khóa. Phương pháp số, phương pháp phần tử hữu hạn, phương pháp sai phân trung tâm, phương pháp gia tốc trung bình không đổi, phương pháp gia tốc tuyến tính. Ngày nhận bài: 21/11/2019 Ngày chấp nhận đăng: 19/03/2020 © 2020 Trường Đại học Công nghiệp Thành phố Hồ Chí Minh