Bài giảng Phương pháp sai phân hữu hạn & phần tử hữu hạn trong truyền nhiệt - Trịnh Văn Quang

1 
PHƯƠNG PHÁP SAI PHÂN HỮU HẠN 
& PHẦN TỬ HỮU HẠN 
TRONG TRUYỀN NHIỆT 
Bài giảng môn Truyền nhiệt 
cho các lớp cao học Cơ khí 
 PGS. TS Trịnh Văn Quang 
Bộ môn Kỹ Thuật nhiệt – Khoa Cơ khí 
ĐẠI HỌC GIAO THÔNG VẬN TẢI HÀ NÔI 
Hà nội -2009 
2 
Mục lục 
Lời nói đầu 3 
PHẦN 1. PHƯƠNG PHÁP SAI PHÂN HỮU HẠN 
2.1 . Bài toán ổn định hai chiều 4 
 1. Phương trình sai phân hữu hạn 4 
 2. Xây dựng hệ phương trình bậc nhất 4 
2.2 . Bài toán dẫn nhiệt không ổn định một chiều 5 
 1. Phương pháp Ma trận nghịch 5 
 2. Phương pháp tính lặp 6 
2.3. Bài toán dẫn nhiệt không ổn định hai chiều 9 
2.4. Giải hệ phương trình tuyến tính của nhiệt độ 13 
 1. Phương pháp định thức 13 
 2. Phương pháp Gauss 13 
 3. Phương pháp Gauss - Jordan 15 
 4. Phương pháp Gauss - Seidel 17 
PHẦN 2. PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 
Giới thiệu khái quát 20 
2.5. Nội dung cơ bản, trình tự giải bài toán nhiệt bằng phương pháp PTHH 20 
2.6. Các phần tử và hàm nội suy 23 
2.6.1. Phần tử một chiều bậc nhất 23 
2.6.2. Phần tử một chiều bậc hai 25 
2.6.3. Phần tử hai chiều tam giác bậc nhất 29 
2.6.4. Phần tử chữ nhật bậc nhất 36 
2.6.5. Các phần tử đẳng tham số 38 
2.7. Thiết lập phương trình đặc trưng phần tử đối với phương trình vi phân dẫn nhiệt 46 
2.7.1. Phương pháp biến phân 47 
2.7.2. Phương pháp Galerkin 53 
2.8. Giải bài toán dẫn nhiệt một chiều bằng phương pháp PTHH 54 
2.9. Dẫn nhiệt qua vách phẳng có nguồn nhiệt bên trong 59 
2.10. Dẫn nhiệt qua vách trụ 67 
2.11. Dẫn nhiệt qua thanh trụ có nguồn trong 71 
2.12. Dẫn nhiệt qua cánh tiết diện thay đổi 75 
2.13. Dân nhiệt ổn định hai chiều dùng phần tử tam giác 80 
2.14. Dẫn nhiệt hai chiều qua phần tử chữ nhật 99 
3 
Lời nói đầu 
Do yêu cầu giải quyết các bài toán thực tế, nhiều năm qua đã có nhiều phương pháp số phát triển. 
Phương pháp phổ biến nhất được sử dụng trong kỹ thuật tính nhiệt là các phương pháp sai phân hữu hạn, 
thể tích hữu hạn và phần tử hữu hạnngoài ra còn có phương pháp phần tử biên giới. ở đây nêu nội dung 
cơ bản của ba phương pháp đầu. 
- Phương pháp Sai phân hữu hạn (SPHH) là phương pháp số tương đối đơn giản và ổn định. Nội dung 
của phương pháp này là biến đổi một cách gần đúng các đạo hàm riêng của phương trình vi phân chủ đạo 
thành thương của các số gia tương ứng. Bằng cách dùng các họ đường song song với các trục toạ độ để 
tạo thành một mạng lưới chia miền nghiệm trong vật thể thành một số hữu hạn các điểm nút, rồi xác định 
nhiệt độ của phẫn tử tại các nút đó thay cho việc tính nhiệt độ trên toàn miền. Như vậy phương pháp 
SPHH đã xấp xỉ các phương trình vi phân đạo hàm riêng thành các phương trình đại số. Kết quả thiết lập 
được hệ phương trình đại số gồm n phương trình tương ứng với giá trị nhiệt độ của n nút cần tìm. 
 Mức độ chính xác của nghiệm trong phương pháp SPHH có thể được cải thiện nhờ việc tăng số điểm nút. 
Phương pháp SPHH rất hữu hiệu trong việc giải nhiều bài toán truyền nhiệt phức tạp mà phương pháp 
giải tích gặp khó khăn. Bởi vậy trong các giáo trình truyền nhiệt hiện đại, phương pháp SPHH được trình 
bày khá kỹ cho chương trình đại học (Holman ..). Tuy nhiên khi gặp phải vật thể có hình dạng bất quy tắc 
hoặc điều kiện biên giới bất thường, phương pháp SPHH cũng có thể khó sử dụng. 
- Phương pháp thể tích hữu hạn (TTHH) có tinh tế hơn phương pháp SPHH và trở nên phổ biến trong kỹ 
thuật tính nhiệt và động học dòng chảy (Patankar 1980). Trong tính nhiệt, phương pháp TTHH dựa trên 
cơ sở cân bằng năng lượng của phân tố thể tích. Kỹ thuật thể tích hữu hạn tập trung vào điểm giữa phân tố 
thể tích rất tương tự với phương pháp SPHH (Malan et al 2002). 
4 
PHẦN 1. PHƯƠNG PHÁP SAI PHÂN HỮU HẠN 
2.1 . Bài toán ổn định hai chiều 
1. Phương trình sai phân hữu hạn 
Phương trình vi phân dẫn nhiệt ổn định hai chiều có dạng : 
 0
