Bài giảng Phương pháp số - Vũ Khắc Bảy

1 1
1 2
2 3
2 4
e
v q
q
q
v q
q
   
          
   
      
Mỗi chuyển vị vi sẽ có 2 thành phần theo phương x' và y' trong tọa độ tổng thể , do vậy 
véc tơ chuyển vị nút phần tử trong tọa độ tổng thể sẽ có 6 thành phần : 
 
11
21
31
42
52
62
e
qu
qv
q
q
qu
qv
q
   
      
             
   
   
      
và {q}e = [T]e . {q'}e với 
 
0 0 0 0
0 0 1 0 0 0
0 0 0 0
0 0 0 0 0 1
y y
e
y y
m
T
m
 
 
   
 
  


 (III.41) 
các y và ym là các cosin chỉ phương của 0y với các trục 
0x' , 0y'. 
nếu gọi  là góc lập bởi trục phần tử dầm với 0x' thì y = - sin ; my = cos , 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
72 
do vậy ta vẫn gọi  = cos ; m = sin và sẽ dẫn đến: 
  
0 0 0 0
0 0 1 0 0 0
0 0 0 0
0 0 0 0 0 1
e
m
T
m
 
 
 
 
 
 


 (III.42) 
sử dụng (III.32) và [K']e = [T]eT [K]e . [T]e ta được : 
 Ma trận cứng của phần tử dầm chịu uốn trong khung phẳng: 
 
2 2
2 2
2 2
3 2
2
2
12 12 6 12 12 6
12 6 12 12 6
4 6 6 2
12 12 6
12 6
4
e
m m am m m am
a m a
a am a aEJ
K
a m m am
a
a
    
 
 
 
     
 
 
  
 
    


 
 (III.43) 
 Mô men uốn tính theo trong tọa độ tổng thể : 
Theo (III.37) và (III.38) ta có 
      1e ee
2
M
M S q
M
 
  
 
 với  
2 2
3 2 2
6 4 6 2
6 2 6 4e
a a a aEJ
S
a a a a a
   
  
  
Thay {q}e = [T]e . {q'}e dẫn đến : 
               e e e ee e e eM S q S T q S q     với      e e eS S T  
 
2 2
3 2 2
6 6 4 6 6 2
6 6 2 6 6 4e
am a a am a aEJ
S
a am a a am a a
     
   
  
 
 
 (III.44) 
3. Phần tử dầm chịu uốn và kéo , nén trong khung phẳng 
Chuyển vị tại mỗi nút của phần tử dầm sẽ theo 2 phương : trục dầm và vuông góc với 
trục dầm . Véc tơ chuyển vị nút của phần tử dầm sẽ là : 
đối xứng 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
73 
  
1
1
1
2
2
2
e
u
v
q
u
v
 
 
 
  
  
 
 
 
  
 và trong tọa độ tổng thể sẽ là :  
11
21
31
42
52
62
e
qu
qv
q
q
qu
qv
q
   
      
             
   
   
      
 (III.45) 
Do tính chất cộng tác dụng nên ta sử dụng ở đây các công thức (III.19) và (III.43) ta 
được: 
Ma trận cứng phần tử dầm chịu uốn và kéo nén dọc trục trong tọa độ tổng thể : 
 (III.46) 
 
2 2 2 2
2 2 2 2
2 2
2 2
2 2
2
12 12 6 12 12 6
12 6 12 12 6
4 6 6 2
12 12 6
12 6
4
e
F. m H F m mH amH F. m H F m mH amH
Fm H a H F m mH Fm H a H
a H amH a H a HE
K
a F. m H F m mH amH
Fm H a H
a H
        
 
     
 
      
  
 
  
     
     

  
 
trong đó :  = cos  ; m = sin ;  - góc lập bởi trục phần tử dầm và 0x' 
 a – độ dài thanh dầm F – diện tích thiết diện ; H = 2
J
a
Ví dụ 
Tính chuyển vị và các nội lực của các 
thanh trong khung phẳng chịu kéo nén dọc 
trục và chịu uốn. Các thanh đều có độ dài 
a, diện tích thiết diện F, mô đun đàn hồi E, 
chịu tác dụng lực ngoài và các điều kiện 
biên như hình vẽ. 
( lực p0 lập với các thanh góc 450) 
đối xứng 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
74 
Giải: 
1) Rời rạc hóa các kết cấu: đánh số phần tử và số nút => tính ma trận chỉ số [b] 
 Phần tử Nút đầu Nút cuối 
(1) 1 2 
(2) 2 3 
=> ma trận chỉ số [b] 
Nút đầu i Nút cuối j Chỉ số đ.ph 
Phần tử 1 2 3 4 5 6 
(1) 1 2 3 4 5 6 
(2) 4 5 6 7 8 9 
2)Thiết lập ma trận cứng tổng thể: 
 Chú ý rằng tại các nút 1 và 3 do bị ngàm chặt nên các chuyển vị và góc xoay bằng 
0 tức là ta có q'1 = q'2 = q'3 = q'7 = q'8 = q'9 = 0 nên trong các ma trận cứng phần tử 
các hàng hoặc các cột là các số 1 ; 2 ; 3 ; 7 ; 8 ; 9 thì được "loại ra" do vậy áp 
dụng (III.46) ta được : 
 Với phần tử (1) : 0 1
2
; m

     
2
4 5 6
12 0 6 4
0 0 5
66 0 4
*
e
H aH
E
K F
a
aH a H
 
      
 
 
 Với phần tử (2) : 0 1 0; m     
2
4 5 6
0 0 4
0 12 6 5
60 6 4
*
e
F
E
K H aH
a
aH a H
 
      
 
 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
75 
 Ma trận cứng tổng thể *K   ( đã tính đến các điều kiện biên) 
2
12 0 6
0 12 6
6 6
*
H F aH
E
K H F aH
a
aH aH 8a H
 
      
 
 
 ; H = 2
J
a
3) Véc tơ tải tổng thể : 
 Có  
