Luận án Dao động của dầm FGM có lỗ rỗng vi mô trong môi trường nhiệt độ chịu tải trọng di động

ác biểu thức năng lượng phần tử. Biểu thức cho ma trận độ cứng, ma 
trận khối lượng được thiết lập chi tiết cho mô hình dựa trên lý thuyết biến dạng 
trượt bậc nhất với các hàm dạng thứ bậc. Mô hình phần tử dầm dựa trên lý thuyết 
biến dạng trượt bậc nhất với các hàm dạng chính xác và mô hình phần tử Euler-
Bernoulli cũng được đề cập. Cuối chương trình bày thuật toán số trên cơ sở phương 
pháp tích phân trực tiếp Newmark dùng để tính toán đáp ứng động lực học của dầm 
FGM có lỗ rỗng trong môi trường nhiệt độ chịu tải trọng di động. 
3.1. Véc tơ chuyển vị nút 
Phần tử dầm phẳng hai nút điển hình dùng trong phân tích dao động của dầm 
FGM được như minh họa trên hình 3.1. Trên hình vẽ, l là chiều dài phần tử, ui, wi, θi 
là chuyến vị dọc trục, chuyển vị ngang và góc quay tại nút i; uj, wj, θj là các đại 
lượng tương ứng ở nút j. 
Hình 3.1. Chuyển vị nút (a) và lực nút (b) của phần tử dầm 
Véc-tơ chuyển vị nút d cho một phần tử dầm đặc trưng (i, j) bao gồm 6 thành phần 
 { }Ti i i j j ju w u w  d (3.1) 
trong đó, chỉ số trên ‘T’ được sử dụng để chỉ chuyển vị của một véc tơ hoặc ma 
trận; Ni; Qi; Mi và Nj, Qj ; Mj tương ứng là lực dọc trục, lực cắt và mô men tại các 
nút i và j. Để thiết lập biểu thức cho ma trận độ cứng và ma trận khối lượng cho 
phần tử ta cần đưa vào các hàm nội suy cho chuyển vị dọc trục, chuyển vị ngang và 
a)
b)
uj
Mi
Qi
Nj
QjMj
i j
ui
wi i
i j
Ni
l
40 
góc quay. Khi đó, chuyển vị u0(x), w0(x) và góc quay(x) của tiết diện ngang cho 
phần tử dầm được nội suy qua các hàm dạng như sau 
 0 0, ,u wu w  N d N d N d (3.2) 
trong đó Nu, Nw và N tương ứng là các hàm nội suy (hàm dạng) cho u0(x), w0(x) và 
(x). Một số phần tử dầm với 6 bậc tự do như trong phương trình (3.2) sử dụng các 
các nội suy khác nhau dùng trong phân tích dao động của dầm FGM đã được một số 
tác giả đề nghị trong thời gian gần đây. 
3.2. Hàm nội suy thứ bậc 
Mô hình phần tử hữu hạn cho phần tử dầm FGM có thể được xây dựng trên 
cơ sở lựa chọn các hàm dạng khác nhau. Trong luận án của Lê Thị Hà [12], ma trận 
độ cứng và ma trận khối lượng được xây dựng từ các hàm dạng chính xác, nhận 
được từ lời giải phương trình tĩnh học của một phần tử dầm. Mô hình phần tử hữu 
hạn xây dựng trên các hàm dạng chính xác với 6 bậc tự do có tốc độ hội tụ cao, tuy 
nhiên các hàm dạng này phải tính toán lại mỗi khi thay đổi lưới phần tử, vì thế tốn 
thời gian tính toán. Để tránh yếu điểm này, Luận án sử dụng hàm dạng nội suy thứ 
bậc (hierachical shap functions) để nội suy các chuyển vị và góc quay của phần tử 
dầm. Cần nhấn mạnh rằng số bậc tự do của mô hình dựa trên các hàm dạng thứ bậc 
thường cao hơn 6, và không phải tất cả các bậc tự do này đều là các chuyển vị hay 
góc quay tại nút phần tử. 
 Phương pháp PTHH sử dụng các hàm dạng thứ bậc được trình bày chi tiết 
