Phân tích động học của dầm timoshenko chịu nhiều tải trọng di động và ảnh hưởng của nhiệt độ bằng phương pháp phần tử hữu hạn

SCIENCE TECHNOLOGY 
Số 51.2019 ● Tạp chí KHOA HỌC & CÔNG NGHỆ 39
PHÂN TÍCH ĐỘNG HỌC CỦA DẦM TIMOSHENKO 
CHỊU NHIỀU TẢI TRỌNG DI ĐỘNG VÀ ẢNH HƯỞNG 
CỦA NHIỆT ĐỘ BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 
DYNAMIC ANALYSIS OF TIMOSHENKO BEAMS UNDER MOVING MANY LOADS AND THE EFFECT 
OF TEMPERATURE BY THE FINITE ELEMENT METHOD (FEM) 
Nguyễn Văn Luật*, 
Trần Thị Thu Thủy, Nguyễn Thị Thu Hường 
TÓM TẮT 
Bài báo nghiên cứu ứng xử động học của dầm Timoshenko dựa trên phương 
pháp phần tử hữu hạn, trong đó xây dựng công thức phần tử hữu hạn với ba 
chuyển vị nút cho dầm chịu nhiều tải trọng di động và ảnh hưởng của nhiệt độ. 
Kết quả tính đưa ra được ứng xử động học của dầm chịu nhiều tải trọng điều hòa 
di động và sự thay đổi của nhiệt độ, trong đó có so sánh với dầm Bernoulli. 
Nghiên cứu dầm Timoshenko bằng phương pháp phần tử hữu hạn đã có nhiều 
tác giả quan tâm nghiên cứu, tuy nhiên điểm khác biệt của bài báo là ở cách xây 
dựng công thức phần tử hữu hạn trực tiếp từ nguyên lý năng lượng cực tiểu với 
cách lựa chọn hàm dạng Kosmatka dựa trên đặc điểm của dầm. 
Từ khóa: Dầm Timoshenko, tải trọng động, biến dạng nhiệt, phương pháp 
phần tử hữu hạn (FEM). 
ABSTRACT 
The article studies the dynamic behavior of Timoshenko beams based on 
finite element method which builds the element formula FEM for Timmoshenko 
beam under moving harmonic load with three nodal displacement. The 
calculated results show dynamic behavior of beam under action of moving 
harmonic many loads, which are compared with Bernoulli beams. Timoshenko 
beams which used FEM method has been researched by many authors, but the 
difference of the article is the way to construct the FEM directly from the 
principle of minimum energy with the selection of Kosmatka shape function
based on the characteristics of the beam. 
Keywords: Timoshenko beam, moving load, temperature deformation, finite 
element methods (FEM). 
Khoa Cơ khí, Trường Đại học Công nghiệp Hà Nội 
*Email: luatnv1980@gmail.com 
Ngày nhận bài: 30/8/2018 
Ngày nhận bài sửa sau phản biện: 06/3/2019 
Ngày chấp nhận đăng: 25/4/2019 
1. MỞ ĐẦU 
Trong nghiên cứu [10] tác giả đã trình bày phương pháp 
phần tử hữu hạn (PTHH) ứng dụng với mô hình dầm dầm 
Timoshenko chịu một tải trọng di dộng, trong nghiên cứu 
này nhóm tác giả tiếp tục nghiên cứu với nhiều tải trọng di 
động đồng thời và có sự thay đổi của nhiệt độ. Điều này 
cho phép nghiên cứu có thể ứng dụng mô phỏng động 
học với các mô hình bài toán đa dạng hơn. Trong thực tế 
ảnh hưởng của nhiệt độ tới biến dạng của vật liệu là đáng 
kể, nhất là trong các chi tiết máy, kết cấu với môi trường 
làm việc ở nhiệt độ cao. Ngoài ra các kết cấu dầm thường 
chịu nhiều tải trọng di động đồng thời như các kết cấu cầu, 
đường ray dưới tác động của các phương tiện giao thông. 
Mô hình dầm Timoshenko bằng các phương pháp số đã 
được nhiều nhà khoa học trong và ngoài nước quan tâm 
nghiên cứu trong đó có sử dụng phương pháp PTHH 
[2,3,4], như ứng xử động học với dầm Bernoulli đã được 
công bố trong [7], với mô hình Timoshenko cũng được 
phân tích trong hầu hết các bài toán về dầm [6] trong đó sử 
dụng phương trình vi phân cân bằng cho một phần tử dầm 
cho ra các hàm dạng có khả năng hội tụ nhanh. Tuy nhiên 
điểm khác biệt trong cách tiếp cận của bài báo là các ma 
trận độ cứng, ma trận khối lượng, véc tơ lực nút và lực nhiệt 
được xây dựng trực tiếp dựa trên nguyên lý năng lượng cực 
tiểu, từ phiếm hàm năng lượng thiết lập được qua chuyển 
vị nút với các hàm dạng được chọn dựa trên các hàm dạng 
của Kosmatka theo đặc trưng của dầm Timoshenko, trong 
đó chứa tham số biến dạng trượt. Thuật toán số để tính 
toán cho bài toán động dựa trên thuật toán lặp Newmark. 
Các kết quả tính được so sánh với mô hình dầm Bernulli để 
thấy được ứng xử động học thay đổi giữa hai mô hình dầm. 
2. XÂY DỰNG CÔNG THỨC PHẦN TỬ HỮU HẠN 
Với sự xuất hiện của biến dạng trượt mô hình dầm 
Timoshenko giả thuyết rằng thiết diện ngang phẳng sau 
biến dạng vẫn phẳng nhưng không còn trực giao với lớp 
trung hòa, biến dạng dọc trục và biến dạng trượt cho bởi [5]: 
x
u θε z
x x
 