1 1
1 1
1
2
1
1
2
30 0
4
5
60 0
y y
x x
R R
R R
p
Q qa
Q qa
    
   
   
           
   
       
      
 Véc tơ tải trên phần tử dầm (2) do tải 
trọng đều q và trọng tải nút P0 tác động. 
Theo (III.13) ta có  
2
2
2
2
12
2
12
q
aq
a q
p
aq
a q
  
 
   
  
 
 
 
 
  
mà    12 2
q qTp [T] p  
áp dụng (III.42) với : 0  nên 1 0  ; m , do đó ta được 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
76 
 
2 2
2
2
2
2
1 0 0 0 02
0 0 0 0
0 1 0 0 12 12
0 0 1 0
2 20 0 0 0
00 0 0 1
12
12
q
aq
aq
a q a q
p
aq aq
a q
a q
 
                             
      
    
    
     
    
 
 
 ; véc tơ tải trọng nút  
2
1
22
2
0
0
nut
( )
x
( )
y
Q
Q
p'
R
R
 
  
    
 
 
 
  
Do vậy véc tơ tải phần tử (2) trong tọa độ tổng thể là : 
     
2
2 1
2 22 1
22 2 2 2 2
2
2 2
2 2 2
3
2 22
0
0 12 1212
2 22
0
0
12 12
     
                                      
      
     
     
      
   
   
   
nut q
( )
( ) ( )x
x x
( )
y ( ) ( )
y y
aq aqaq Q
Q Q aq
Q a q a qa q
p p p
aq aqRaq R R
R
R R
a q a q a
4
5
6
7
8
912
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
q
Vậy có véc tơ tải tổng thể : 
 
1
1
2
2
2
2
1
2
30
5
42
52
612
72
8
912
 
 
 
 
 
 
 
 
 
 
  
 
 
 
 
 
 
 
 
 
 
( )
y
( )
x
( )
x
( )
y
R
R
aq
qa
p ' a q
aq
R
R
a q
 =>  
2
5
2
2
12
*
aq
p ' qa
a q
 
 
  
  
 
 
  
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
77 
và ta có hệ phương trình : 
 
4
5
2 2
6
5
12 0 6 2
0 12 6 2
6 6
12
* *
aq
H F aH q
E
K q H F aH q qa
a
qaH aH 8a H a q
 
    
                 
         
  
 với H = 2
J
a
Giải hệ phương trình trên ta tính được các giá trị 4 5 6q , q ,q 
Chú ý : Để giải hệ phương trình trên ta nên đưa các ẩn cần tìm về cùng một thứ nguyên: 
4 2
5
6
12 0 6 30
0 12 6 24
12
6 6 1
H F H q
a q
H F H q
E
H H 8H aq
     
          
          