trong các tài liệu [132, 133]. Điểm chính của các đa thức thứ bậc là các đa thức bậc 
cao luôn chứa tất cả các số hạng của đa thức bậc thấp hơn. Như vậy, để nhận được 
một đa thức mới bậc cao hơn, ta chỉ cần thêm vào các đa thức cũ một vài số hạng. 
Chẳng hạn, từ phép nội suy bậc nhất cho hàm f(x) trên miền [0,l ] 
1 1 2 2 1 2
1 1
2 2
f N f N f f f
  
 (3.3) 
với 2 1x
l
 là tọa độ tự nhiên. Phép nội suy bậc 2 cho hàm f(x) có thể được viết 
dưới dạng 
 1 1 2 2 3 3f N f N f N f (3.4) 
41 
trong đó N3 là đa thức bậc hai, có thể biểu diễn dưới dạng 
 23 0 1 2N c c c  (3.5) 
Các hằng số c0, c1, c2 được tìm sao cho N3 triệt tiêu tại hai đầu nút phần tử, = 1. 
Như vậy có vô số cách chọn N3. Chẳng hạn c2 được chọn sao cho N3 = 1 ở giữa 
phần tử ( = 0) ta nhận được biểu thức cho N3 
 23 1N  (3.6) 
 Tương tự, thêm đa thức bậc ba vào vế phải của (3.4) ta nhận được phép nội 
suy bậc 3, trong đó hàm N4 có thể biểu diễn dưới dạng 
 2 34 0 1 2 3N d d d d   (3.7) 
Các hằng số d0, d1, d2, d3 được tìm sao cho N4 = 0 tại  = 1. Tương tự như trên ta 
có vô số cách chọn N4. Một trong số cách chọn là N4 = 0 tại giữa phần tử trong khi 
đạo hàm của nó bằng 1. Trong trường hợp này ta nhận được 
 24 1N   (3.8) 
Tương tự ta có thể xây dựng được các đa thức bậc cao hơn cho các hàm dạng thứ 
bậc. Các hàm dạng N1, N2, N3, N4 cho bởi phương trình (3.3), (3.6) và (3.8) được 
minh họa trên Hình 3.2a. 
Hình 3.2. a) Hàm dạng thứ bậc; (b) chi tiết về chuyển vị và góc quay 
42 
Sử dụng các hàm nội suy N1, N2, N3, N4 nêu trên, chuyển vị dọc trục u0, 
chuyển vị ngang w0 và góc xoay  được nội suy như sau 
0 1 1 2 2
1 1 2 2 3 3
0 1 1 2 2 3 3 4 4
u N u N u
N N N
w N w N w N w N w
   
 (3.9) 
trong đó u1, u2, 1,, w4 là giá trị của các biến ở các nút và ở trong phần tử. Hình 
3.2(b) minh họa giá trị của chuyển vị và góc xoay theo phương trình (3.9). 
3.3. Trường chuyển vị với ràng buộc 
Phần tử dầm dùng để phân tích động lực học của kết cấu có thể được xây 
dựng từ 9 bậc tự do như trên hình 3.2. Tuy nhiên, phần tử sẽ hiệu quả hơn khi số 
bậc tự do ít hơn. Việc này có thể nhận được bằng cách hạn chế biến dạng trượt là 
một hằng số [133]. Biến dạng trượt trong (2.25) nhận được từ phép nội suy (3.9) có 
dạng 
2
z 4 3 3 1 2
2 1 4 1 2 3
6 4 1 1
2 2
1 12 2
2
x w wl l
w w w
l
     
  
 (3.10) 
Để xz = const, ta cần có 
4 3
3 1 2
6 0
4 1 1 0
2 2
w
l
w
l

 
 (3.11) 
Từ (3.11) ta rút ra 
 3 1 2 4 3,8 6
l lw w   (3.12) 
Sử dụng công thức (3.12) ta có thể viết lại (3.9) như sau 
1 2
2
1 2 3
2 2
1 2 1 2 3
1 11 1
2 2
1 11 1 1
2 2
1 11 1 1 1
2 2 8 6
u u u
l lw w w
 
      
       
 (3.13) 