= 
  
ωψ θ
x

= 
 
(1) 
Xem xét phần tử dầm mặt cắt hình chữ nhật với mỗi nút 
có 3 chuyển vị nút u, ω, θ lần lượt là chuyển vị dọc trục theo 
phương x, chuyển vị theo phương z, góc xoay của mặt cắt 
 CÔNG NGHỆ 
 Tạp chí KHOA HỌC & CÔNG NGHỆ ● Số 51.2019 40
KHOA HỌC
ngang, Do có biến dạng trượt nên u, ω, θ độc lập với nhau 
và không phụ thuộc vào z. 
Hình 1. Mô hình dầm Timoshenko 
Các thành phần ứng suất theo định luật Húc được xác 
định bởi 
x xσ Eε= , ( )
ωτ κG θ
x

= 
 
2) 
Khi chịu tác động của nhiệt độ: Giả sử nhiệt độ của dầm 
tăng thêm T độ, sự thay đổi nhiệt độ ảnh hưởng đến chiều 
ngang của dầm rất nhỏ nên có thể bỏ qua biến dạng này. 
Khi nhiệt độ thay đổi nhỏ và không ảnh hưởng đến sự thay 
đổi hệ số giãn nở nhiệt thì chiều dài của dầm thay đổi một 
khoảng δT (δT có thể dương hoặc âm khi nhiệt độ tăng hoặc 
giảm tương ứng) δT = αTL. Trong đó, α là hệ số giãn nở 
nhiệt của vật liệu là một trong những đặc trưng cơ học của 
vật liệu, với thép α = 12.10-6 (1/0C), với aluminum α = 23.10-6 
(1/0C), L là chiều dài của phần tử dầm. 
Biến dạng dọc trục do nhiệt là: εT = αT (3) 
Khi đó ứng suất tại điểm nào đó của phần tử dầm được 
xác định bởi: 
σ = E(εx - εT) (4) 
Hàm năng lượng biến dạng của một phần tử dầm có 
dạng [3]: 
[ ]
V
1U σε τψ dV
2
= + = [( ) ( ) ]x T x T
V
1 ε ε E ε ε τψ dV
2
 + 