4) Tính các nội lực: mô men uốn 
Chú ý: 
khi tính mô men uốn trên phần tử thanh (2) cần cộng thêm mô men M0(x) theo (III.39) 
III.3 Hệ khung không gian 
1. Ma trận cứng phần tử : 
Ta xét phần tử khung không gian là dầm thẳng có thiết diện không đổi, trên mặt cắt 
ngang tồn tại lực dọc, mô men uồn tròn cả hai mặt phẳng quán tính chính và mô men 
xoắn. Các bậc tự do chuyển vị đặc trưng cho trạng thái chuyển vị - biến dạng của phần tử 
dầm 2 nút được biểu diễn trên hình. Hệ tọa độ địa phương xyz với trục x là trục dầm, y 
và z là hai trục chính của mặt phẳng cắt ngang 
 Véc tơ chuyển vị phần tử hai nút là : 
{q}(e) = { q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12}T 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
78 
 Trong đó : 
q1 và q7 chuyển vị dọc trục dầm gây ra biến dạng dọc trục dầm 
q2 và q8 chuyển vị thẳng theo trục y 
q6 và q12 : góc xoay trong mặt phẳng xy 
gây ra biến dạng uốn trong mặt phẳng xy 
q3 và q9 chuyển vị thẳng theo trục z 
q5 và q11 : góc xoay trong mặt phẳng xz 
gây ra biến dạng uốn trong mặt phẳng xz 
q4 và q10 góc xoắn quanh trục x 
gây ra biến dạng xoắn của thanh 
Như vậy trong 12 bậc tự do chuyển vị sẽ gây ra 4 nhóm biến dạng độc lập nhau. Vì vậy 
ma trận cứng phần tử [ K ](e) có kích thước 12 × 12 được thành lập từ 4 ma trận con như 
sau : 
a) Biến dạng dọc trục ( q1 và q7) 
1 7
1
(e)
7
q q
q1 1E.F[K ]
q1 1L
 
   
 (III.47) 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
79 
b) Biến dạng xoắn theo trục x ( q4 và q10) 
Thấy rằng đa thức xấp xỉ của hàm góc xoắn 
θ (x) là hàm bậc nhất : 
 θ (x) = a1 + a2x ( 0 ≤ x ≤ L ) 
 hay : 
 1
2
a
θ(x) [1 x ] .
a
 
  
 
 = [P(x)] . {a} , theo kết quả của ví dụ 1 phần (II.2) ta có : 
biểu diễn xấp xỉ hàm góc xoắn θ (x) theo các góc xoắn tại các nút: 
 θ (x) = [N]e . { q}xoắn = 
x1
L
  
 
q4 + 
x
L
q10 (III.48) 
ở đây [N]e = 
x x1
L L
      
Nếu là thanh tròn thì biến dạng xoắn sẽ là xyz
dr
dx

  và ứng suất tiếp xy xyG.   
( theo Hooke và G là mô đun trượt , r là khoảng cách từ tâm đến điểm tính toán) 
áp dụng  ε xoắn =   θ xoắn =  [N]{q} xoắn = [B] {q} xoắn 
  σ xoắn = { yzτ }=  [T] q xoắn (III.49) 
trong đó : [T] = [D]. [B] là ma trận tính ứng suất 
 gọi [Kxoắn ]e - Ma trận cứng phần tử 
[Kxoắn ]e 
e
T
V
[B] [D][B]dv  (III.50) 
ở đây ta có { } = { yz } và { } = [ B ] . {q }(xoắn) 
với B = 
r r
L L
  
 
 ; [D ] = [ G ] 
Do đó ta được : 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
80 
 [Kxoắn ]e = 
L
2
0 F
1
1 1LG dx r dF. .
1 L L
L
         
  
  (III.51) 
=> [Kxoắn ]e = 
4 10
4x
10
q q
q1 1GJ
q1 1L
 
  
 (III.52) 
ở đây Jx là mô men quán tính độc cực của mặt cắt ngang:
4
x
RJ
2

 ; R – bán kính dầm 
Trong trường hợp nếu mặt cắt ngang dầm là hình chữ nhật có các cạnh a × b thì : 
 Jx = c ab3 với c là hằng số được lấy theo bảng sau ( phụ thuộc vào tỷ số a/b) 