2
2
2
2




y
T
x
T
 (2.1) 
Xây dựng phương trình sai phân hữu hạn (SPHH) như sau : 
Chia vật thể bởi một mạng các đường vuông góc có bước mạng x, y, 
ứng với hai chiều x,y. Khi đó tại điểm nút i,j các đạo hàm bậc nhất và 
bậc hai của nhiệt độ viết dạng sai phân như sau (hình 2.1) : 
x
TT
x
T
x
T jiji

 ,1,
y
TT
y
T
y
T jiji

 1,,
 Hình 2.1. Mạng các điểm nút 
2
,1,,1
22
2
)(
)()(
)(
)(
x
TTTT
x
T
x
T jijijiji

 (2.2) 
2
1,,,1,
22
2
)(
)()(
)(
)(
y
TTTT
y
T
y
T jijijiji

 (2.3) 
Thay (2.2) và (2.3) vào phương trình vi phân (2.1) sẽ được : 
 0
)(
)()(
)(
)()(
2
1,,,1,
2
,1,,1 
y
TTTT
x
TTTT jijijijijijijiji (2.4) 
(2.4) là phương trình SPHH dẫn nhiệt viết cho điểm nút (i,j) 
2. Xây dựng hệ phương trình bậc nhất 
Để giải (2.4) , có thể chọn x = y. Khi đó sẽ được : 
)(
4
1
1,1,,1,1, jijijijiji TTTTT (2.5) 
Vậy nhiệt độ tại điểm nút bằng trung bình cộng của bốn điểm nút xung quanh . 
5 
Từ (2.5) viết lần lượt cho các điểm, rồi chuyển các nhiệt độ đã biết sang vế phải, các nhiệt độ chưa biết 
sang vế trái, sắp xếp lại sẽ được n phương trình cho n điểm nút chưa biết nhiệt độ bên trong vật, tạo thành 
hệ phương trình bậc nhất : 
nnnnnn
n
nn
CTaTaTa
CTaTaTa
CTaTaTa
...
............
...
...
2211
221222121
11212111
 (2.6) 
Từ đó có thể giải ra các nhiệt độ cần tìm bằng các phương pháp: Gauss, Gauss Seidel, Gauss Jordan, Ma 
trận nghịch đảo ... 
2.2 . Bài toán dẫn nhiệt không ổn định một chiều 
1. Phương pháp Ma trận nghịch đảo 
Phương trình vi phân dẫn nhiệt không ổn định 1 chiều : 
2
2
x
T
a
T





 (2.7) 
a. Các điểm bên trong vật 
Gọi p là thời điểm trước, (p+1) là thời điểm sau. Phương trình (2.7) được sai phân hoá như sau : 
 

 pi
p
i TTTT
1
 (2.8) 
Vế phải của (2.8) viết cho thời điểm sau (p+1) 
2
1
1
111
1
22
2
)(
)()(
)(
)(
x
TTTT
x
T
x
T pi
p
i
p
i
p
i

 
 (2.9) 
thay (2.8) và (2.9) vào (2.7): 
2
1
1
111
1
1
)(
)()(
x
TTTT
a
TT pi
p
i
p
i
p
i
p
i
p
i

 (2.10) 
(2.10) là phương trình SPHH dẫn nhiệt không ổn định 1 chiều, để giải (2.10) cần biến đổi: 
 pi
p
i
p
i
p
i
p
i TTTTT
x
a
11
1
11
12
)2(
)(
. 
 (2.11) 
Đặt 
2)(
.
x
a
Fo

 sẽ được 
 )2.( 11
11
1
1 
 pi
p
i
p
i
p
i
p
i TTTFoTT (2.12) 
vậy : 
 pi
p
i
p
i
p
i TFo.T Fo)T(-FoT 
1
1
11
1 21 (2.13) 
Phương trình (2.13) biểu thị các nhiệt độ tại thời điểm sau theo nhiệt độ tại thời điểm trước. 
b. Các điểm trên biên 
6 
Các điểm trên biên có i = 1. Phân tố bề mặt vật có bề dày x/2, diện tích y. z = 1m 1m, nhận nhiệt từ 
môi trường và nhiệt từ phân tố liền kề phía trong (i = 2) 
- Dòng toả nhiệt từ môi trường bên ngoài tới sau thời gian  : 
  111 ppKh TThq (2.14) 
- Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian : 
  
 )( 11
1
2
pp
k TT
x
k
q (2.15) 
Độ tăng nội năng dU phân tố sau thời gian  : 
 )(
2
)(. 1
1
11
1
1
pppp TT
x
cTTVcdU 
 (2.16) 
Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : 
 )(
2
)( 1
1
1
1
1
1
2
1
1
1 pppppp
K TT
x
cTT
x
k
TTh 
  (2.17) 
 ppppppK TTTT
xc
k
TT
xk
xh
c
k
1
1
1
1
1
1
22
1
1
1
2
)(
)(
2
)(
2 
 

 (2.18) 
Đặt Fo =
2)(
.
x
a
 
, 
k
xh
Bi
.
, 
 c
k
a ; Fo là tiêu chuẩn Phuriê, Bi là tiêu chuẩn Biô, a là hệ số khuyếch 
tán nhiệt độ sẽ được : 
 ppppppK TT)TFo(TTTBi.Fo 1111112111 22 
Chuyển nhiệt độ và các đại lượng đã biết sang vế phải 
 ppK
pp TTBi.FoTFoTFoBi.Fo 1
11
2
1
1 .2.2.122 
 (2.19) 
(2.13) và (2.19) là các phương trình dạng hàm ẩn đối với nhiệt độ cần tìm các điểm ở thời điểm sau theo 
nhiệt độ thời điểm trước và nhiệt độ môi trường. Từ đó có thể thành lập hệ phương trình tuyến tính các 
nhiệt độ cần tìm sau : 
nnnnnn
n
nn
CTaTaTa
CTaTaTa
CTaTaTa
...
............
...
...
2211
221222121
11212111
 (2.20) 
trong đó: 
 aij là các hệ số của nhiệt độ phải tìm, 
 Ti là nhiệt độ cần tìm ở thời điểm (p+1), viết gọn của 
1 P
iT 
 Ci là các hệ số chính là nhiệt độ đã biết ở thời điểm trước 
Hệ trên viết dạng ma trận như sau : 
    iiij CTa (2.21) 
trong đó: 
  ija là ma trận vuông gồm các hệ số của nhiệt độ phải tìm, 
 iT là ma trận cột gồm nhiệt độ cần tìm ở thời điểm (p+1) 
 iC là ma trận cột gồm các hệ số chính là nhiệt độ đã biết ở thời điểm trước 
Từ đó giải ra các nhiệt độ cần tìm tại thời điểm (p+1): 
7 
    iiji CaT 1 (2.22) 
   1 ija là ma trận nghịch đảo của [aii], 
Sau khi giải ra các nhiệt độ tại thời điểm nào đó, thì các nhiệt độ đã biết này trở thành hệ số [Ci] trong 
phương trình (2.22) để tính các nhiệt độ ở thời điểm tiếp theo 
2. Phương pháp tính lặp 
a. Các điểm bên trong vật 
Gọi p là thời điểm trước, (p+1) là thời điểm sau. Phương trình (2.7) được sai phân hoá như sau : 
 

 pi
p
i TTTT
1
 (2.23) 
Vế phải của (2.7) viết cho thời điểm trước p : 
2
11
22
2
)(
)()(
)(
)(
x
TTTT
x
T
x
T pi
p
i
p
i
p
i

 (2.24) 
thay (2.23)và (2.24)vào (2.7) : 
2
11
1
)(
)()(
x
TTTT
a
TT pi
p
i
p
i
p
i
p
i
p
i

 (2.25) 
Để giải (2.25) cần biến đổi như sau : 
 )2(
)(
.
112
1 p
i
p
i
p
i
p
i
p
i TTT
x
a
TT 

 ( 2.26) 
Đặt 
2)(
.
x
a
Fo

 sẽ được : 
 )2.( 11
1 p
i
p
i
p
i
p
i
p
i TTTFoTT 
Vậy : 
 pi
p
i
p
i
p
i TFoTFoTFoT 11
1 ..)21(. 
 (2.27) 
Phương trình (2.27) cho biết mỗi nhiệt độ tại vị trí i ở thời điểm sau (p+1) được tính theo các nhiệt độ ở 
thời điểm trước. Phương trình có dạng hàm tường, bởi vậy không thể lập ma trận được mà phải tính dần. 
Có thể áp dụng phương pháp tính lặp. 
8 
Để các nghiệm hội tụ cần điều kiện : 
 (1- 2Fo) 0 (2.28) 
 tức là : 
 Fo 
2
1
hay phải chọn bước thời gian đủ nhỏ : 
a
x
2
)( 2 
  (2.29) 
b. Các điểm trên biên 
Phân tố bề mặt vật có bề dày x/2, diện tích y. z = 1 1 nhận nhiệt từ môi trường và nhiệt từ phân tố 
phía trong 
- Dòng toả nhiệt từ môi trường bên ngoài tới sau thời gian  : 
  ppKh TThq 1 (2.30) 
- Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian : 
  
 )( 12
pp
k TT
x
k
q (2.31) 
- Độ tăng nội năng dU phân tố dày x/2 sau thời gian  : 
 )(
2
)(. 1
1
11
1
1
pppp TT
x
cTTVcdU 
 (2.32) 
- Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : 
 )(
2
)( 1
1
1121
pppppp
K TT
x
cTT
x
k
TTh 
  (2.33) 
Hay 
 ppppppK TTTT
xc
k
TT
xk
xh
c
k
1
1
112212
)(
)(
2
)(
2 
 

 (2.34) 
Với Fo =
2)(
.
x
a
 
, Bi = 
k
xh .
, a = 
 c
k
 sẽ được : 
 ppppppK TT)TFo(TTTBi.Fo 111121 22 
Chuyển nhiệt độ tại thời điểm sau cần tính sang vế trái, chuyển các đại lượng đã biết và nhiệt độ thời 
điểm trước sang vế phải 
 pppK
p FoTTBi.FoFoTBi.FoT 21
1
1 2)221(.2 
 (2.35) 
9 
(2.35) là phương trình dạng hàm tường cho biết nhiệt độ tại biên thời điểm sau 11
 pt theo nhiệt độ các 
điểm thời điểm trước. 
Điều kiện để xác định 11
 pT , tức nghiệm hội tụ cần phải thoả mãn : 
 (1- 2Fo -2Bi.Fo) 0 (2.36) 
2.3. Bài toán dẫn nhiệt không ổn định hai chiều 
Bài toán dẫn nhiệt không ổn định hai chiều, với điều kiện biên hỗn hợp loại 2 và loại 3 được mô tả bởi 
- Phương trình vi phân dẫn nhiệt ổn định hai chiều: 




  


2
2
2
2
2.
y
T
x
T
aTa
T

 (2.37) 
- Điều kiên biên loại 2 : với một biên giả sử là chữ nhật có x = 0  a; y = 0  b 
 qx = 0 = q1() ; qx = a = q2() (2.38) 
 qy = 0 = q3() ; qy = b = q4() 
- Điều kiện biên loại 3 : 
 T
k
h
x
T
x


1
0
; T
k
h
x
T
ax


2 
 T
k
h
x
T
y


3
0
; T
k
h
x
T
by


4 (2.39) 
Đối với các hình phức tạp không thể giải bằng phương pháp giải tích, nên phải dùng phương pháp số . Một 
trong các phương pháp số là PP SPHH được xây dựng như sau : 
Chia vật thể bởi một mạng các đường vuông góc 
có bước mạng x , y, ứng với hai chiều x,y. 
Khi đó tại điểm nút i,j các đạo hàm bậc nhất và 
bậc hai của nhiệt độ viết dạng sai phân như sau ( 
hình 2.2) : 
a. Các điểm bên trong vật 
Hình 2.2. Mạng các điểm nút 
Tại nút i, j , ở mỗi thời điểm các số hạng có thể viết 
10 
2
,1,,1
2
,1,,1
22
2
)(
.2
)(
)()(
)(
)(
x
TTT
x
TTTT
x
T
x
T jijijijijijiji

 (2.40) 
2
1,,1,
2
1,,,1,
22
2
)(
.2
)(
)()(
)(
)(
y
TTT
y
TTTT
y
T
y
T jijijijijijiji

 (2.41) 
Riêng đạo hàm theo thời gian luôn có 
 


 p
ji
p
ji TTTT ,
1
,
 (2.42) 
Viết (2.40), (2.41) ở thời điểm p rồi cùng với (2.42) thay vào phương trình vi phân (2.37) sẽ được : 
2
1,,1,
2
,1,,1,
1
,
)(
.2
)(
.2
. y
TTT
x
TTT
c
kTT
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
 
 (2.43) 
 Viết (2.40), (2.41) ở thời điểm (p+1) rồi cùng với (2.42) thay vào phương trình vi phân (2.37) sẽ được : 
2
1
1,
1
,
1
1,
2
1
,1
1
,
1
,1,
1
,
)(
.2
)(
.2
. y
TTT
x
TTT
c
kTT
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
 
 (2.44) 
(2.43) và (2.44) sẽ dẫn tới các hệ phương trình nhiệt độ tại các điểm nút bên trong vật, giải theo phương 
pháp khác nhau. 
- Từ (2.43) sẽ có: 
 pji
p
ji
p
ji
p
ji
p
ji
p
ji
p
jip
ji T
c
k
y
TTT
x
TTT
T ,2
1,,1,
2
,1,,11
, .
.)(
.2
)(
.2
 
 (2.45) 
 (2.45) là dạng hàm tường vì vế trái chưá một nhiệt độ tại điểm i,j ở thời điểm (p+1), phải giải bằng phương 
pháp tính thế dần. 
- Từ (2.44) sẽ có: 
 pji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji
p
ji tt
c
k
y
ttt
x
ttt
,
1
,2
1
1,
1
,
1
1,
2
1
,1
1
,
1
,1 .
.
.
)(
.2
)(
.2
 
 (2.46) 
(2.46) là dạng hàm ẩn vì chưá nhiệt độ các điểm ở thời điểm (p+1). (2.46) tạo thành hệ n phương trình bậc 
nhất, giải bằng phương pháp ma trận nghịch đảo, có thể chọn bước thời gian  tuỳ ý. 
Từ (2.45) và (2.46) có thể tìm được nhiệt độ tại các điểm bên trong vật. 
b. Các điểm trên biên 
11 
Các điểm trên biên phải áp dụng phương pháp cân bằng năng lượng trên phân tố thể tích . 
Tại bề mặt điều kiên loại 2 được quy về điều kiện loại 3 tại thời điểm p như sau : 
- Điều kiên loại 2 : 
 Dòng bức xạ là PI. )(qR  , với  là hệ số hấp thụ của vật, I
P là năng suất bức xạ chiếu tới 
- Điều kiên loại 3 : 
 Dòng đối lưu từ không khí là )( )(qK
P
m
P
K TTh  
- Dòng nhiệt tổng : 
 )(
.
. )( )(q Pm
P
K
P
m
P
P
K
PP
m
P
K TThT
h
I
ThITTh 
 

 (2.47) 
trong đó : 
P
m
P
K TT , là nhiệt độ không khí và nhiệt độ bề mặt của kết cấu 
 h ,  là hệ số toả nhiệt và hệ số hấp thụ của bề mặt 
h
I P.
 là nhiệt độ quy đổi của bức xạ 
h
I
TT
P
P
K
P
K
.
  là nhiệt độ tương đương của không khí có kể đến bức xạ 
Theo nguyên lý bảo toàn năng lượng thì tại phần tử thuộc nút (i,j) tổng các dòng nhiệt nhận dẫn đến phần 
tử từ xung quanh sau thời gian  bằng độ tăng nội năng của phần tử . Bởi vậy phương trình cân bằng 
năng lượng viết cho các phần tử (được giới hạn bởi đường nét đứt trong hình) như sau : 
 Hình 2.3 a Hình 2.3 b Hình 2.3 c Hình 2.3 d. 
+ Các phần tử bên trong mặt cắt , hình 2.3 a : Phần tử (i,j) rộng x , cao x, dài 1m : 
 x
y
k
TTx
y
k
TTy
x
k
TTy
x
k
TT pji
p
ji
p
ji
p
ji
p
ji
p
j
p
ji
p
ji
1
,
1
1,
1
,
1
1,
1
,
1
,11
1
,
1
,1 
 pjipji TTyx ,1,..c (2.4 ... và 22: 
 (-T3 + 2T6 – T5 ) = 0 (11) 
 ( -T5 + T6 + 0T8 ) = 0 (20) -T3 – 2T5 + 4T6 – T9 = 0 (30) 
 ( T6 – T9 + 0T8 ) = 0 (22) 
Nút 7: có các phương trình 15 và 18: 
 ( -T4 + 0T5 + T7) = 0 (15) 
 ( 0T5 – T8 + T7 ) = 0 (18) -T4 + 2T7 – T8 = 0 (31) 
Nút 8: có các phương trình 17,21 và 24: 
 (-T5 + 2T8 – T7 ) = 0 (17) 
 ( -T5 + 0T6 + T8 ) = 0 (21) - 2T5 - T7 + 4T8 – T9 = 0 (32) 
( 0T6 – T9 + T8 ) = 0 (24) 
Nút 9: chỉ có phương trình 23 : -T6 + 2T9 – T8 = 0 (33) 
. Chuyển 9 phương trình (25)(33) trên sang dạng ma trận 






0
0
0
0
0
0
0
0
0
210100000
141020000
012001000
100420100
020282020
001024001
000100210
000020141
000001012
9
8
7
6
5
4
3
2
1
T
T
T
T
T
T
T
T
T
90 
3.2. Phương pháp lắp ghép thực tế 
Thực tế khi bài toán có số nút lớn, việc chuyển các ma trận của từng phần tử sang dạng đại số rất mất 
công và thời gian, bởi vậy các bước trên được tóm lược cho gọn và đơn giản hơn như sau 
a. Đánh số các phương trình tại các nút có nhiệt độ tương ứng 
 - Viết phương trình ma trận đặc trưng của 9 phần tử và đánh số ở bên phải từ 1 đến 24 như sau 
Phần tử 1: Phần tử 2: 
 6
5
4
0
0
0
110
121
011
3
2
1
0
0
0
101
011
112
4
5
2
4
2
1












T
T
T
T
T
T
 Phần tử 3: Phần tử 4: 
 12
11
10
0
0
0
110
121
011
9
8
7
0
0
0
101
011
112
5
6
3
5
3
2












T
T
T
T
T
T
Phần tử 5: Phần tử 6: 
 18
17
16
0
0
0
110
121
011
15
14
13
0
0
0
101
011
112
7
8
5
7
5
4












T
T
T
T
T
T
 Phần tử 7: Phần tử 8: 
 24
23
22
0
0
0
110
121
011
21
20
19
0
0
0
101
011
112
8
9
6
8
6
5












T
T
T
T
T
T
24 số trên cũng biểu thị 24 phương trình đại số triển khai từ 9 phương trình ma trận của 9 phần tử tương 
ứng. Mỗi chữ số biểu thị một phương trình ở một nút có nhiệt độ tương ứng. Thí dụ: 
Phương trình (5) biểu thị phương trình nhiệt độ tại nút 5, đó là : -T2 + 2T5 – T4 = 0 
Phương trình (21) biểu thị phương trình nhiệt độ tại nút 8, đó là : -T5 + 0T6 + T8 = 0 
b. Lập bảng nguyên tắc lắp ghép 
Bậc tự do, số nút và số của phương trình được lập bảng để chỉ dẫn cho việc cộng các phương trình trong 
cùng nút như sau 
 Bảng 2.12 
Phần tử 1 2 3 4 5 6 7 8 
Thứ tự nút 
trong phần tử 
1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 
Số nút toàn cục 1 2 4 2 5 4 2 3 5 3 6 5 4 5 7 5 8 7 5 6 8 6 9 8 
Phương trình số 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
c. Bảng lắp ghép: 
91 
Để thể hiện mỗi nút có phương trình nào, hệ số của nhiệt độ trong mỗi phương trình, cần sắp xếp lại Số 
của phương trình theo Số nút và Hệ số của nhiệt độ có mặt trong các phương trình đó để tạo thành Bảng 
lắp ghép như sau 
 Bảng 2.13 
Nút 
số 
Số của 
phương 
trình tại 
mỗi nút 
Hệ số của nhiệt độ có mặt trong phương trình 
T1 T2 T3 T4 T5 T6 T7 T8 T9 
1 1 2 -1 -1 
2 2 -1 +1 
4 +1 -1 
7 +2 -1 -1 
3 8 -1 +1 
10 +1 -1 
4 3 -1 +1 
6 +1 -1 
13 +2 -1 -1 
5 5 -1 -1 +2 
9 -1 +1 
12 +1 -1 
14 -1 +1 
16 +1 -1 
19 +2 -1 -1 
6 11 -1 -1 +2 
20 -1 +1 
22 +1 -1 
7 15 -1 +1 
18 +1 -1 
8 17 -1 -1 +2 
21 -1 +1 
24 +1 -1 
9 23 -1 -1 +2 
d. Lập ma trận độ cứng tổng thể 
Trong mỗi cột nhiệt độ, cộng các hệ số của nhiệt độ tại mỗi nút lại sẽ được ma trận độ cứng tổng thể. 
Trong đó các cột của ma trận ứng với các nhiệt độ có các hệ số vừa được cộng ở trên, còn các hàng ứng 
với thứ tự các nút của các phần tử. Nghĩa là bảng lắp ghép trên cho ngay hình ảnh của ma trận hệ số nhiệt 
độ. 
92 






0
0
0
0
0
0
0
0
0
210100000
141020000
012001000
100420100
020282020
001024001
000100210
000020141
000001012
9
8
7
6
5
4
3
2
1
T
T
T
T
T
T
T
T
T
 (2.308) 
4. Giải phương trình 
Trong bài toán trên, các nhiệt độ tại biên đã biết chỉ có nhiệt độ T5 là chưa biết, nên chỉ cần từ phương 
trình của nút 5 giải ra: 
 - 2T2 – 2T4 + 8T5 – 2T6 – 2T8 = 0 (2.309) 
Thay các giá trị T2 = T4 = T6 = 100 
0C, T8 = 500 
0C vào sẽ được T5 = 1600/8 = 200 
0C 
Có thể so sánh kết quả với nghiệm chính xác bằng phương pháp giải tích của Holman: 
1
1
1
12
sinh
sinh
sin
1)1(2
)(),( T
H
b
n
y
b
n
x
b
n
n
TTyxT
n
n
 
 (2.310) 
trong đó T1 nhiệt độ cạnh trên , T2 nhiệt độ 3 cạnh còn lại của chữ nhật ; b, H bề rộng và cao chữ nhật. 
Với x = 0,5; y = 0,5 giải ra 
 T5(0,5;0,5) = 200,11 
0C 
 Nếu dùng phương pháp sai phân hữu hạn dễ dàng suy ra 
 CTTTTT 086425 200500100100100
4
1
4
1
 (2.311) 
Thấy rằng phương pháp phần tử hữu hạn cho kết quả đồng nhất với phương pháp giải tích và sai phân hữu 
hạn. 
Có thể xác định được nhiệt độ tại các điểm khác trong hình phẳng, tùy theo yêu cầu mức chính xác mà 
phải dùng mạng lưới tam giác mịn hơn. Nghiệm của mạng lưới tam giác có cấu trúc đều luôn chính xác 
hơn mạng lưới chia không đều. 
Thí dụ 2.13. Xác định nhiệt độ tại điểm 4(40;40) trong tam giác như trên hình 2.27. Biết các nút 1,2 và 3 
có nhiệt độ tương ứng là 1000C, 2000C và 1000C. Tọa độ các điểm đó tương ứng là (50;0), (50;50), và 
(0;50) 
93 
 Hình 2.27. 
Nhiệt độ tại các vị trí bên trong phần tử tam giác được xác định theo nhiệt độ ba nút: 
 T = N1T1 + N2T2 + N3T3 
Trong đó cần xác định các hàm hình dạng N1, N2 và N3 xác định theo 
 ycxba
A
N iiii 
2
1
 ; i = 1,2,3 (2.312) 
Toạ độ các nút: 
Nút 1 2 3 
Tọa độ 
(m) 
x 50 50 0 
y 0 50 50 
Xác định các hệ số ai, bi , ci : 
50500;05050;250050.050.50 23132123321 xxcyybyxyxa
50050;50050;250050.500 31213231132 xxcyybyxyxa
05050;50500;2500050.50 32321312213 xxcyybyxyxa 
Tính D = 2A: 
 2500
5001
50501
0501
2 
 A 
Tính các hàm nội suy N1, N2 và N3 
  yyxycxba
A
N 50
50
1
50.02500
2500
1
2
1
1111
  yxyxycxba
A
N 50
50
1
50502500
2500
1
2
1
2222
  )50(
50
1
50502500
2500
1
2
1
3333 xyxycxba
A
N 
   
3(0;50) 2(50;50) 
4(40,40) 
1(50;0) 
94 
Tại điểm 4 có x = 40; y = 40 
5
1
;
5
3
;
5
1
321 NNN (2.313) 
Tính nhiệt độ tại điểm 4: 
C03322114 1602040.320100
5
1
200
5
3
100
5
1
TN TN TN T 
Mật độ dòng nhiệt: 
   332211332211
2
TbTbTb
A
k
TNTNTN
x
kTN
x
k
x
T
kqx 






 2/20100.50200.50100.0
2500
10
cmW 
   332211332211
2
TcTcTc
A
k
TNTNTN
x
kTN
y
k
y
T
kqy 






 2/20100.0200.50100.50
2500
10
cmW 
Như vậy mật độ dòng nhiêt là không đổi trên toàn miền tam giác. 
Thí dụ 2.14. 
Cho hình vuông phẳng dày 1 m, cạnh 5 cm như trên hình 2.28. Hệ số dẫn nhiệt là 2W/cm0C. Tam giác 
nửa phía trên có nguồn sinh nhiệt trong 1,2W/cm2, nửa dưới tại điểm (1;1) cm có nguồn điểm q* = 
5W/cm theo hướng bề dày. Cạnh đáy được cách nhiệt, cạnh dọc đứng bên phải có nhiệt độ 1000C, cạnh 
đỉnh có đối lưu trong môi trường có nhiệt độ Ta = 30
0C, với hệ số tỏa nhiệt 1,2 W/m2 K. Cạnh đứng bên 
trái có dòng nhiệt q = 2W/cm2. 
Xác định nhiệt độ tại các điểm góc còn lại trong hình. 
Tách hình vuông thành 2 phần tử tam giác như trên hình 2.24. Các phương trình đặc trưng của mỗi phần 
tử có thể thành lập riêng theo các công thức đã biết . 
 Hình 2.28 
 Bảng 2.14 
Nút 1 2 3 4 
Tọa độ 
(m) 
x 0 5 0 5 
y 0 0 5 5 
Phần tử 1 1 2 3 
Phần tử 2 1 3 2 
1 2
3 4 
 
5 cm 
5 cm 
 (1,1) 
q* (5W/cm) 
qV =1,2W/cm
2 
T=100
0
C 
q = 2W/cm2 
h = 1,2 W/cm2 
Ta = 30
0C 
Cách nhiệt 
95 
Công thức chung 
            

  
S
TT
dSNNhdBDBK (2.314) 
         

  
S
T
S
TT
V dSNThdSNqdNqf . (2.315) 
Phần tử 1: Tam giác 1 2 3 
-Tính A, các hệ số ai,bi, ci : 
25
501
051
001
1
1
1
2
33
22
11
yx
yx
yx
A 
 550;550;250.05.5 23132123321 xxcyybyxyxa
000;505;05.00.0 31213231132 xxcyybyxyxa
505;000;00.50.0 32321312213 xxcyybyxyxa 
a1 = 25,0 b1=-5,0 c1 = -5,0 
a2 = 0,0 b2 = 5,0 c2 = 0,0 
a3 = 0,0 b3 = 0,0 c3 = 5,0 
+ Tính [K]1 , do không có đối lưu nên 
  



2
33231
32
2
221
3121
2
1
2
33231
32
2
221
3121
2
1
1
4
1
ccccc
ccccc
ccccc
k
bbbbb
bbbbb
bbbbb
k
A
K yx 



101
011
112
25025
000
25025
000
02525
02525
25.2
2
+ Tính [f]1 ,do có nguồn điểm và bức xạ cạnh bên:        

S
TT
dSNqdNqf 1 
 - Số hạng nguồn điểm trong 
)1;1(



k
j
i
N
N
N
q  , tại điểm (1;1) có 
5
3
25
15
5525
2
1
2
1
1111 yx
A
ycxba
A
N 
96 
5
1
25
5
)050(
2
1
2
1
2222 yx
A
ycxba
A
N 
5
1
25
5
500
2
1
2
1
3333 yx
A
ycxba
A
N 
 Vậy nguồn điểm trong  =1 









1
1
3
1
1
3
5
1
5*
)1;1(k
j
i
N
N
N
q  
- Số hạng bức xạ cạnh bên: 
   dl
N
N
N
qdSNq
S
T  
1
3
3
2
1
Hình 2.29. Phần tử 1 
Trên cạnh 31 có N2 = 0; còn N3 và N1 thay đổi giữa 0 và 1, nên áp dụng công thức tích phân 
)!1(
!!
 ba
ba
dlNN
l
b
j
a
i đối với N1 và N3 thì đều có 
2
1
)!101(
!0!101
31 
l
ji
ll
dlNNdlNdlN 
 Vậy: 
  









5
0
5
1
0
1
2
5.2
1
0
1
2
31qldSNq
S
T 
Véc tơ tải 
      









  
 4
1
2
5
0
5
1
1
3
.*1
S
TT
dSNqdNqf 
- Phương trình đặc trưng của phần tử 1 






4
1
2
101
011
112
3
2
1
T
T
T
 (2.316) 
Phần tử 2: 
1 2 
 3 
q 
97 
+ Tính [K]2 : 
            

  
S
TT
dSNNhdBDBK 2 
Phần tử 2 1 2 3 
Nút 2 4 3 
Tọa độ 
(m) 
x 5 5 0 
y 0 5 5 
Diện tích 
25252525
501
551
051
1
1
1
2
33
22
11
yx
yx
yx
A 
Các hệ số 
550;055;255.05.5 23132123321 xxcyybyxyxa 
505;505;255.50.0 31213231132 xxcyybyxyxa
055;550;250.55.5 12321312213 xxcyybyxyxa 
Tính [K]2 
  



210
120
000
6
.
4
1 23
2
33231
32
2
221
3121
2
1
2
33231
32
2
221
3121
2
1
2
lh
ccccc
ccccc
ccccc
k
bbbbb
bbbbb
bbbbb
k
A
K yx

- Số hạng đầu 






110
121
011
000
02525
02525
25250
25250
000
25.2
2
4 2
33231
32
2
221
3121
2
1
2
33231
32
2
221
3121
2
1
ccccc
ccccc
ccccc
bbbbb
bbbbb
bbbbb
A
k 
- Số hạng sau 
210
120
000
210
120
000
6
5.1.2,1
210
120
000
6
. 23lh 
 Vậy ma trận độ cứng 
98 
  
300
041
011
210
120
000
110
121
011
2K
- Tính véc tơ phụ tải : 
Theo công thức tổng quát, trong tam giác tiêu biểu có 
nguồn trong, bức xạ và đối lưu, hình 2.30, thì đã có 
 









1
1
0
2
.
0
1
1
2
.
1
1
1
3
2312 lhTlqAqf aVe

Hình 2.30. Tamgiác tiêu biểu 
Phần tử 2 có nguồn trong, có đối lưu và không có bức xạ nên rút ra ngay được véc tơ lực sẽ là 
 






1
1
0
2
.
1
1
1
3
43
2
lhTAq
f aV

Thay số với qV =1,2; A=25/2; =1; h=1,2; Ta=30; l34=5 
 















95
95
5
90
90
0
5
5
5
1
1
0
2
1.5.30.2,1
1
1
1
2.3
1.25.2,1
2f 
- Phương trình đặc trưng của phần tử 2 






95
95
5
300
041
011
3
4
2
T
T
T
 (2.317) 
Lắp ghép các phần tử 
Lắp ghép các phương trình tại 4 nút trên như sau 
- Đánh số phương trình ma trận các phần tử 
)3(
)2(
)1(
4
1
2
101
011
112
3
2
1






T
T
T
; ghép với 
)6(
)5(
)4(
95
95
5
300
041
011
3
4
2






T
T
T
- Bảng lắp ghép 
 Bảng 2.15 
2 
4 3 
qV 
Ta = 30 
T = 100 
Hình 2.31. Phầ n tử 2 
99 
Nút 
số 
Phươn
g trình 
Hệ số nhiệt độ Phụ tải 
T1 T2 T3 T4 
1 1 2 -1 -1 -2 
2 2 -1 1 1 
4 1 -1 5 
3 3 -1 1 -4 
6 3 95 
4 5 -1 4 95 
- Lập phương trình ma trận đặc trưng tổng 






95
91
6
2
4010
0401
1021
0112
4
3
2
1
T
T
T
T
- Áp đặt điều kiện biên : biên 24 có nhiệt độ 1000C, tức T2 = 100; T4 =100 thay vào phương trình nút 2 và 
nút 4, phương trình 1 trở thành : 2T1 – T3 = -2 + 100 
Nên dạng ma trận là 






100
91
100
98
1000
0401
0010
0102
4
3
2
1
T
T
T
T
 giải ra 






100
40
100
69
4
3
2
1
T
T
T
T
 (2.318) 
2.14. Dẫn nhiệt hai chiều qua phần tử chữ nhật 
Phần tử chữ nhật có điều kiện biên hỗn hợp thể hiện trên hình 2.32. 
 Hình 2.32. Phần tử chữ nhật 
Phân bố nhiệt độ trong phần tử chữ nhật được viết dạng 
100 
 44332211 TNTNTNTNT (2.319) 
Lấy gốc tại điểm 3 (k), các hàm nội suy sẽ là 
b
x
a
y
N
ab
xy
N
a
y
b
x
N
a
y
b
x
N
2
1
2
4
2
1
2
2
1
2
1
4
3
2
1
 (2.320) 
Ma trận đạo hàm của hàm nội suy là 
   
















)2()2(
)2()2(
4
1
4321
4321
ybxxyb
yyyaya
ab
y
N
y
N
y
N
y
N
x
N
x
N
x
N
x
N
B (2.321) 
Ma trận dộ cứng sẽ là 
            
 S
TT
dSNNhdVBDBK 
trong đó   
y
x
k
k
D
0
0
, và [B]T[D][B] là 
      
)2()2(
)2()2(
4
1
0
0
)2(
)2(
)2()2(
4
1
ybxxyb
yyyaya
abk
k
yby
xy
xya
ybya
ab
BDB
y
xT
Thay vào phương trình trên sẽ được [K] là ma trận 4 4. 
Số hạng điển hình trong ma trận là 
dxdy
ab
xy
dxdyxb
ba
k
dxdyya
ba
k b ab a yb a x
2
0
2
0
2
0
2
0
2
22
2
0
2
0
2
22 4
)2(
16
)2(
16
 (2.322) 
Sau khi tích phân sẽ được 
101 
  
4200
2400
0000
0000
12
2211
2211
1122
1122
6
2211
2211
1122
1122
6
hl
a
bk
b
ak
K
yx (2.323) 
Véc tơ tải là 
   







1
1
1
1
4
2
0
2
0
4
3
2
1
GA
dxdy
N
N
N
N
GdANGf
b aT
 (2.324) 
Tích phân các dòng nhiệt đối lưu và bức xạ tại biên giới được xác định như các phần tử tam giác. 
Thí dụ 2.15. Xác định phân bố nhiệt độ trong tấm phẳng vuông trong ví dụ 2.13. Sử dụng phần tử chữ 
nhật, hình 2.33. 
 Hình 2.33. 
Ma trận độ cứng 
 Ma trận độ cứng theo (2.323) 
  
4200
2400
0000
0000
12
2211
2211
1122
1122
6
2211
2211
1122
1122
6
hl
a
bk
b
ak
K
yx 
Thay số vào được 
102 
  
4200
2400
0000
0000
2211
2211
1122
1122
15
5
2211
2211
1122
1122
15
5
K 
sau khi biến đổi được 
  
20442
42024
4282
2428
6
1
K (2.325) 
Véc tơ phụ tải 
 












1
1
0
0
2
.
1
0
0
1
2
1
1
1
1
4
6 3114
4
3
2
1
* lhTlq
N
N
N
N
qf a



thay các giá trị vào được 
 



3,93
7,97
3,8
7,5
f (2.326) 
Phương trình ma trận đặc trưng cuả phần tử 






3,93
7,97
3,8
7,5
20442
42024
4282
2428
6
1
4
3
2
1
T
T
T
T
 (2.327) 
Áp điều kiện biên, với T2 = T3 = 100
0C, sẽ được phương trình đặc trưng và giải ra nghiệm 






8559
0100
0100
2634
20002
0100
0010
2008
4
3
2
1
, 
, 
, 
, 
T
T
T
T
 -
 - 
 giải ra 



838536
0000100
0000100
484688
, 
, 
, 
, 
T (2.328) 
103 
Trên đây là toàn bộ Phần Phương pháp số trong truyền nhiệt của chương trình cao học cơ khí. 
Các bài toán dẫn nhiệt không ổn định không được dạy trong chương trình cao học. 
Bạn đọc cần phần này để áp dụng tính toán trong nghiên cứu, có thể xem trong cuốn 
 “ Cơ sở Phương pháp Phần tử hữu hạn trong truyền nhiệt “– Nxb Thế giới 
có tại mạng Vinabook. 
Hoặc liên hệ với thày Trịnh Văn Quang tại địa chỉ sau: 
 quangnhiet@yahoo.com.vn