43 
Và biến dạng trượt bây giờ có dạng 
 2 1 1 2 31 1 22 3xz w wl    (3.14) 
Phần tử dầm trong luận án này được xây dựng từ trường chuyển vị theo công 
thức (3.13) và biến dạng trượt theo công thức (3.14). Véc tơ chuyển vị nút cho một 
phần tử có 7 bậc tự do là 
 1 1 1 3 2 2 2 Tu w u w   d (3.15) 
Sử dụng công thức (3.15), trường chuyển vị và góc xoay (3.13) có thể được 
viết dưới dạng ma trận 
 0 0, ,u wu w  N d N d N d (3.16) 
trong đó 
  1 2 1 3 2
1 3 4 2 3
0 0 0 0 0 , 0 0 0 0 ,
0 0
8 6 8
TT
u
T
w
N N N N N
l l lN N N N N
 
  
 
N N
N
 (3.17) 
tương ứng là các ma trận hàm dạng của u0, θ và w0. 
3.4. Ma trận độ cứng phần tử 
Ma trận độ cứng của dầm được xây dựng dựa trên biểu thức năng lượng biến 
dạng. Sử dụng phép nội suy (3.16), (3.17), và với lưu ý (.),x = (.),ξ 2/l và dx = l/2dξ, 
ta có thể viết biểu thức cho năng lượng biến dạng đàn hồi cho một phần tử dầm, Ue, 
dưới dạng sau 
1
22 2
11 0, 12 0, , 22 , 33 0,
1
1
11 12 , 22
1 ,
, , , ,
3
,
3 ,( ( )
1 2
21
)
1
2
e
T
T T T
u u
w
u
w
T
T T
U A u A u A A w d
l
A A A
d
l A
 

    
      
  
    


N N N N N N
N
d
N
d
N
d
N
kd
 (3.18) 
trong đó 
 uu u   k k k k k (3.19) 
với k là ma trận độ cứng phần tử, và 
44 
11
12
, 3
1
, ,
1
1
, ,
1
1
, ,
1
1
3 ,
1
2
2
2
2 2
2
T
uu u u
T
u u
T
w
T T
w
d
l
d
l
d
l
l d
l
A
A
l
A
 
   
    
    



 
k N N
k N N
k N N
k N N N N
 (3.20) 
tương ứng là các ma trận độ cứng sinh ra do biến dạng dọc trục, tương hỗ giữa dọc 
trục và uốn, biến dạng uốn và biến dạng trượt. Các ma trận này có dạng toán học 
đơn giản và có thể viết dưới dạng tường minh dưới đây 
 11 12
0 1 0 0 1
1 0 0 1 0
1 1
, 0 0 0 0 0
1 1
0 1 0 0 1
1 0 0 1 0
uu u
A A
l l
kk 
2 2 2
2 2
3322 2
2 2 2
21 1
2 3 2
1 0 1 2 4 3 2 4
16 2 4 20 0 ,
3 3 3 9 3 3
1 0 1 21 1
2 3 2
2 4 3 2 4
l l l
l l l l l
AA l l l lll l
l l l
l l l l l
 

k k (3.21) 
Để tương thích trong tính toán ma trận dưới đây, các ma trận độ cứng trong phương 
trình (3.21) cần được mở rộng thành ma trận có kích thước (7x7) bằng cách thêm 
vào các hàng và cột với các hệ số bằng 0 tương ứng với các hệ số 0 trong ma trận 
các hàm dạng (3.17) 
3.5. Ma trận độ cứng do ứng suất nhiệt ban đầu 
Ma trận độ cứng phần tử do ứng suất nhiệt ban đầu được xây dựng dựa trên 
năng lượng biến dạng do tăng nhiệt độ phương trình (2.30) cho một phần tử có dạng 
 , ,
1 1
2
0,
1 1
1 1 1
2w w
T T T
Te T T TU N w N dl l l
d   
 d N N d d k d (3.22) 
45 
với kT là ma trận độ cứng do tăng nhiệt độ 
1
, ,
1
2 T
T w T wN dl  

 k N N (3.23) 