a/b 1 1,5 2 3 5 10 
c 0,141 0,196 0,229 0,263 0,291 0,312 
c) Biến dạng uốn trong mặt phẳng xz ( do { q}xz = { q3 , q5 , q9 , q11}T ) 
áp dụng 
 
3 5 9 11
3
2 2
y 5
xz 3e
9
2 2
11
q q q q
12 6L 12 6L q
EJ q6L 4L 6L 2L
K
q12 6L 12 6LL
q6L 2L 6L 4L
 
 
     
 
  
 (III.53) 
với 2y
F
J z dF  - mô men quán tính của mặt cắt ngang đối với đường trung hòa. 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
81 
 Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jy = 
3ab
12
 Nếu thiết diện tròn bán kính R thì : 
4
y
RJ
4

 
d) Biến dạng uốn trong mặt phẳng xy ( do { q}xz = { q2 , q6 , q8 , q12}T ) 
áp dụng 
2 6 8 12
2
2 2
6z
xy 3e
8
2 2
12
q q q q
12 6L 12 6L q
q6L 4L 6L 2LEJK
q12 6L 12 6LL
q6L 2L 6L 4L
 
 
        
 
  
 (III.54) 
với 2z
F
J y dF  - mô men quán tính của mặt cắt ngang đối với đường trung hòa. 
 Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jz = 
3a b
12
 Nếu thiết diện tròn bán kính R thì : 
4
z
RJ
4

 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
82 
 Từ các kết quả trên ta xây dựng ma trận cứng phần tử [ K](e) là ma trận cỡ 12 × 12 