( )2 2x T x T
V V
1 1E ε 2ε ε ε dV τψdV
2 2
= + + 
2 2
x T x T
V V V V
1 1 1Eε dV τψdV Eε ε dV Eε dV
2 2 2
= + + 
(5) 
Thay (1), (2) vào biểu thức năng lượng (5) thu được: 
/
, , , ,
/
/
, , ,
/
( )
( ) ( )
L h 2
2 2 2
x x x x
0 h 2
L L h 2
2
x x x T
0 0 h 2
T
V
1U Eb u 2zu θ z θ dxdz
2
κ GA θ w dx Eb u zθ ε dxdz
2
1 Eε dV
2
= +
+ 
+
(6) 
Do u,x, θ,x 
không phụ thuộc vào z nên 
/
, ,
/
h 2
x x
h 2
zu θ dV 0
= 
từ (6) nhận được biểu thức năng lượng: 
, , ,
,
( ) ( )
L L
2 2 2
x x x
0 0
L
x T T
0 V
1 κU EAu EIθ dx GA θ w dx
2 2
1EAu ε dx Eε dV
2
= + + 
 +
(7) 
Trong đó, 
/
/
, . ,
h 2
2
h 2
I z dz A b h dV Adx
= = = , κ: hệ số điều 
chỉnh cho độ vênh của thiết diện ngang, với mặt cắt hình 
chữ nhật theo tài liệu [2] thì κ = 5/6. 
Nội suy các chuyển vị của dầm qua các chuyển vị nút d: 
. ,uu=N d . ,ωω =N d .θθ=N d (8) 
đặt vào biểu thức năng lượng thu được 
, , , ,
, ,( ) ( )
L L
T T T T
u x u x θ x θ x
0 0
L
T T
ω x θ ω x θ
0
1 1U EA dx EI dx
2 2
κ GA dx
2
= +
+ 
d N N d d N N d
d N N N N d
,
L
T
u x T T
0 V
1EA ε dx Eε dV
2
 + d N
(9) 
Giả sử dầm chịu tác động của nhiều tải trọng điều hòa 
i iF P cos t=  (i = 1,..,Nf) chuyển động với cùng vận tốc v 
trên dầm. 
Sử dụng nguyên lý năng lượng cực tiểu dẫn đến 
phương trình cân bằng giữa công của nội lực với công của 
ngoại lực và lực quán tính trên chuyển dịch khả dĩ động 
học thu được phương trình phần tử hữu hạn cho kết cấu: 
. . T= +
..
M D+ K D F F (10) 
Trong đó, M là ma trận khối lượng kết cấu, mật độ khối 
lượng ρ. 
1
( )
eN
i
i=
= M m , 
L
T
0
ρA dx= m N N , ( )
eN
i
i 1=
= D d
(11) 
K là ma trận độ cứng kết cấu, FT là véc tơ lực nhiệt, Ne là số phần tử 
( )
eN
i
i 1=
= K k , b s= +K K K , ( )
eN
T T i
i 1=
= F f , 
eN
T
i ω
i 1
P .cos( t).
=
= F N (12) 
Với 
SCIENCE TECHNOLOGY 
Số 51.2019 ● Tạp chí KHOA HỌC & CÔNG NGHỆ 41
, , , ,
L L
T T
b u x u x θ x θ x
0 0
EA dx EI dx= + K N N N N 
, ,( ) ( )
L
T
s ω x θ ω x θ
0
κ GA dx= K N N N N
(13) 
L
T u,x T
0
EA dx=  F N 
Kb, Ks có thể gọi là ma trận độ cứng uốn, độ cứng trượt 
tương ứng của phần tử. 
Do u, ω, θ độc lập trong lý thuyết dầm Timoshenko nên 
trong phân tích PTHH thường hay sử dụng hàm dạng nội 
suy tuyến tính cho các chuyển vị nút này, nhưng khi áp 
dụng tính toán lại hay xảy ra hiện tượng nghẽn trượt do các 
hàm dạng không đặc trưng đúng cho cấu hình biến dạng 
của dầm. Để tránh hiện tượng này thì trong tính toán PTHH 
thường hay sử dụng tính tích phân hàm năng lượng dựa 
trên phép cầu phương Gauss [3]. Cách tiếp cận của bài báo 
để tránh hiện tượng nghẽn trượt, ngoài ra theo (1) giữa 
chuyển vị ngang và góc quay vẫn có ảnh hưởng lẫn nhau, 
có thể sử dụng hàm dạng tuyến tính cho chuyển vị dọc trục 
u và góc quay θ, hàm dạng Kosmatka [6] chứa tham số biến 
dạng trượt  cho chuyển vị ngang ω: 
u
ω
θ
u
ω
θ
= 
N
N d
N
(14) 
Trong đó: 
 T1 1 1 2 2 2u ω θ u ω θ=d 