Dạng tường minh của kT có dạng đơn giản sau 
2 2
2
2 2
1 0 0 1 0
1 10 0 0
12 12
40 0 0 0
45
1 0 0 1 0
1 10 0 0
12 12
T
T
l l
N
l
l
l l
k (3.24) 
3.6. Ma trận khối lượng phần tử 
Ma trận khối lượng được xây dựng trên các hàm nội suy cho trường chuyển 
vị. Vì vậy, ta có thể viết biểu thức động năng cho phần tử dầm công thức (2.32) 
dưới dạng sau 
1
2 2 2
11 0 12 0 22
1
1
11 11 12 22
e
1
0( ) 24
( 2 )
4
1( )
2
1
2
T T T T T
u u w w u
T T
uu ww u
l I u I u I d
l I I
w
I I d 
 


  

 
 
 d N N N N N N N N d
d m m m m d d md
  
 
  



 (3.25) 
trong đó m = muu + mww + mu + m là ma trận khối lượng của phần tử dầm 
1
11
1
1
11
1
1
12
1
1
22
1
2
2
2
2
T
uu u u
T
ww w w
T
u u
T
l I
l I
l I
d
d
d
l dI

 





m N N
m N N
m N N
m N N
 (3.26) 
46 
tương ứng là các ma trận khối lượng sinh ra do chuyển vị theo phương: dọc trục; 
phương ngang; tương hỗ giữa chuyển vị dọc trục và sự quay của tiết diện ngang; sự 
quay của tiết diện ngang. Thay thế các hàm dạng vào (3.26) và tính tích phân ta 
nhận được dạng tường minh cho ma trận khối lượng như sau 
2 2
2
22 11
2 2
w
1 1 11
8 30 2 8
1 01 1
8 40 8 402
8 21 1 , 0 03 5 3 30 315 30
1 11 1 12 2 8 30 8
0
8 40 8 40
w
l
l l l l
lI lI l l l
l l l
l l l l

m m 
 11 12
10 1 1 0
2
11 0 0 011 2
2 , 1 0 0 1 0
13 31 10 1 0 12
2
1 0 0 1 0
2
uu u
lI lI

m m (3.27) 
Các ma trận trong (3.24) và (3.27) cũng cần được mở rộng thành các ma trận có 
kích thước (7x7) để bảo đảm sự tương thích trong tính toán ma trận. 
3.7. Phần tử dựa trên các hàm nội suy chính xác 
Để so sánh mô hình phần tử dầm dựa trên lý thuyết biến dạng trượt bậc nhất 
và các hàm dạng thứ bậc được xây dựng ở trên. Luận án cũng tiến hành xây dựng 
phần tử dầm dựa trên các hàm dạng chính xác. 
Kosmatka [65] xây dựng các hàm nội suy chính xác cho phần tử dầm 
Timoshenko thuần nhất bằng cách giải phương trình vi phân cân bằng cho một phần 
tử và chỉ ra rằng các hàm nội suy này có nhiều ưu điểm, đặc biệt là tốc độ hội tụ 
nhanh. Dựa trên ý tưởng này, luận án tiến sỹ của Lê Thị Hà [12] đã tiến hành xây 
dựng hàm nội suy chính xác cho phần tử dầm Timoshenko làm từ vật liệu FGM có 
47 
cơ tính biến đổi ngang. Các hàm dạng chính xác cho phần tử dầm FGM dựa trên lý 
thuyết biến dạng trượt có dạng sau 
 
 
 
1 2 3 4 5 6
1 2 3 4 5 6
1 2 3 4 5 6
N
N
N
u u u u u u u
w w w w w w w
N N N N N N
N N N N N N
N N N N N N      
 (3.28) 
2
2
1 2
2
2
3 4
2 2
2 2
5 6
6 ( )
,
(1 )
3 ( )
,
(1 )
6 ( ) 3 ( )
,
(1 ) (1 )
a
u u
a
u u
a a
u u
x x
l x l lN N
l l
x x
xl lN N
l l
x x x x
l l l lN N
l l


 
 (3.29) 
