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

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 
 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 
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 
 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 
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 
 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 
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 
 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 
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 
 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