2. Chuyển ma trận cứng phần tử [K](e) về ma trận cứng [ K' ](e) trong hệ tọa độ tổng 
thể 
Để ghép nối các phần tử ta phải chuyển ma trận cứng phần tử về ma trận cứng [ K' ](e) 
trong hệ tọa độ tổng thể 
Ta có [ K' ](e) =      T(e) (e) (e)T K T 
trong đó [T](e) là ma trận chuyển hệ trục tọa độ , là ma trận vuông cấp 12, nó có dạng : 
  (e)
[n] [0] [0] [0]
[0] [n] [0] [0]
T
[0] [0] [n] [0]
[0] [0] [0] [n]
 
 
 
 
 
 
 (III.55) 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
83 
với : 
 [0] = 
0 0 0
0 0 0
0 0 0
 
 
 
  
 và [n] = 
x x x
y y y
y z z
m n
m n
m n
 
 
 
  



 (III.56) 
trong đó  x , mx , nx là cosin chỉ phương của trục x 
  y , my , ny là cosin chỉ phương của trục y 
  z , mz , nz là cosin chỉ phương của trục z trong hệ tổng thể x'y'z' 
khi đó : 
x x '
y [n]. y '
z z '
   
      
      
3. Các bước thiết lập chương trình tính cho hệ khung có Ne thanh với R nút: 
Bước 1. Lập bảng đánh số thanh- nút : có Ne thanh với R nút 
Nút T.T 
thanh Đầu : i Cuối : j 
1 1 2 
... ... ... 
k ki kj 
... ... ... 
Ne Nei Nej 
Ở bảng này đánh số và nhập bằng tay ( nếu số lượng ít) hoặc viết một đoạn chương trình 
để tự động đánh số thứ tự các thanh và các nút. Số liệu ghi vào mảng được khai báo : 
 Thanh_nut( 1 to Ne, 1 to 2) as integer 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
84 
Bước 2. Lập bảng nhập các thông số của các thanh dầm : 
 (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) 
TT 
thanh 
L E G Jx Jy Jz x mx nx y my ny z mz nz 
1 
... 
k 
... 
Ne 
Số liệu ở bảng này được ghi vào mảng : 
 Thanh(1 to Ne, 1 to 15) as double 
 Chú ý : Các trục chính x , y , z ở các thanh có gốc là nút đầu và phải lập thành 
tam giác thuận. 
Bước 3. Lập bảng ma trận chỉ số [b] : Thanh – Nút - Bậc tự do : 
Nút đầu : i Nút cuối : J T.T 
thanh q6i -5 q6i -4 q6i -3 q6i -2 q6i -1 q6i q6j-5 q6j-4 q6j-3 q6j-2 q6j-1 q6j 
1 
... 
k 
... 
Ne 
Số liệu sẽ tự tính và ghi vào mảng : 
 b(1 to Ne , 1 to 12) as integer 
Bước 4. Tính các giá trị của [T](e) được ghi vào mảng: 
 T(1 to 12 , 1 to 12) as double 
Bước 5. Tính [ K](e) được ghi vào mảng : 
 Ke(1 to 12, 1 to 12) as double 
Bước 6. Tính [K'](e) =  T(e)T . [K](e) . [T](e) được ghi vào mảng : 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
85 
 KeT(1 to 12, 1 to 12) as double 
Bước 7. Ghép [K'](e) vào ma trận cứng tổng thể [K'] được ghi vào mảng: 
 KT ( 1 to N, 1 to N) as double 
 Ở đây N = 6 . R tổng số bậc tự do 
Bước 8. Tính véc tơ tải phần tử => véc tơ tải tổng thể 
Bước 9. Áp các điều kiện biên 
Chú ý rằng sau khi có được ma trận cứng tổng thể [K'] ( mà định thức của nó = 0) , ta sẽ 
tính cho véc tơ tải tổng thể trên cơ sở tính các véc tơ tải phần tử. Nhờ các điều kiện biên 
ta sẽ bỏ đi được một số phương trình và một số biến ( số lượng của chúng bằng nhau), 
điều đó có nghĩa là ma trận cứng tổng thể [K'] sẽ bị bỏ đi một số hàng và một số cột và 
sẽ thành ma trận cứng tổng thể [K'* ] . Ma trận [K'* ] là ma trận vuông có cấp nhỏ hơn 
N và định thức của nó  0 , do vậy sẽ dẫn về một hệ Cramer để giải các biến là các bậc 
tự do ma ta cần tìm. 
Bước 10. Giải hệ phương trình tính các chuyển vị kq ( giá trị các bậc tự do theo hệ tọa 
độ tổng), từ đó tính {q}e = [T]e. {q’}e 
Bước 11. Tính các đại lượng : biến dạng , ứng suất . 
3. Một số thủ tục chương trình 
 Chương trình được viết với các khai báo sau : 
Dim Thanh_nut( 1 to Ne, 1 to 2) as integer ' đánh số các thanh và các nút 
Dim Thanh(1 to Ne, 1 to 15) as double ' các thông số về các thanh 
Dim b(1 to Ne , 1 to 12) as integer 'Ma trận chỉ số 
Dim T(1 to 12 , 1 to 12) as double ' Ma trận chuyển tọa độ 
Dim Ke(1 to 12, 1 to 12) as double ' Ma trận cứng phần tử 
Dim KeT(1 to 12, 1 to 12) as double ' Ma trận cứng phần tử trong hệ tọa độ tổng thể 
Dim KT ( 1 to N, 1 to N) as double ' Ma trận cứng tổng thể 
Dim i , j , k as integer 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
86 
Với các khai báo như trên, ta có được một số thủ tục sau: 
1. Thủ tục tạo lập ma trận chỉ số : [b] 
 Private sub TaoChiSo() 
 dim kk as integer 
 dim s , dau , cuoi as integer 
 For kk =1 to Ne 
 dau = Thanh_nut(kk,1) : cuoi = Thanh_nut(kk,2) 
 for s =1 to 6 
 b(kk,s) = 6*(dau -1) + s 
 b(kk,6+s) = 6*(cuoi -1) +s 
 next 
 Next 
 end sub 
2. Thủ tục tính Te 
 Private sub Te(kk as integer) 
 dim r , s as integer 
 for r =1 to 12 
 for s = 1 to 12 
 T(r,s) = 0 
 next 
 next 
 For r =1 to 3 
 for s = 1 to 3 
 T(r,s) = Thanh(kk, 7 + (r-1)*3 + s) 
 T(3+r,3+s) = T(r,s) 
 T(6+r,6+s) = T(r,s) 
 T(9+r,9+s) = T(r,s) 
 next 
 next 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
87 
 end sub 
3. Thủ tục tính ma trận cứng phần tử : 
 Private Sub MaTranCungPT(kk) 
 dim i, j as integer 
 dim E , G , L, F , Jx , Jy , Jz as double 
 for i=1 to 12 
 for j=1 to 12 
 Ke(i,j) = 0 
 next 
 Next 
L = Thanh(kk,1) : F = Thanh(kk,2) : E = Thanh(kk,3) : G = Thanh(kk,4) 
Jx = Thanh(kk,5) : Jy = Thanh(kk,6) : Jz = Thanh(kk,7) 
Ke(1,1)= E*F/L : Ke(1,7) = -E*F/L 
Ke(2,2)= 12*E*Jz/L^3 : Ke(2,6)= 6*E*Jz/L^2 : Ke(2,8)= -Ke(2,2) : Ke(2,12)= Ke(2,6) 
Ke(3,3)= 12*E*Jy/L^3 : Ke(3,5)= 6*E*Jy/L^2 : Ke(3,9)= -Ke(3,3) : Ke(3,11)= Ke(3,5) 
Ke(4,4)= G*Jx/L : Ke(4,10) = - Ke(4,4) 
Ke(5,5) = 4*E*Jy/L : Ke(5,9) = - 6*E * Jy/L^2 : Ke(5,11) = 2*E*Jy/L 
Ke(6,6)= 4*E*Jz/L : Ke(6,8) = - 6*E * Jz/L^2 : Ke(6,12) = 2*E*Jz/L 
Ke(7,7) = E*F/L 
Ke(8,8) = 12*E*Jz/L^3 : Ke(8,12) = - 6*E*Jz/L^2 
Ke(9,9) = 12*E*Jy/L^3 : Ke(9,11) = - 6*E*Jy/L^2 
Ke(10,10) = G*Jx/L : Ke(11,11) = 4*E*Jy/L : Ke(12, 12) = 4*E*Jz/L 
for i =2 to 12 
 for j=1 to i -1 
 Ke(i,j) = Ke(j,i) 
 next 
Next 
End sub 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
88 
4. Thủ tục tính [K'](e) : Ma trận cứng phần tử trong hệ tọa độ tổng thể , 
 [K'](e) = [T]T.[K](e).[T] 
Private Sub MaTranCungPT_Tong 
 Dim MaTranTG(1 to 12, 1 to 12) as double 
 dim i, j , s as integer 
 dim tg as double 
for i =1 to 12 
 for j =1 to 12 
 MaTranTG(i,j) = 0 
 KeT(i,j) = 0 
 next 
Next 
For i =1 to 12 
 for j = 1 to 12 
 tg = 0 
 for s = 1 to 12 
 tg = tg + Ke(i,s) * T(s,j) 
 next 
 MaTranTG(i,j) = tg 
 next 
Next 
For i =1 to 12 
 for j = 1 to 12 
 tg = 0 
 for s = 1 to 12 
 tg = tg + T(s,i) * MaTranTG(s,j) 
 next 
 KeT(i,j) = tg 
 next 
Next 
End sub 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
89 
5. Thủ tục ghép nối ma trận cứng phần tử vào ma trận cứng tổng thể. 
 Private sub GhepNoi(k as integer) 
 dim i , j as integer 
 for i = 1 to 12 
 for j=1 to 12 
 KT(b(k,i) , b(k,j)) = KT(b(k,i) , b(k,j)) + Ke(i,j) 
 next 
 Next 
 end sub 
6. Thủ tục tính toán chính , Thủ tục này được dùng sau khi đã nhập các số liệu vào các 
mảng Thanh, Thanh_nut 
 Private Sub Tinh 
 Dim K , Kj as integer 
call TaoChiSo 
For k = 1 to N 
 for kj =1 to N 
 KT (k, kj) = 0 
 next 
Next 
For k = 1 to N 
 call Te(k) 
 call MaTranCungPT(k) 
 call MaTranCungPT_Tong 
 call GhepNoi(k) 
Next 
end sub 
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 
Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 
90 
 Tài liệu tham khảo 
1. Rao S.S. The Finite Element Method in Enginieering, Second Edition, Pegamon 
Press, 1989. 
2. Zienkiewicz O.C and Taylor R.L. The Finite Element Method, volum 1, 2. 
4th Edition McGraw – Hill Book Co. , 1991 
3. Dũng N.L. Giáo trình phương pháp phần tử hữu hạn trong cơ học,Trường Đại học 
Bách khoa TP Hồ Chí Minh, 1993 
2. Chu Quốc Thắng, Phương pháp phần tử hữu hạn, NXB Khoa học và kỹ thuật, 1997