1 4
3 2
2
3 2
3
3 2
5
3
6
0
1 2 3 1
1
2 1
1 2 2
1 2 3
1
1
1 2
w w
w
w
w
w
N N
x x xN
l l l
l x x xN
l l l
x x xN
l l l
l xN
l
 

 





2
2
x x
l l
 
 (3.30) 
1 4
2
2
2
3
2
5
2
6
0
6
1
3 4 1
1
6
1
3 2
1
N N
x xN
l l l
l x xN
l l
x xN
l l l
l x xN
l l
 





 




 (3.31) 
48 
Trong các phương trình (3.30) và (3.31), 12
11
a
A
A
 và 222
33
12A
l A


 là tham số biến 
dạng trượt, phụ thuộc vào tính chất vật liệu dầm. Chú ý rằng các hàm dạng cho 
chuyển vị ngang w0(x) và góc xoay θ(x) theo công thức (3.30) và (3.31) có dạng 
giống với hàm nội suy do Kosmatka xây dựng trong [65] ngoại trừ cách định nghĩa 
tham số biến dạng trượt  . Nếu lấy 2
12EI
l GA


 thì các hàm dạng này quay về hàm 
Kosmatka. Biểu thức cho ma trận độ cứng và ma trận khối lượng cho một phần tử 
vẫn có dạng (3.20) và (3.26) nhưng với các hàm dạng cho bởi (3.29)-(3.31). 
3.8. Phần tử dầm Euler-Bernoulli 
Phần tử dầm dựa trên lý thuyết dầm Euler-Bernoulli có thể xây dựng từ 
trường nội suy tuyến tính cho u0 và đa thức Hermit cho w0. Ma trận độ cứng và ma 
trận khối lượng cho phần tử dầm Euler-Bernoulli có dạng dưới đây 
Ma trận độ cứng phần tử 
Sử dụng phép nội suy, biểu thức năng lượng biến dạng đàn hồi và biến dạng do tăng 
nhiệt độ cho một phần tử dầm Euler-Becnoulli, được xác định như sau 
2 2
11 0, 12 0, 0, 22 0,
0
11 12 22, , , , , ,
0
1 ( 2 )
2
1 ( 2 )
2
1 1( )
2 2
e x x xx xx
l
T T T T
u u u w w wx x x xx xx xx
T T
uu uw ww
l
U A u A u w A w dx
A A A dx
 
 
 d N N N N N N d
d k k k d d k d
 (3.32) 
trong đó 
 uu uw ww k k k k (3.33) 
và 
, 11 ,
0
, 12 ,
0
, 22 ,
0
2
l
T
uu u x u x
l
T
uw u x w xx
l
T
ww w xx w xx
A dx
A dx
A dx
k N N
k N N
k N N
 (3.34) 
49 
tương ứng là ma trận độ cứng: dọc trục, tương hỗ giữa dọc trục - chống uốn và độ 
cứng chống uốn. 
Ma trận độ cứng do ứng suất nhiệt ban đầu được xác định giống với phần tử 
dầm Timoshenko được xác định qua phương trình (3.23). 
Cũng giống như đối với phần tử dầm dựa trên lý tuyết biến dạng trượt bậc 
nhất, ma trận khối lượng của phần tử dầm Euler-Becnoulli được xây dựng trên các 
hàm nội suy cho trường chuyển vị. Vì vậy, ta có thể viết biểu thức động năng cho 
phần tử dầm 
2 2 2
11 12 0, 22 0,
0
11 11 12 , ,
0 0 0
22 ,
0
1 ( ( ) 2 )
2
1 ( 2 )
2
1 1( )
2 2
x x
l
T T T T T
u u w w u w x w x w x
T T
u
l
u w u
e
w
I u w I u w I w dx
I I I I dx
 
 
 
 d N N N N N N N N d
d m m m m d d md
    
 
   

 (3.35) 
trong đó m = muu + mww + mu + m là ma trận khối lượng của phần tử dầm 
110
110
12 ,0
, 22 ,0
l T
uu u u
l T
ww w w
l T
u u w x
l T
w x w x
I dx
I dx
I dx
I dx