u
L x x0 0 0 0
L L
= 
N 
1 2 3 4ω ω ω ω ω0 N N 0 N N = N 
θ
L x x0 0 0 0
L L
= 
N
(15) 
( )
1
3 2
3 2
1 x x xN 2 3 1
1 LL L
 

= + +
+
ω , 
( ( ) ( ) )
2
3 2
ω 3 2
L x x xN 2 1
1 2 2 LL L
 

= + + +
+
( )
2
3 2
ω 3 2
1 x x xN 2 3
1 LL L


= 
+
, 
( ( ) )
4
3 2
ω 3 2
L x x xN 1
1 2 2 LL L
 

= 
+
Đặt các hàm dạng trong (12) vào (11) thu được các ma 
trận độ cứng, véc tơ lực nhiệt của phần tử dầm 
Timoshenko. 
 TT EAαT 1 0 0 1 0 0= f 
b
EA EA0 0 0 0
L L
0 0 0 0 0 0
EI EI0 0 0 0
L L
EA EA0 0 0 0
L L
0 0 0 0 0 0
EI EI0 0 0 0
L L
 =
K (16) 
( )2
s
0 0 0 0 0 0
1 1 1 10 0
L 2 L 2
1 L 1 L6 0 0κGA 2
2 4 2 45
0 0 0 0 0 01
1 1 1 10 0
L 2 L 2
1 L 1 L0 0
2 4 2 4
 

 + +
 =
 +
K , 
2
12EI
κAGL
 =
(17) 
3. ÁP DỤNG VÀ SO SÁNH 
Áp dụng các kết quả xây dựng được ở trên với ma trận 
độ cứng, ma trận khối lượng, véctơ lực nút và lực nhiệt để 
phân tích ứng xử động học cho mô hình dầm Timoshenko 
chịu nhiều tải trọng di động và sự thay đổi của nhiệt độ 
trong đó có so sánh với dầm Bernoulli. Sử dụng phương 
pháp gia tốc trung bình trong họ các phương pháp 
Newmark [3] để xây dựng thuật toán số, với bước thời gian 
đảm bảo thuật toán ổn định thỏa mãn Δt ≤ 2/ωmax, ωmax là tần số riêng lớn nhất. Chương trình tính được viết trên 
phần mềm Matlab. 
Dầm Bernoulli chịu tải trọng di động với 3 chuyển vị nút 
thì ma trận độ cứng, ma trận khối lượng được kết hợp từ 
phần tử dầm Bernoulli 2 nút sử dụng hàm dạng Hermite và 
phần tử thanh được cho dưới dạng [3]: 
3 2 3 2
2 2
3 2 3 2
2 2
EA EA0 0 0 0
L L
12EI 6EI 12EI 6EI0 0
L L L L
6EI 4EI 6EI 2EI0 0
L LL L
EA EA0 0 0 0
L L
12EI 6EI 12EI 6EI0 0
L L L L
6EI 2EI 6EI 4EI0 0
L LL L
= 
K
(18) 
 CÔNG NGHỆ 
 Tạp chí KHOA HỌC & CÔNG NGHỆ ● Số 51.2019 42