m N N
m N N
m N N
m N N
 (3.36) 
tương ứng là các ma trận khối lượng sinh ra do chuyển động theo phương dọc trục, 
phương ngang, tương hỗ giữa chuyển vị dọc trục và sự quay của tiết diện ngang, sự 
quay của tiết diện ngang. Giống như phần tử TBHi, các ma trận độ cứng và ma trận 
khối lượng trong các phương trình (3.34) và (3.36) cần được bổ sung các hệ số bằng 
0 tạo thành ma trận có kích thước (6x6) để đảm bảo sự tương thích cho việc cộng 
các ma trận trong các phương trình (3.33) và (3.35). 
3.9. Phương trình chuyển động rời rạc 
Các biểu thức cho ma trận độ cứng và ma trận khối lượng của phần tử được 
nối ghép để tạo thành các ma trận độ cứng và ma trận khối lượng tổng thể cho dầm. 
Bỏ qua ảnh hưởng cản của vật liệu dầm, phương trình chuyển động cho dầm FGM 
50 
có lỗ rỗng vi mô chịu tải trọng di động có thể viết dưới dạng ngôn ngữ phần tử hữu 
hạn như sau 
 B
ex
T( ) MD K K D F (3.37) 
trong đó 
nE nE
1 1
,
i i 
  D d D d  (3.38) 
tương ứng là các véc-tơ chuyển vị và gia tốc tổng thể tại các điểm nút của dầm 
nE
e
nE nE nE
B T
1 1 1
x
T
1
,, ,
i i i i 
    K k K k M m F f (3.39) 
tương ứng là các ma trận độ cứng, ma trận khối lượng và véc - tơ tải trọng nút tổng 
thể. Trong các phương trình (3.38) và (3.39), nE là tổng số phần tử dầm được rời 
rạc và ký hiệu
nE
1i 
 được hiểu theo nghĩa ghép nối các ma trận và véc-tơ lực nút phần 
tử thành các ma trận độ cứng, ma trận khối lượng và véc-tơ lực nút tổng thể theo 
phương pháp phần tử hữu hạn. 
3.10. Thuật toán Newmark 
3.10.1. Họ các phương pháp Newmark 
 Phương pháp tích phân trực tiếp thay các đạo hàm riêng trong phương trình 
(3.37) bằng các sai phân hữu hạn, tức là thay thế ,D D bằng sai phân của chuyển vị 
nút D tại một số thời điểm khác nhau. Nhiều phương pháp tích phân trực tiếp khác 
nhau được trình bày trong các sách chuyên khảo và việc lựa một phương pháp phụ 
thuộc vào bài toán cũng như kinh nghiệm của người phân tích. 
Ý tưởng trung tâm của phương pháp tích phân trực tiếp là chia tổng thời gian 
t ra thành các phần nhỏ, gọi là bước thời gian. Phương trình chuyển động được 
viết tại một thời điểm có dạng sau 
 exi i i MD KD F (3.40) 
trong đó chỉ số dưới i dùng để chỉ thời điểm i t với t là bước thời gian. Xuất phát 
từ thời điểm ban đầu i = 0 (với các chuyển vị và vận tốc ban đầu đã biết, thường là 
bằng không) ta sẽ xác định được chuyển vị, vận tốc và gia tốc nút tại thời điểm tiếp 
51 
theo (i+1) t. Tập hợp các giá trị của chuyển vị, vận tốc và gia tốc tại tất cả các điểm 
thời gian cho ta bức tranh đáp ứng động lực học kết cấu của dầm theo thời gian. 
Phương pháp tích phân trực tiếp có thể chia làm hai nhóm: phương pháp tích 
phân trực tiếp hiển và phương pháp tích phân trực tiếp ẩn. Trong phương pháp tích 
phân trực tiếp hiển, véc-tơ chuyển vị nút tại thời điểm mới (i+1) t được xác định 
hoàn toàn qua các thông tin đã biết ở các thời điểm trước đó, tức là chuyển vị, vận 
tốc và gia tốc tại các thời điểm i t, (i-1) t 
 1 1 1 1( , , , , , ...)i i i i i i if D D D D D D D    (3.41) 
Chuyển vị tại thời điểm mới (i+1) t trong phương pháp tích phân trực tiếp ẩn 
không chỉ được xác định qua chuyển vị, vận tốc và gia tốc ở các thời điểm trước đó 
mà còn qua vận tốc và gia tốc tại chính thời điểm hiện tại 
 1 1 1( , , , , )i i i i i if D D D D D D    (3.42) 
 Bởi vì các véc-tơ vận tốc nút 1i D và gia tốc 1i D trong phương trình (3.42) 
là các đại lượng chưa biết, vì thế trong nhiều trường hợp thuật toán lặp cần được sử 
dụng trong phương pháp tích phân trực tiếp ẩn [134]. 
Các phương pháp tích phân trực tiếp khác nhau được biết tới trong lĩnh vực 
động lực học kết cấu dưới tên gọi họ các phương pháp Newmark (Newmark family 
of methods). Véc-tơ chuyển vị nút và vec-tơ vận tốc nút thời điểm mới (i+1) t trong 
họ các phương pháp Newmark có thể viết dưới dạng tổng quát như sau 
2
1 1
1 1
1 2 2
2
1
i i i i i
i i i i
tt
t
 
 
D D D D D
D D D D
  
   
 (3.43) 
trong đó  và  là các hằng số được lựa chọn để kiểm soát tính hội tụ của thuật toán 
số và độ chính xác của lời giải số. Tùy theo giá trị của  và  ta có các thuật toán 
tích phân trực tiếp khác nhau. 
Một đặc tính quan trọng trong việc lựa chọn phương pháp tích phân trực tiếp 
đó là tính ổn định của phương pháp trên quan điểm số. Người ta đã chứng minh 
được rằng, phương pháp sẽ là ổn định không điều kiện nếu các hằng số  và  thỏa 
mãn điều kiện [131] 
52 
 12
2
  (3.44) 
3.10.2. Phương pháp gia tốc trung bình 
Trên quan điểm toán học, véc-tơ chuyển vị và vận tốc tại các điểm nút trong 
phương pháp gia tốc trung bình có thể nhận được từ khai triển Taylor của véc-tơ 
chuyển vị nút quanh thời điểm i t và (i+1) t. Cụ thể 
2
1
2
1 1 1
2
2
i i i i
i i i i
tt
tt
D D D D
D D D D
 
 
 (3.45) 
Cộng và trừ các phương trình (3.45) cho nhau và bỏ qua các số hạng bé bậc cao ta 
nhận được 
1 1
1 1
( )
2
( )
2
i i i i
i i i i
t
t
D D D D
D D D D
 
   
 (3.46) 
từ phương trình (3.46) ta cũng có thể giải để tính các véc-tơ vận tốc và gia tốc nút 
tại thời điểm (i+1) t như sau 
1 1
1 12
2 ( )
4 4( )
i i i i
i i i i i
t
t t
D D D
D D D D D
D 
 


 (3.47) 
Phương trình chuyển động (3.40) viết tại thời điểm (i+1) t là 
 1
ex
1 1i i i MD KD F (3.48) 
Kết hợp phương trình (3.47) với phương trình (3.48) ta nhận được phương 
trình để xác định véc-tơ chuyển vị nút tại thời điểm mới (i+1) t như sau 
 ef ef1 1i i K D F (3.49) 
trong đó Kef, Fef tương ứng là các ma trận độ cứng và véc-tơ lực nút hữu hiệu 
(effective stiffness matrix and effective load vector) với biểu thức cụ thể như sau 
ex
ef
2
. ..
ef
1 1 2
4
4 4
i ii i i
t
t t 
K M K
F F M D D D
 (3.50) 
53 
Như vậy các phương trình (3.48), (3.49) và (3.50) cho phép ta hoàn toàn xác 
định được các véc-tơ chuyển vị, vận tốc và gia tốc tại thời điểm mới (i+1) t. 
3.11. Véc-tơ lực nút 
Véc-tơ lực nút Fex trong phương trình (3.37) nhận được bằng cách nối ghép 
véc-tơ lực nú