KHOA HỌC
/ /
/ /
2 2
2 2
2 6 0 0 1 6 0 0
156 22L 54 13L0 0
420 420 420 420
22L 4L 13L 3L0 0
420 420 420 420ρAL
1 6 0 0 2 6 0 0
54 13L 156 22L0 0
420 420 420 420
13L 3L 22L 4L0 0
420 420 420 420
= 
m
(19) 
Xét mô hình dầm Timoshenko chịu số tải trọng Nf = 5 tải 
trọng điều hòa di động cách đều nhau, có cùng tần số và 
vận tốc v (hình 2): Fi = Picos(Ωt), i = 15. Nhiệt độ thay đổi 
là ΔT = 200C, α = 12.10-6 (1/0C). Dầm thép có chiều dài 
L = 20m, mặt cắt ngang hình chữ nhật có các thông số kích 
thước, vật liệu và tải trọng: h = 2b = 0,2m; E = 2.1011N/m2, 
G = 8.1010N/m2, ρ = 7860 (kg/m3), κ = 5/6, P1 = 2000N, 
P2 = 3000N, P3 = 4000N, P4 = 5000N, P5 = 6000N. 
Các kết quả tính trong bảng 1 cho thấy, biến dạng nhiệt 
của dầm phụ thuộc vào nhiệt độ và vận tốc của tải trọng 
nhưng không phụ thuộc vào tần số của tải trọng. Sự phụ 
thuộc của biến dạng nhiệt vào vận tốc hoàn toàn không 
theo quy luật tuyến tính. So sánh ứng xử động học của mô 
hình dầm Timoshenko và dầm Bernoulli khi vận tốc và tần số 
của tải trọng động thay đổi được thể hiện trên các hình 3, 4, 
5 (để thấy rõ ảnh hưởng của sự thay đổi nhiệt độ đến biến 
dạng của dầm thì dịch chuyển theo trục của dầm do nhiệt 
độ được nhân thêm 103m). Từ các kết quả thu được cho thấy 
ảnh hưởng lẫn nhau giữa vận tốc của tải trọng, tần số lực 
kích động và sự thay đổi nhiệt độ đến ứng xử động học và 
biến dạng của dầm. Cụ thể ở hình 3 khi tần số lực kích động 
nhỏ thì ứng xử động học của hai mô hình dầm không khác 
nhau nhiều nhưng vẫn có sự khác biệt khi vận tốc của các tải 
trọng nhỏ (v = 20m/s). Ở các hình 4 và 5 khi tăng tần số lực 
kích động thì ứng xử động học của hai mô hình khác nhau rõ 
rệt khi vận tốc tải trọng nhỏ (v = 20m/s), khi vận tốc tải trọng 
lớn (v = 80m/s) thì ứng xử động học của hai mô hinh dầm 
cũng tương tự nhau. Như vậy việc sử dụng mô hình dầm nào 
tùy thuộc vào từng bài toán cụ thể như bài toán dao động, 
tải trọng di động có tần số lực kích động nhỏ hoặc vận tốc tải 
trọng lớn có thể sử dụng mô hình dầm Bernoulli, khi lực kích 
động có tần số lớn hoặc tải trọng di chuyển với vận tốc nhỏ 
thì cần phải tính tới ảnh hưởng biến dạng trượt nên mô hình 
dầm Timoshenko sẽ phù hợp hơn. Qua các hình 3, 4, 5 có thể 
thấy biên độ dao động của dầm phụ thuộc vào tần số của 
lực kích động khá rõ, tần số tải trọng càng tăng thì biên độ 
càng giảm. 
Hình 2. Mô hình dầm chịu nhiều tải trọng di động 
Bảng 1. Sự thay đổi của nhiệt độ, vận tốc đến biến dạng nhiệt của dầm 
ΔT 
v = 20m/s v = 40m/s v = 60m/s v = 80m/s 
ΔL (m) ΔL (m) ΔL (m) ΔL (m) 
10 0,0041 0,00044 0,0035 0.00034 
20 0,0082 0,00088 0,0071 0,00069 
30 0,0122 0,0013 0,0106 0,0010 
Hình 3. Ứng xử động học dầm Timoshenko và dầm Bernoulli trong trường 
hợp Ω = 20 (rad/s) 
SCIENCE TECHNOLOGY 
Số 51.2019 ● Tạp chí KHOA HỌC & CÔNG NGHỆ 43
Hình 4. Ứng xử động học dầm Timoshenko và dầm Bernoulli trong trường 
hợp Ω = 40 (rad/s) 
Hình 5. Ứng xử động học dầm Timoshenko và dầm Bernoulli trong trường 
hợp z = 60 (rad/s) 
4. KẾT LUẬN 
Dưới tác động của nhiều tài trọng động đồng thời và sự 
thay đổi của nhiệt độ bài báo đã đưa ra được ứng xử động 
học cho mô hình dầm Timoshenko. Cụ thể xây dựng được 
thuật toán số cho bài toán dựa trên phương pháp 
Newmark, chương trình số với thuật toán và công thức 
PTHH được viết bằng ngôn ngữ Matlab để phân tích ứng xử 
động học của dầm chịu nhiều tải trọng di động, nhiệt độ. 
Kết quả nghiên cứu này có thể giúp cho việc mô phỏng 
ửng xử động học của dầm đa dạng và sát với các bài toán 
thực tế hơn. Qua các kết quả tính thấy được sự ảnh hưởng 
rõ rệt của vận tốc, tần số của các tải trọng và sự thay đổi 
nhiệt độ đến ứng xử động học và biến dạng nhiệt của dầm. 
So sánh ứng xử động học của dầm Timoshenko với mô 
hình dầm Bernoulli thì kết quả sai lệch không đáng kể khi 
tần số kích động của tải trọng nhỏ hoặc vận tốc tải trọng 
lớn và khác nhau rõ rệt khi lực kích động có tần số lớn hơn 
do ảnh hưởng của biến dạng trượt bên trong. Từ đó có thể 
giúp ích cho việc lựa chọn mô hình dầm thích hợp với từng 
bài toán thực tế. 
TÀI LIỆU THAM KHẢO 
[1]. T.J.R. Hughes, 2000. The finite element method. Linear static and dynamic 
finite element analysis. Dover publication, Inc., Mineola. 
[2]. I.H. Shames and C.L. Dym, 1985. Energy and finite element methods in 
structural mechanics. McGraw-Hill, New York. 
[3]. Daryl L. Logan, 2007. A first course in the Element Finite Method, 4th 
Edition, Thomson Canada Limited. 
[4]. Y.W. Kwon and H. Bang, 2000. Finite element method using Matlab. CRC 
Press, New York, 2nd edition. 
[5]. D. C. Pham, 2013. Essential solid mechanics. Institute of Mechanics, 
Hanoi. 
[6]. Kosmatka. J.B., 1995. An improved two-node finite element for stability 
and natural frequencies of axial-loaded Timoshenko beams, Computers and 
Structures. 57, pp. 141-149. 
[7]. Nguyen Dinh Kien, Tran Thanh Hai, 2006. Dynamic analysis of prestressed 
Bernoulli beams resting on two-parameter foundation under moving harmonic 
load. Vietnam Journal of Mechanics, 28, pp. 176-188. 
[8]. M. Simsek, T. Kocaturk, D. Akbas, 2012. Dynamic behavior of an axially 
functionally graded beam under action of a moving harmonic load. Composite 
Structures, 94, pp. 2358-2364. 
[9]. K. Rajabi, M.H. Kargarnovin, M. Gharini, 2013. Dynamic analysis of a 
functionally graded simply supported Euler - Bernoulli beam subjected to a moving 
oscillator. Acta Mechanica, 224, pp. 425-446. 
[10]. Nguyễn Văn Luật, Khuất Đức Dương, Nguyễn Thi Thu Hường, 2018. 
Phương pháp phần tử hữu hạn cho dầm Timoshenko chịu tải trọng di động. Tạp chí 
Khoa học & Công nghệ, Trường Đại học Công nghiệp Hà Nội, số 47, 10/2018. 
AUTHORS INFORMATION 
Nguyen Van Luat, Tran Thi Thu Thuy, Nguyen Thi Thu Huong 
Faculty of Mechanical Engineering, Hanoi University of Industry