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 - PGS. TS Trịnh Văn Quang

PGS. TS Trịnh Văn Quang

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í

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

1

Mục lục

Lời nói đầu 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 2. Xây dựng hệ phương trình bậc nhất 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 2. Phương pháp tính lặp 2.3. Bài toán dẫn nhiệt không ổn định hai chiều 2.4. Giải hệ phương trình tuyến tính của nhiệt độ 1. Phương pháp định thức 2. Phương pháp Gauss 3. Phương pháp Gauss - Jordan 4. Phương pháp Gauss - Seidel PHẦN 2. PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN

3 4 4 4 5 5 6 9 13 13 13 15 17 20 20 23 23 25 29 36 38 46 47 53 54 59 67 71 75 80 99

2

Giới thiệu khái quát 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 2.6. Các phần tử và hàm nội suy 2.6.1. Phần tử một chiều bậc nhất 2.6.2. Phần tử một chiều bậc hai 2.6.3. Phần tử hai chiều tam giác bậc nhất 2.6.4. Phần tử chữ nhật bậc nhất 2.6.5. Các phần tử đẳng tham số 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 2.7.1. Phương pháp biến phân 2.7.2. Phương pháp Galerkin 2.8. Giải bài toán dẫn nhiệt một chiều bằng phương pháp PTHH 2.9. Dẫn nhiệt qua vách phẳng có nguồn nhiệt bên trong 2.10. Dẫn nhiệt qua vách trụ 2.11. Dẫn nhiệt qua thanh trụ có nguồn trong 2.12. Dẫn nhiệt qua cánh tiết diện thay đổi 2.13. Dân nhiệt ổn định hai chiều dùng phần tử tam giác 2.14. Dẫn nhiệt hai chiều qua phần tử chữ nhật

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ạn…ngoà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).

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

2

2

(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) :

T



T

i

,

j

i

 ,1

j





 T  x

 

T x



x



T

T i

,

j

i

,

j

1





 T  y

 T  y



y

Hình 2.1. Mạng các điểm nút

2

T

T





(

)

(



T

)

)



1

i

i

j

i

,

,

j

i



,1

j

(2.2)





T 2

 ( (  x

T 2 )

 T 2  x 2

(

T



T



T

)

j ( )

)

i

,

j

 1

i

,

i

,

j

i

,

j

 1

(2.3)





x )   T ( 2

 

T 2 y

 ( (  y

T 2 )



y

)

j (

Thay (2.2) và (2.3) vào phương trình vi phân (2.1) sẽ được :



T

)





)

T (



T

)





T

)

T ( i

 1

j

T i

 ,1

j

ji ,

ji ,

 1

ji ,

 1

ji ,

(2.4)





0

T ( 2

T ( 2

ji , (



y

)

ji , (



x

)

 0    T 2 x T  2 y 

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

(2.5)

T







T



T

)

ji ,

T ( i

 ,1

j

T i



,1

j

ji ,

 1

ji ,

 1

1 4

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 .

4

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 :

  ...  C 1

2

2

(2.6)

  ...  C

Ta 11 1 Ta 21 1 ... Ta 12 2 Ta 22 ... ... Ta 1 nn Ta 21 n ... 

n

n

  ...  C Ta 11 n Ta 22 n Ta nn

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



a

T   

 T 2  x

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 :

p

1 

p

T i

T i

(2.8)





 

T 

T   



 1

 1

p

p

2

  Vế phải của (2.8) viết cho thời điểm sau (p+1)  1 p T ( T  1 i i

 1 p  1

(2.9)

 1

p

p

 1

p

 1

  ) T i    ) x  T ( i 2 ) (   T 2 x

 1 p  1

 1 p  1

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

p

1 

p

1 

p

(2.11)





)





( T i

1  p 1 

2 T i

T i

1  p 1 

T i

T i

.   a 2 x ( ) 

Đặt

sẽ được

Fo



  . a 2)  ( x p  1

p

p

 1

(2.12)









)

p  1 1 

T i

T i

TFo .( i

T 2 i

T i

p  1  1

1 

p

p

(2.13)

-FoT

Fo.T

21 







(

Fo)T i

p 1  1 

p 1  1 

T i

i

i

  )  ( T ) 2 (  ) x thay (2.8) và (2.9) vào (2.7): p T i T ( i T i T i T i  a     ) x  T ( i 2 ) (

vậy : 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

5





1

1

p

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  :

(2.14)





q

 p Th K

T 1

h

  - Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian :



1

p



1

(2.15)

q





)



k

p T ( 2

T 1

k  x

Độ tăng nội năng dU phân tố sau thời gian  :

p

 1

p

p

 1

p

(2.16)

dU



 .  c



)



 c



)

TV ( 1

T 1

T ( 1

T 1

 x 2









1

1

1

1

p

p

p



1

p

(2.17)

 c







)



)

  

 p Th K

p T ( 2

T 1

T 1

T ( 1

T 1

Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : k  x  x 2

p

p

p

p

 1

 1

 1

 1

 1

)

2

2











(2.18)



T 1

T 1

T 1

T 1

p ( T 2

 p T K

(

k  c

k  c

  2 x  )

a 

,

,

Đặt Fo =

; Fo là tiêu chuẩn Phuriê, Bi là tiêu chuẩn Biô, a là hệ số khuyếch



Bi

k c

.  xh k

p

 1

 1

 1

 1

 1

p

p

p

p K

p 2

 1

p

p

p

(2.19)

2  Bi.Fo Fo(T   2  T) 1 T 1 T 1

 .1 T 1

Bi.Fo Fo      2 2 2 . TFo 2 T 1

 ...   C 1

2

n

2

(2.20)

  xh 2 ( k x  )  a . 2)  x ( tán nhiệt độ sẽ được :    T T 1 Chuyển nhiệt độ và các đại lượng đã biết sang vế phải   1  1 p 2 . T Bi.Fo K (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 : Ta Ta 11 1 12 2 Ta Ta 21 1 22 ... ...

 ... C  

Ta 1 nn Ta 21 ... ... 

n

1P

(2.21)

Ta ij i

   i  C

trong đó: aij là các hệ số của nhiệt độ phải tìm, iT Ti là nhiệt độ cần tìm ở thời điểm (p+1), viết gọn của 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 :  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):

6

  ...  C Ta 11 n Ta 22 n Ta nn n

1

(2.22)

ij

  i C

là ma trận nghịch đảo của [aii],

 1

    T a i  ija 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 :

p

1 

p

T i

T i

(2.23)





T   

 

T 

 



Vế phải của (2.7) viết cho thời điểm trước p :

p

p

2

p  1

p  1

(2.24)

thay (2.23)và (2.24)vào (2.7) :

p

 1

p

p

p

  ) T ( i T i T i    T 2  x  ( (  x T ) 2 )  ) x  T ( i 2 ) (

p  1

p  1

(2.25)

Để giải (2.25) cần biến đổi như sau :

p

 1

p

p

( 2.26)









)

T i

T i

T ( i

p  1

T 2 i

T i

p  1

  . a 2 )  ( x

Đặt

sẽ được :

Fo



  . a 2)  ( x

p

p

p

 1

)









p  1

p  1

T i

T i

T i

T 2 i

TFo .( i

 1

p

p

(2.27)



21(





.) TFo i

. TFo i

. TFo i

T i

p  1

p  1

Vậy : 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.

7

  ) T i T i T ( i T i T i  a     ) x  T ( i 2 ) (

Để các nghiệm hội tụ cần điều kiện : (1- 2Fo)  0 (2.28) tức là :

Fo 

1 2

hay phải chọn bước thời gian đủ nhỏ :

(

(2.29)



2 ) x a 2

p

(2.30)





q

 

 p Th K

T 1

h

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  : - Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian :

p

(2.31)

q





)



k

p T ( 2

T 1

k  x

- Độ tăng nội năng dU phân tố dày x/2 sau thời gian  :

p

 1

p

p

 1

p

(2.32)

dU



 .  c



)



 c



)

TV ( 1

T 1

T ( 1

T 1

 x 2

- Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên :

p

p

p



1

p

(2.33)





)



 c



)

  

 p Th K

T 1

p T ( 2

T 1

T ( 1

T 1

k  x

 x 2

Hay

p

p

p

 1

p

2





2



)





(2.34)



 p T K

T 1

p ( T 2

T 1

T 1

T 1

k  c

  xh 2 x  ( ) k

k  c

  2 x  )

(

Với Fo =

, Bi =

, a =

sẽ được :

k c

xh . k

 . a 2)  ( x

p

p

 1

p

p



 Bi.Fo T

p K

p 2

p

p

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

(2.35)

2 2    Fo(T   T 1 T)T 1 1 T 1

p 2

p K

8

. T Bi.Fo Bi.Fo  21(  1 FoT Fo   2 2 2 ) T 1 T 1

theo nhiệt độ các

pt 1 1

, tức nghiệm hội tụ cần phải thoả mãn :

pT 1 1

(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 điểm thời điểm trước. Điều kiện để xác định (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.37)

2 

a .

T

a



T   

T  2 x 

 

T 2 y

  

  

- Đ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()

;

T





T

 ax



0

x

;

(2.39)

T





T

 by



0

y

- Điều kiện biên loại 3 :  T  x  T  x

h 1  k h 3  k

 T  x  T  x

h 2 k h 4 k

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

9

2

(

T



T

)



(



T

)

T



T .2



T

)

i

 1

j

i

,

i

,

j

i



,1

j

i



,1

j

i



,1

j

(2.40)







T 2

j 2

j (



x

)

i , x )

(



 

T 2 x

 ( (  x

T 2 )

2



T

)

T



T .2



T

(

T

)





T

)

i

,

j

 1

i

,

j

i

,

j

 1

i

,

j

 1

,

,

i

i

j

 1

(2.41)







j 2

T ( 2

i , y )

(



 

T 2 y

 ( (  y

j (



y

)

T 2 ) Riêng đạo hàm theo thời gian luôn có



T

T i

p 1  , j

p , ji

(2.42)





T   

T   





1











 p j ,

p  ,1

p  ,1

T i

T i

T i

T i

T i

T i

p j ,

p j ,

p j ,





1

1

j

j

(2.43)





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 : p T .2 i j , 2 x )

T .2 i y

p j , 2 )









(

(

 k    . c 

   

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 :

1

1

1











T i

 p j ,

T i

p j ,

T i

p   ,1

1 j

T i

p   ,1

1 j

T i

 p 1 j , 1 

T i

 p 1 j , 1 

(2.44)





2

2





 p T .2 i j , x  )

(

 p T .2 i j , y  )

(

 k    . c 

   

(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ó:









T i

p  ,1

j

T i

p  ,1

j

T i

p , j

 1

T i

p , j



1

1

(2.45)





T i

 p j ,

T i

p , j

.2 T i y



(

p j , 2 )

p .2 T i , j 2 ) x



(

k .    . c

   

   

(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ó:

1

1

t





t

t



.2 t



t

.2 t

p  i

 1 ,1 j

p  1  ,1 j i

p , i

1   j 1

p , i

1   j 1

1

(2.46)



t



t

p , i

 j

p , i

j

2

2

p , i y

 j )

(



p , i  x

(

 j )

k   .  . c

   

  .  

(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

10

, với  là hệ số hấp thụ của vật, IP là năng suất bức xạ chiếu tới

PI. )(q R  



)

 )(q K

P Th ( K

P T m

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à - Điều kiên loại 3 : Dòng đối lưu từ không khí là - Dòng nhiệt tổng :

P

P

(2.47)

 )(q





)



 . I







)





P Th ( K

P T m

P T m

P Th (  K

P T m

 . I h

 P  Th K 

  

trong đó :

P

là nhiệt độ không khí và nhiệt độ bề mặt của kết cấu

K TT ,

P m

h ,  là hệ số toả nhiệt và hệ số hấp thụ của bề mặt

là nhiệt độ quy đổi của bức xạ

I P. h

P

là nhiệt độ tương đương của không khí có kể đến bức xạ



P  T K

P T K

. I h

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.

y 

y 

x 



















  1 p T ,11 j 

p 1  , 1  j

p 1  , 1  j

 T i

 T i

p 1  , j

p 1  , j

p 1  , j

p 1  , j

 p ,1 

T i

T i

T i

T i

1 j

k  x

k  x

k  y

k  y

+ 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   (2.48)

 .c

T





p

 T i   . Tyx i

 p 1 j ,

, ji

1 

1 

1 

1 















 







 1  p Thx K 

  1 p T 1, ji 

p T , ji

p T , ji

p T , ji

p T , ji

 T i

 T i

 1 p ,1 j 

 1 p ,1 j 

k  y

k  x

 y 2

+ Tại biên giới tiết diện, phần tử rộng x, cao y/2, hình 2.3b, có bức xạ và đối lưu tại mặt trên:   

 x  

 c  . x



(2.49)

k  y x  2 p

T , ji

  p 1 T , ji

 y 2

11

1 

1 

1 





















  p 1 Th K 

  p 1 Th K 

 p  1 T ji 1, 

 p 1 T ji ,

p T ji ,

p T ji ,

p T ji ,

1 p  j  ,1

y  2

x  2

x  2

k  y

+ Các phần tử tại góc lồi, hình 2.3c : phần tử rộng x/2, cao y/2, có bức xạ, đối lưu tại 2 mặt lồi ngoài :    T  i 

   



 c



(2.50)

p

y  2   p 1 T ji ,

T , ji

k  x  y x . 2 2

1

1

1

1

y  .









.







  p Th K 

 T i

 T i

 p j ,

 p j ,

p   ,1

 p j ,

p   ,1

T i

T i

T i

1 j

1 j

y  2

k  x

y  2

+ Các phần tử tại góc khuyết trong, hình 2.3d : rộng x, cao y, có đối lưu, bức xạ tại hai mặt khuyết :  

1

1

1

1

1



p



k x  

1  p  j , 1

 p  Th  K

 p j ,

 T i

1  p  j , 1

 p j ,

 Tyx i

 p j ,

j

,

 p j ,

     

(2.51)

Bi



Sau khi lấy x = y , và đặt

, thay vào các phương trình trên sẽ được :

,

Fo





k  c

.  xh k

  2x  

Phương trình tại các phần tử thuộc nút bên trong :

(2.52)

4)Fo

(1

T

T













)



T

( TFo i

p  1  ,1 j

T i

p  1  ,1 j

p , ji

1   1

p , ji

1   1

T i

 p 1 j ,

p , ji

 1

Fo







.2





.2

(2.53)

  41

 TFoBi

  .      c    x . T i T i T i T i k  y  x 2  x 2 3 4 k  y   T  i 

 Fo

 T i

p T 2 ji ,

p  1 j ,1 

p  1 j  ,1

 1  1

T i

p ji ,

p T ji ,

p  1 TFoBi .  K

Phương trình tại các phần tử thuộc nút trên biên :  Phương trình tại các phần tử thuộc nút ở góc lồi :

 1

)





Bi







.4

(2.54)

 Fo

 4

1 p  j ,1 

p 1  T ji 1, 

  p T 1 ji ,

p T ji ,

 p 1 TFoBi .  K

TFo (  2 i Phương trình tại các phần tử thuộc nút ở góc lõm :

)

4

1



















Fo

FoBi .

p j ,

T i

T i

 p 1 j ,

1 p  j ,1 

1 p  j ,1 

T 2 i

T 2 i

 p 1 j ,  1

p 1  j ,  1

TFo ( i

 p 1 TFoBi . .  K

4 3

4 3

  

 T  i 

2 3 (2.55) (2.52), (2.53), (2.54) và (2.55) là các phương trình đặc trưng để tính nhiệt độ tại các nút trong bài toán dẫn nhiệt không ổn định hai chiều, tuỳ thuộc vị trí nút cụ thể trong hình mặt cắt mà các chỉ số i,j được lấy giá trị tương ứng. Từ đó viết lần lượt cho các nút, lập thành hệ phương trình bậc nhất của nhiệt độ. 2.4. Giải hệ phương trình tuyến tính của nhiệt độ Khi nhiệt độ viết dạng hàm ẩn được biểu thị bởi hệ phương trình

12

  ...  C 1

2

2

(2.56)

  ...  C

Ta 11 1 Ta 21 1 ... Ta 12 2 Ta 22 ... ... Ta 1 nn Ta 21 n ... 

n

n

Viết dạng ma trận :

n

a 13 a

23

2

(2.57)



a 11 a 21 ... a

a 12 a 22 ... a

a

... ... ... ...

a 1 a 21 ... a

1 n

n

2

n

3

nn

n

n

     

     

T  1  T  2  ...   T 

      

C  1  C     C 

      

  ...  C Ta 11 n Ta 22 n Ta nn

Các phương pháp giải thông dụng 1. Phương pháp định thức

a 13 a 23

a 13 a 23

D





D 1

... ... ... ...

... ... ... ...

a 12 a 22 ... a n

2

a n

3

a 1 n a n 2 ... a nn

2

a n

3

a 1 n a n 2 ... a nn

a  11  a  21  ...  a  n 1

   ;   

aC  12 1  aC  2 22  ... ...  aC  n n

   ;   

23

a 13 a 23





D 2

D n

... ... ... ...

... ... ... ...

aC 1 13 aC 2 ... aC n n

3

a 1 n a n 2 ... a nn

a 12 a 22 ... a n

2

a n

3

C 1 C 2 ... C n

a  11  a  21  ...  a  n 1

   ;....,   

a  11  a  21  ...  a  n 1

   ;   

Nghiệm sẽ là

(2.58)



;

T



;...;

;

T 1

2

n  T

D 1 D

D 2 D

D n D

2. Phương pháp Gauss Biến ma trận vuông aij thành ma trận “tam giác”. Phép biến đổi ma trận dựa trên nguyên tắc biến đổi hệ phương trình cơ bản quen thuộc sau:

1. Nhân (hay chia) một phương trình với một hằng số thì phương trình đó không đổi 2. Cộng (hay trừ) một phương trình với một phương trình khác trong hệ sẽ được phương trình mới

tương đương với tương với phương trình ban đầu

Thí dụ 2.1 : Cho hệ phương trình (a1), (b1) Hệ ban đầu: hệ 1 2x + 2y = 4 (a1)

áp dụng tính chất 1 với (a1) (a1)/2 x + y = 2 (a2) hệ 2  hệ 1

áp dụng tính chất 2 với (b1) x + y = 2 (a2) hệ 3  hệ 2

13

x + 4y = 3 (b1)

(b1)-(a2) 0 + 3y = 1 (b2)

x + 4y = 3 (b1) (b2)  y = 1/3 ; (a2)  x = 2 - y = 2 – 1/3 = 5/3. Thử lại : (a1) : 2.(5/3) + 2.(1/3) = 12/3 = 4 (b1): 5/3 + 4.(1/3) = 9/3 = 3

Các bước của phương pháp Gauss Hệ ban đầu

n



(1)

n

... ...

1 a 11 1 a 21 1 a 31 1 a n 1

1 a 12 1 a 22 1 a 32 1 a n

2

... ... 1 a 33 1 a n

3

1 a n 1 1 a 2 1 a 3 1 a nn

T  1  T  2    T  n

      

      

      

1  C 1  1 C  2  ...   1 C  n

      

a. Làm các số hạng đầu của mỗi hàng thành 1, bằng cách chia từng hàng cho số hạng đầu tiên của mỗi hàng đó:

... ...

/ /

/ /

/ /

n





n

a a

a a

/ /

/ /

... ...

/ /

a a

1 a 11 1 a 21 1 31 1 n 1

1 a 11 1 a 21 1 31 1 n 1

1 a 12 1 a 22 1 a 32 1 a n

1 a 11 1 a 21 1 a 31 1 a n 1

3

1 a 31 1 a n 1

1 a 33 1 a n

1 a n 1 1 a 2 1 3 1 nn

1 a 11 1 a 21 1 31 1 n 1

T  1  T  2    T  n

      

      

      

1 1  aC / 1 11  1 1 aC /  2 21  1 a ... /  31  1 1 / aC  n n 1

      

...

a

C

1

T T

C

2

n

2 1 2 2



(2)

n

T

a a a

a a a

... ...

a a

C

n

2 n 1 2 2 2 3 2 nn

2 12 2 22 2 32 2 n

... 2 33 2 n

3

2

2 n

      

      

a / a / 2  a      

 1  1   1  1  

      

      

1

...

...

...

a

a

0

)







(3)

)...

...







a

a

0

a

)

)...

...







a

0

a

a

2 a 12 2 a  12 2 a 12 2 a 12

2 a n 1 2 a  n 1 2 a n 1 2 a n 1

2 a 1 n 3 n 2 3 n 3 3 nn

2 a ( 33 2 a () n 3

2 a 12 3 22 3 32 3 n

2 a ( n 2 2 a ( n 3 2 a ( nn

... 3 33 3 n 3

2 a 13 2 a 13

3 n

2

T  1  T  2    T  n

      

      

b. Từ hàng thứ 2, làm các số hạng đầu của các hàng bằng 0, bằng cách lấy các hàng 2, 3...n trừ đi hàng 1 :   2 a (11   22  2 a (11  32  2 a (11    n 2

2  C 1  2 2 CC (   1 2    2 2 CC (   n 1

T   1   T )   2   )    T )    n

2  C 1  3 C  2    C 

 1      

  )     ) 

      

      

c. Từ hàng 2 trở đi , làm các số hạng thứ 2 của mỗi hàng thành 1, bằng cách chia mỗi hàng cho số hạng thứ 2 của hàng đó (tức lập lại bước 1 với hàng 2 trở đi)

14

2 a 12 / a

a

... ...

2 a n 1 / a

a

3 22

2





(4)



/ /

/ /

a a

... ...

/ /

a a

1 0 0 0

2 a 12 1 1 1

... ... ...

a a a

n

3 22 3 a 32 3 a n

2

3 22 3 a 32 3 a n 2

3 a 33 3 a n 3

3 32 3 n 2

3 22 3 32 3 n 2

3 n 2 3 a 3 n 3 a nn

3 n 2

... 4 23 4 33 4 n 3

2 a 1 n 4 a n 2 4 a n 3 4 a nn

4 n

T  1  T     T 

      

T  1  T  2    T  n

      

      

      

2  C 1  4 C  2    C 

      

 1  0   0  0  

      

2  C 1  3 aC /  2    3 aC /  n

      

2

2 C 1 4 C 2



(5)





a a

a a

... ...

1 0 0 0

2 a 12 1 0 0

2 a 12 1  

2 a n 1 4 a n 2 a  a 

... ... ...

a a a

C



n

4 33 4 n 3

4 n 3 4 nn

4 n 2 4 n 2

... 4 a ... 23 4  a 23 4 a  23

2 a 1 n 4 a n 2 5 a n 3 5 a nn

... 4 23 5 33 5 n

4 n

4 2

5 n

3

      

T  1  T  2    T 

T  1  T     T  n

      

      

      

      

      C 

2  C 1  4 C  2    C 

      

      

2

2

(6)







... 4 a ... 23 ...1 ...1

2 a 1 n 4 a n 2 a / / a

2 a 12 1 0 0

2 a 12 1 0 0

a a a a

1 0 0 0

... ...

a a

a a

n

... 4 a ... 23 5 / a 33 5 / a n

2 n 1 4 n 2 6 n 3 6 nn

5 n 3 5 nn

5 33 5 n 3

5 33 5 n 3

5 33 5 n

3

3

T  1  T     T 

T  1  T     T  n

      

      

2  C 1  4 C  2  4 aC /  3  5 aC /  n

2  C 1  4 C  2  6 C ..  3  6 C  n

      

      

      

      

      

, thì sẽ được tam giác sau

1k

nna

d. Làm các số hạng thứ hàng thứ 3 trở đị bằng 0, bằng cách lấy hàng 3, 4.. n trừ đi hàng 2 (tức lập lại bước 2 với hàng thứ 3 trở đi)  1  0   110  110   e. Lập lại bước 1 đối với hàng 3 trở đi …để các số hàng thứ 3 của mỗi hàng trở thành 1  1  0   0  0   g. Tiếp tục như vậy cho đến khi số hạng

1

(7)



0 0 0

2 a 12 1 0 0

... 4 ... a 23 1... 0

2 a 1 n 4 a 2 n 6 a 3 n 1

T  1  T  2  T  3  T  n

      

      

      

2   C 1   4 C   2   6 C ..   3   k C   n

h. Giải ra tính ngược từ dưới lên: hàng dưới cùng :

;

k n





C



T  n C ,.6 3

6 3

T 3

T 3

6 Ta nn 3

6 CTa nn 3

hàng chứa T3 có : 3. Phương pháp Gauss - Jordan Là phương pháp biến ma trận [aij ] thành ma trận đơn vị. Giả sử đã có hệ phương trình ban đầu là ma trận tam giác là

1

...

n



0 0 0

1 a 12 1 0 0

1 a 13 1 a ... 23 1... 0

1 a n 1 1 a 2 1 a n 3 1

T  1  T  2  T  3  T  n

      

C 1     1 C   2   1 C ..   3   1 C   n

      

      

15

1

...

0

a

a

a. Lấy hàng 2 làm gốc, nhân hàng 2 với     T C

12a sẽ được: 2

1 a 12

2 n 2

2 23

2

2

Lấy hàng 1 trừ đi hàng vừa có ở trên



...

a

a

01

...

2 2

1 a 12

1 a 12

1 a 13

1 a 1 n

2 n 2

n

n







2  23 1 a ... 23 1... 0

 1 a 2 1 a n 3 1

10 00 00

2 a 13 1 a ... 23 1... 0

2 a n 1 1 a 2 1 a n 3 1

1 0 0

0 0 0

n

C  1 C 2 1 C .. 3 1 C n

T  1  T  2  T  3  T 

      

T  1  T  2  T  3  T  n

      

      

1  C 1      

      

      

      

2   C 1   1 C   2   1 C ..   3   1 C   n

 01       

(1)

1

b. Lấy hàng 3 làm gốc, nhân hàng 3 với

; Lấy hàng 2 trừ đi

00

a

...

a

2

1 23

2 n 3

    T C 3

3

23a sẽ được:

01

...

a

.

a

a

2 3

1 23

1 2

2 n 3

... 1 23



(2)





2 a 1 n  n 1 a n 3 1

10 00 00

2 a 13 .0 1... 0

2 a 1 n 2 a n 2 1 a n 3 1

10 00 00

2 a 13 a  1... 0

2 C 1 1  C 2 1 .. C 3 1 C n

T  1  T  2  T  3  T  n

      

T  1  T  2  T  3  T  n

      

      

  C     

      

      

      

2   C 1   2 C   2   1 C ..   3   1 C   n

hàng vừa có  01      

c. Tiếp tục như vậy sẽ được

...

(4)



10 00 00

2 a 13 ...0 1... 0

2 a 1 n 2 a n 2 0 1

T  1  T  2  T  3  T  n

      

 01      

      

2   C 1   2 C   2   2 C ..   3   1 C   n

2

2

13a của hàng 1, lấy hàng 3 làm gốc, nhân hàng 3 với

13a , rồi lấy hàng 1 trừ đi kết quả mới

(5)



10 00 00

3 a 1 n 2 a n 2 0 1

n

T  1  T  2  T  3  T 

      

d. Để triệt tiêu có. .. 3  a ...001 14  ...0   1...  0  

      

3   C 1   2 C   2   2 C ..   3   1 C   n

3

3

14a của hàng 1, lấy hàng 4 làm gốc, nhân hàng 4 với

14a rồi lấy hàng 1 trừ đi kết quả mới

e. Để triệt tiêu có.

(6)



0...0 2 a ...010 24 1... 00 00 0

4 a 1 n 2 a n 2 0 1

T  1  T  2  T  3  T  n

      

 01      

      

4   C 1   2 C   2   2 C ..   3   1 C   n

16

Cứ như vậy đến khi hàng 1 chỉ còn số hạng đầu , các số hạng khác đều bằng 0. Tiếp tục làm với hàng 2, 3, ..n g. Cuối cùng có ma trận đơn vị như sau, và có ngay các nghiệm cần tìm



(2.59)





00

0

1

n

n

00...001   00...010   01...000  

     

T  1  T  2  T  3  T 

      

T   1   T   2   .. T   3   T  

k   C 1   k C   2   k .. C   3   k C   n

k   C 1   k C   2   k ... C   3   k C   n

4. Phương pháp Gauss - Seidel Nội dung cơ bản của phương pháp này là cách tính lặp. Phương pháp Gauss- Seidel bao gồm các bước sau. Ban đầu chuyển hệ phương trình nhiệt độ dạng hàm tường cho các nút dạng như sau

   .. )1(;

2 ......

T 1 T    .. )2(: Ta 21 2 Ta 1 12 Ta 3 31 Ta 3 32 Ta n 1 n Ta n n 2

n



n

.1

 1

Lần 1: - Bước 1. Trừ một nhiệt độ tại nút 1 (hoặc nút m nào đó định tính trước tiên), tất cả nhiệt độ còn lại cho

bằng không, thay vào (1) tính ra T1

T   .. a n )(; Ta n 1 1 Ta n 2 2 T nn

- Bước 2. Thay các giá trị T1 mới và T3 = 0, ..,Tn = 0 vào (2) tính ra T2 - Bước 3. Thay các giá trị T1 , T2 mới và T4 = 0, ..,Tn = 0 vào (3) tính ra T3. ...... - Bước n. Thay các giá trị T1 , T2 , ..., Tn-1 mới vào (n) tính ra Tn. Như vậy khi tính được một giá trị nhiệt độ mới phải sử dụng ngay trong các phương trình còn lại . Nghĩa là mọi phương trình luôn phải nhận được giá trị mới nhất nếu có, cho đến phương trình cuối cùng. Lần 2: Lặp lại từ đầu - Bước 1. Thay các giá trị T2, T3, .., Tn vừa có ở lần 1 vào (1) để tính T1 mới. - Bước 2. Thay các giá trị T3, .., Tn của lần 1 đã có và T1 mới vào (2) để tính T2 mới... Tiếp tục như lần 1 đến Tn. Quá trình tính được tính lặp lại lần 3 , lần 4 ...với các giá trị nhiệt độ mới nhất, cho đến khi nào chênh lệch nhiệt độ tại mọi điểm ở hai lần tính sát nhau nhỏ tới mức đủ chấp nhận thì dừng. Thí dụ 2.2 Giải bài toán ổn định hai chiều điều kiện biên loại 1: Một dầm bêtông , tiết diện ngang có hình dạng như hình bên có x=y. Biết nhiệt độ tại các cạnh và góc của tiết diện như trên hình 2.4 . Xác định nhiệt độ tại các điểm bên trong.1,2,3,4,5,6 . Giải : Do x=y , theo (4) các nhiệt trở thành phần của mọi phân tố đều bằng nhau là Rịj =1/ , nên sẽ có :

17

,

T









ji ,

 T i 1

T i

2

T i

3

T i

4

Hình 2.4. Chia mạng tiết diện ngang dầm bêtông

1 4

Tại các điểm 1,2,3,4,5,6 viết được 6 phương trình nhiệt độ dạng hàm tường sau :

T1 = (T2 + 60 + 100 + 50)/ 4 (1) T2 = (T1 + T3 + T5 + 100)/4 (2) T3 = (T2 + T4 + T6 + 100)/4 (3)

T4 = (T3 +100 + 80 +70 )/ 4 (4) T5 = (T2 + T6 + 50 + 40 )/ 4 (5) T6 = (T3 + T5 + 70 + 40 )/ 4 (6)

Bước 1: Thay T2 = 0; T3 = 0; T4 = 0; T5 = 0; T6 = 0 vào (1) tính được T1 = 52,50 Bước 2: Thay T1 =52,5 (giá trị mới) và T3 = 0; T5 = 0 vào (2) tính được T2 = 38,125 Bước 3: Thay T2 = 38,125 vào (3) tính được T3 = 34.5313 Bước 4: .....tiếp tục như vậy sẽ tính được T 4, T 5 , T 6 thứ tự như sau : 52.5000 38.1250 34.5313 71.1328 32.0313 44.1406 Các lần sau : Kết quả tính lặp sau 8 lần viết theo ma trận hàng T = [T1 T2 T3 T4 T5 T6] như sau

(1) 52.5000 38.1250 34.5313 71.1328 32.0313 44.1406 (2) 62.0313 57.1484 68.1055 79.5264 47.8223 56.4819 (3) 66.7871 70.6787 76.6718 81.6679 54.2902 60.2405 (4) 70.1697 75.2829 79.2978 82.3245 56.3808 61.4197 (5) 71.3207 76.7498 80.1235 82.5309 57.0424 61.7915 (6) 71.6875 77.2133 80.3839 82.5960 57.2512 61.9088 (7) 71.8033 77.3596 80.4661 82.6165 57.3171 61.9458 (8) 71.8399 77.4058 80.4920 82.6230 57.3379 61.9575

Bước 6 : Sai số tuyệt đối 2 lần cuối tương ứng là :

0.0366 0.0462 0.0259 0.0065 0.0208 0.0117

là quá nhỏ nên có thể dừng phép tính lặp . Nếu tính theo phương pháp ma trận nghịch đảo , nhiệt độ các điểm tương ứng sẽ là : 71.8630 77.4380 80.5120 82.6310 57.3340 61.9500 Các bài toán thực tế có số nhiệt độ phải tìm lên tới hàng trăm thì phương pháp Gauss -Seidel tỏ rõ rất ưu thế. 5. Phương pháp Ma trận nghịch đảo Hệ phương trình tuyến tính nhiệt độ dạng ma trận :

n

a 13 a

23

2

(2.60)



a 11 a 21 ... a

a 12 a 22 ... a

a

... ... ... ...

a 1 a 21 ... a

1 n

n

2

n

3

nn

n

n

     

     

T  1  T  2  ...   T 

      

C  1  C     C 

      

18

Hay ở dạng gọn sau : [aij]. [Ti] = [Ci] (2.61) Từ đó sẽ rút ra được : [Ti] = [Ci] [aij] - 1 (2.62) trong đó [aij] - 1 là ma trận nghịch đảo của [aij] có dạng :

b 13 a

23

(2.63)



  a

1 ij

... ... ... ...

1

b 12 a 22 ... b n

2

b n

3

b 1 n b 21 ... b nn

b  11  b  21  ...  b  n

     

Các phần tử bịj của ma trận nghịch đảo là phần bù của ma trận chuyển vị của [aịj] . Khi đó nhiệt độ phải tìm sẽ là :

2

1

12

3 + ....

 +.... C+ bC+ bC = b 11

n

23

22

(2.64)



+ ....

2 21 T = b 3 31 1 ......

n

T 1 13 C+ bC+ bC=bT 1 2 3 C+ bC+ bC 32 2 33 3 ...... .... .. ... Cb 1 n Cb 2 n n  Cb 3 n

n

n

3

2

n

Ngày nay nhờ công cụ tính toán hiện đại và các phần mềm tiên tiến nên phương pháp ma trận nghịch đảo được giải rất thuận tiện .

19

 + ....  T n = b C 1 n 1 C+ bC b 2 3 Cb nn

Phần 2. PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN

Giới thiệu khái quát

Phương pháp phần tử hữu hạn (PTHH) là một công cụ số để xác định nghiệm xấp xỉ đối với một lớp rất rộng các bài toán kỹ thuật. Phương pháp PTHH rất được chú ý trong đào tạo kỹ thuật và công nghệ bởi vì nó là một công cụ phân tích có tính đa dạng và mềm dẻo cao. Phương pháp PTHH bắt đầu được hình thành từ nhu cầu giải các bài toán phân tích kết cấu trong lý thuyết đàn hồi trong kỹ thuật công trình và kỹ thuật hàng không. Những người đầu tiên đưa ra phương pháp này là Alexander Hrennikoff (1941) và Richard Courant (1942). Sau Courant đã có nhiều tác giả sử dụng ph- ương pháp rời rạc hoá như Polya, Hersch,Weinberger... tập trung vào nghiên cứu các bài toán giá trị riêng. Từ nửa cuối năm 1950, các tác giả đã phát triển dần hoàn chỉnh phương pháp PTHH. Năm 1959 Greestadt sử dụng nguyên lý biến phân để xác định hàm xấp xỉ trong từng phần tử, và xây dựng các nội dung cơ bản của phương pháp và sau này trở thành lý thuyết toán học của phương pháp PTHH. Các nhà vật lý cũng đã phát triển phương pháp PTHH để áp dụng trong các bài toán vật lý, kỹ thuật như Prager, Synge. Besselinh, Melosh, Fraeijs de Veubeke và Jones đã coi phương pháp PTHH là một dạng của phương pháp Ritz, và là một phương pháp tổng quát nhất để nghiên cứu các bài toán đàn hồi. Họ đã áp dụng cho các bài toán biến phân trong cơ học chất rắn và đã đạt được kết quả khá chính xác. Năm 1965, Zienkiewicz và Cheung đã chứng minh rằng Phương pháp PTHH có thể áp dụng cho tất cả các bài toán của lý thuyết trường, và được công nhận là một phương pháp nội suy rộng. Năm 1973, Fix và Strang đã xây dựng những lý luận toán học chặt chẽ cho phương pháp PTHH, và từ đó nó trở thành một lĩnh vực toán học ứng dụng và được phổ biến và ứng dụng rộng rãi trong kỹ thuật, để xây dựng mô hình dạng số cho các hiện tượng vật lý như trường điện từ và động học chất lỏng… 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 Việc giải các bài toán liên tục bằng phương pháp PTHH luôn được thực hiện theo một trình tự gồm các bước nối tiếp nhau như sau: Bước 1: Rời rạc hóa bài toán , chọn phần tử hữu hạn Miền nghiệm của bài toán, tức vật thể, được chia thành các phần tử có kích thước nhỏ gọi là các phần tử hữu hạn sao cho không có kẽ hở cũng như sự chồng lên nhau giữa các phần tử để bảo đảm tính liên tục của bài toán. Kết quả tạo nên một mạng các phần tử hữu hạn. Tùy thuộc tính chất của bài toán mà chọn phần tử có hình dạng khác nhau: - Với bài toán một chiều, các phần tử được chọn là các đoạn thẳng. - Với bài toán hai chiều, các phần tử được chọn là các hình phẳng như tam giác, tứ giác, chữ nhật… - Với bài toán ba chiều, phần tử được chọn là các hình khối, như khối tứ diện, lập phương, hình hộp, lăng trụ … Mỗi loại phần tử có thể chọn là bậc nhất, bậc hai hoặc bậc ba…tùy theo nhiệt độ phụ thuộc vào toạ độ là hàm bậc mấy. Đặc biệt là trong một loại bài toán có thể dùng các phần tử có dạng khác nhau. Giữa các phần tử ngăn cách nhau bởi biên giới là các nút, đoạn thẳng, hay bề mặt.

20

Hình 2.5. Các dạng phần tử hữu hạn

Tuỳ thuộc loại phần tử mà mỗi phần tử có hai hay nhiều nút. Sau khi rời rạc, nhiệt độ cần phải tìm trong miền liên tục của vật thể được xấp xỉ tại các nút của các phần tử. Bước 2: Chọn hàm nội suy Mối quan hệ giữa nhiệt độ T bên trong phần tử với giá trị nhiệt độ tại các nút Ti được gọi là hàm nội suy Ni (hay hàm hình dạng).

k

(2.65)





... 



TNTNT 11 22

TN k k

TN i i



i

1 

Hoặc ở dạng ma trận

T



..

N



 TN 

(2.66)

 NN 1

2

k

T  1  T   2  ...   T  k

      

ở đây: 1, 2, i…, k là các chỉ số thứ tự các nút trong một phần tử N1 , N2 …Nk là hàm nội suy tại các nút 1, 2…k, và [N] là ma trận hàm nội suy T là nhiệt độ tại điểm bất kỳ trong phần tử T1, T2, Tk tương ứng là nhiệt độ cần tìm tại các nút 1, 2…k ,và [T] là véc tơ nhiệt độ cần tìm. Các hàm nội suy N thường được chọn là các đa thức đại số vì có thể dễ dàng tính đạo hàm và tích phân chúng trong mỗi phần tử. Bậc của đa thức được chọn phụ thuộc vào số các điểm nút của phần tử, đặc điểm và số lượng các ẩn của một nút cũng như yêu cầu liên tục cần có trên biên của phần tử. Bước 3: Thiết lập phương trình đặc trưng của phần tử Phương trình đặc trưng của phần tử biểu thị đặc tính cá thể của các phần tử riêng lẻ, đó là mối quan hệ giữa nhiệt độ chưa biết tại các nút với các phụ tải nhiệt. Để thiết lập phương trình đặc trưng của phần tử, cần thực hiện xấp xỉ hàm cần tìm là nhiệt độ với một số lượng hữu hạn các biến số tại các nút, hình thành một phương trình ma trận của phần tử ở dạng

21

(2.67)



   TK e

e

 e f

là véc tơ phụ tải nhiệt hoặc nhiệt độ cho trước tại nút biên nào đó.

eK là ma trận các hệ số của nhiệt độ, được gọi là ma trận độ cứng của phần tử. f

 ở đây: e là chỉ số biểu thị cho phần tử  eT là nhiệt độ phải tìm tại các nút.   e Một số phương pháp có thể sử dụng để xác định nghiệm xấp xỉ đối với bài toán đã cho là

1. Phương pháp Ritz (tích phân cân bằng nhiệt) 2. Phương pháp Rayleigh Ritz (Biến phân) 3. Phương pháp số dư trọng số

(2.70)

TK

   f 

Nhờ áp dụng một số lý thuyết trong toán học như lý thuyết biến phân, tích phân từng phần, tích phân số và các phép tính ma trận, có thể đưa các phương trình vi phân của bài toán về dạng xấp xỉ (2.67) đối với mỗi phần tử. Bước 4: Lắp ghép các phương trình phần tử để nhận được phương trình tương thích của hệ Để tìm đặc tính của toàn cục của hệ thống, chúng ta bắt buộc phải kết hợp tất cả các phương trình ma trận của các phần tử riêng lẻ, thủ tục đó gọi là lắp ghép các phần tử. Đó là việc tổ hợp các phương trình ma trận của các phần tử riêng lẻ một cách thích hợp để tạo được ma trận đặc trưng trạng thái của toàn bộ khu vực nghiệm của bài toán. Nói cách khác là tập hợp các phương trình vi phân liên tục theo ẩn Te cần tìm ở tất cả các nút của tất cả các phần tử  eT dạng ma trận (2.67) ở trên thành hệ (n phần tử) cũng dưới dạng: 

là tải nhiệt tại các nút của cả hệ

[ K ] là ma trận các hệ số của cả hệ  T là véctơ ẩn của cả hệ  f

Phương trình cho cả hệ (2.70) cũng giống phương trình cho một phần tử chỉ khác là nó có kích thước lớn hơn nhiều. Bước 5: Giải hệ phương trình (2.70) Hệ phương trình (2.70) được giải bằng các phương pháp chuẩn như: Lặp, khử, Gauss, ma trận nghịch đảo... tương tự như giải hệ phương trình trong phương pháp SPHH. Bước 6: Tính các đại lượng thứ cấp. Trong bài toán nhiệt, từ nhiệt độ các nút đã tìm được, có thể tính gradient nhiệt độ, dòng nhiệt theo các hướng, biến dạng nhiệt …

22

2.6. Các phần tử và hàm nội suy 2.6.1. Phần tử một chiều bậc nhất

(2.71)

T



x

1  

2

Trong phần tử bậc nhất, nhiệt độ là hàm bậc nhất của toạ độ: Trong đó 1, 2 là hai tham số cần xác định nên mỗi phần tử cần có hai nút. Gọi hai nút của phần tử là i

và j, có toạ độ là

ix và jx thì nhiệt độ tương ứng tại đó là :

;

(2.72)

T



x

j

1  

2

j

1  

2

i

 x T i

1. Hàm nội suy Mặc dù nhiệt độ tại hai nút vẫn còn là ẩn số phải tìm, nhưng nhiệt độ tại các điểm bên trong phần tử được nội suy theo nhiệt độ hai nút như sau:

(2.73)

T









N



 TN

 N

TN i i

TN j

j

i

j

j

T   i  T 

  

ở đây

jN gọi là các hàm nội suy tại hai nút.

iN và

Từ hai phương trình trong (2.73) giải ra 1, 2 rồi thay vào (2.74), sắp xếp lại sẽ được :

x



x

j

(2.74)

T





j

x



x

x x

 

x i x

j

i

j

i

   

  T  i 

   

  T  

Theo (2.74) các hàm nội suy

là

i NN ,

j

x



x

j

(2.75)

N





N





 N



i

j

x



x

x x

 

x i x

j

i

j

i

   

   

Nếu lấy

, thì

x

 ;0

x



l

i

j

(2.76)

 N



Lấy tổng hàm nội suy

(2.77)

1

N

i N 

j

23

 1  x l x l                  

Từ (2.76) , (2.77) có thể thấy hàm nội suy

jN là hàm bậc nhất theo x, biến đổi ngược chiều nhau có

Nút i 1 0 1

Nút j 0 1 1

x Giữa 0 và 1 Giữa 0 và 1 1

iN và giá trị tại các vị trí khác nhau như sau, bảng 2.1 Bảng 2.1 Hàm Ni Nj Ni +Nj

sẽ được như sau:

 1

N

N

i

j

Từ các kết quả khảo sát trên cho thấy hàm nội suy có hai đặc điểm quan trọng sau: - Hàm nội suy nhận giá trị 1 tại một nút xác định và nhận giá trị 0 tại nút khác. - Tổng của hai nội suy trong phân tố bằng 1 ở mọi vị trí bên trong phần tử, kể cả ở trên biên. 2. Quan hệ giữa biến x với các toạ độ nút Từ (2.76) rút ra toạ độ x ứng với hàm iN rồi thay

i

(2.78)

N

xNx







 N

i

j

i

i

xN j

j

x x

j

   

  

Quan hệ giữa biến x với các toạ độ nút cũng được biểu thị qua hàm nội suy, giống như nhiệt độ. 3. Đạo hàm của hàm nội suy

(2.79)







 N



 

 11

 B

dN dx

d dx

1 l

Đạo hàm của hàm nội suy trong phần tử bậc nhất là hằng số không phụ thuộc x. 4. Gradient nhiệt độ Tuy rằng

là ẩn số chưa biết phải tìm, nhưng trong một phân tố

có giá trị không đổi, nên nhiệt

i TT ,

j

i TT ,

j

độ T trong phân tố chỉ phụ thuộc vào x , vậy gradient nhiệt độ ký hiệu g sẽ là

dN

dN

j

j

(2.80)

g







T





  TB

T i

j

dT dx

dN i dx

dx

dN i dx

dx

j

  

  

T  i  T 

  

Sự thay đổi của các hàm nội suy, nhiệt độ và các đạo hàm bên trong phần tử tuyến tính trên được thể hiện trên hình 2.6. Có thể thấy thay đổi điển hình của nhiệt độ là tuyến tính, đạo hàm của các hàm nội suy là hằng số bên trong mỗi phần tử.

24

Hình 2.6. Sự thay đổi của các hàm nội suy, nhiệt độ và các đạo hàm bên trong phần tử tuyến tính.

Ma trận hàm nội suy [N] và ma trận đạo hàm [B] là hai ma trận rất quan trọng được sử dụng để xác định các đặc tính của phần tử sau này. Thí dụ 2.3. Một thanh dài 12 cm có nhiệt độ tại đầu là 1000C và tại cuối thanh là 1600C. Biết rằng nhiệt độ trong thanh thay đổi tuyến tính. Xác định nhiệt độ tại vị trí cách 8 cm từ đầu thanh. Từ phương trình (2.74)





TNTNT i

j

i

j

0

0

Với :



100

TC ;



160

xC ;



;0

x



12

cm ;

x



cm 8 .

T i

j

i

j

x



x

j

Và :

N







i

x



x

12 12

 

8 0

4 12

j

i

j

j

       

       

thay các giá trị trên vào (2.45) có :

0

= 156,6660C







100



160



666,156

C

TNTNT i i

j

j

4 12

8 12

N    x x    08  0 12 8 12 x i x i

2.6.2. Phần tử một chiều bậc hai

Phần tử một chiều bậc hai nhiệt độ thay đổi theo một chiều, nhưng tỷ lệ với tọa độ theo hàm bậc hai. T(x) = 1 + 2x + 3x2 (2.81)

25

Có 3 tham số 1, 2 và 3 cần xác định nên mỗi phần tử cần 3 điểm là các nút i, j và k phân bố đều trên

phần tử. Trong mỗi phần tử có độ dài

, nếu lấy

thì

.

x

x



l

j

k

k

l  ; 2

Nhiệt độ tại ba nút ứng với các toạ độ là

2

2

;

;

(2.82)







l

l







Tk

 2 3

1

T j

 3 2

1

1iT

l 2

l 2

  

  

Từ (2.82) giải ra các hằng số 1, 2, và 3 rồi thay vào (2.81), sắp xếp lại sẽ được :

2

2

2

(2.83)

)( xT







4



4

2

j

3 x l

2 x 2 l

x l

x 2 l

x l

x 2 l

 T  1 

  

  T    

 T  k 

 1  

l  x  x i 0ix

1. Hàm nội suy Nhiệt độ tại các điểm bên trong phần tử được nội suy theo nhiệt độ tại ba nút như sau:

T i

(2.84)

xT

)(

 TNTN





N

N



 TN 

 N

i

i

j

j

TN k k

j

k

i

j

k

   T   T 

    

Trong đó

và

i NN ,

kN là ba hàm nội suy của phần tử một chiều bậc hai. Từ trên ta thấy các hàm nội

j suy của phần tử một chiều bậc hai là

2

2

2

(2.85)

N



N



N





4



4

2

 N



i

j

k

x 3 l

x l

x  l

x 2 2 l

x 2 l

x 2 l

 1  

  

Từ (2.85) thấy rằng các hàm nội suy là hàm số bậc hai của x, giá trị của mỗi hàm thay đổi phụ thuộc vào toạ độ. Ta có thể lập bảng giá trị các hàm nội suy theo phương trình (2.85) như sau : Bảng 2.2

Nút Hàm nội suy

Ni Nj Nk

Giá trị của hàm nội suy j 0 1 0 1

i 1 0 0 1

k 0 0 1 1

Tổng: Ni+ Nj + Nk

Thay đổi của nhiệt độ và thay đổi hàm nội suy của phần tử bậc hai điển hình được thể hiện trên hình 2.7 như sau:

26

Hình 2.7. Thay đổi nhiệt độ và hàm hình dạng của phần tử một chiều bậc hai.

2. Đạo hàm của hàm nội suy

dN

j

(2.86)





  B

dN dx

dN i dx

dx

dN k dx

3  l

4 l

1  l

x 4 2 l

x 8 2 l

x 4 2 l

 

 

  

  

  

  

  

     

  

  

  

Như vậy đạo hàm của hàm nội suy trong phần tử bậc hai là các hàm bậc nhất của x. 3. Gradient nhiệt độ

Ký hiệu gradient nhiệt độ:

 , g

dT dx

T i

dN

j

(2.87)





j

dT dx

dN i dx

dx

dN k dx

dN dx

 

     TB  T 

  

  

  T   T  k

    

Gradient nhiệt độ là

dN

j

i

(2.88)





T



T i

j

T k

dT dx

dN dx

dx

dN k dx

đạo hàm các biểu thức trong (2.70) thay vào (2.73) sẽ có

(2.89)





8

4

j

x 2

dT dx

3 l

4 l

4 l

x 2 l

1 l

x 2 l

  

 T  i 

  

  T    

 T  k 

Như vậy gradient nhiệt độ cũng như dòng nhiệt phụ thuộc vào toạ độ x. Đạo hàm của hàm nội suy trong phần tử bậc hai là các hàm số phụ thuộc vào biến độc lập x. Ta có thể thấy dùng phần tử một chiều bậc hai sẽ có nhiệt độ xấp xỉ chính xác hơn bậc nhất.

27

Hình 2.8

Thí dụ 2.4 : Xác định giá trị các hàm nội suy của phần tử một chiều bậc hai tại vị trí có toạ độ x = l/4, x = l/3. Tại vị trí có x = l/4 các hàm nội suy là

2

2

= 1-3/4 + 1/8 = 0,3750











Ni

x 3 l

x 2 2 l

l )4/(3 l

l )4/(2 2 l

  

 1  

  

2

2

= 0,75

N



4



4

4



4

4/11

j

x 2 l

l )4/( l

l )4/( 2 l

x l

  

  

  

 1     

k

2

2

l

= -0,1250

 8/14/1



2



2

Nk

x l

x 2 l

4/ l

l )4/( 2 l

  

     

  

Thấy rằng Ni + Nj + Nk = 0,3750 + 0,75 - 0,125 = 1 Tại vị trí có x = l/3 , giá trị của các hàm hình dạng là :

2

2

0,2222











9/211



Ni

x 3 l

x 2 2 l

l )3/(3 l

l )3/(2 2 l

  

 1  

  

2

2

= 0,8889

4



4

 9/43/4

N



4



4

j

x 2 l

l )3/( l

l )3/( 2 l

x l

  

  

  

 1     

2

2

l

= -0,1111



2

 9/23/1



2

Nk

x 2 l

3/ l

l )3/( 2 l

x l

k      

  

  

Cũng thấy ngay rằng Ni + Nj + Nk = 0,2222 + 0,8889 – 0,1111 = 1

28

2.6.3. Phần tử hai chiều tam giác bậc nhất

y

Tk k (xk,yk) 

Ti

 i (x ,y )

 Tj j (xj,yj)

Phần tử hai chiều tam giác bậc nhất là phần tử tam giác có nhiệt độ bên trong phần tử phụ thuộc bậc nhất vào hai chiều tọa độ x và y, được biểu thị bởi T(x,y) = 1 + 2x + 3y (2.90) Phần tử tam giác bậc nhất là phần tử hai chiều đơn giản nhất, nhiệt độ có chứa 3 hệ số.

x

Hình 2.9. Phần tử tam giác bậc nhất trong toạ độ gốc

Do tam giác bậc nhất có 3 nút (hình 2.9), các giá trị của 1, 2 và 3 được xác định từ quan hệ

(2.91)



x



y

;



x



x



y

;



y

T

T

  2

1

j

 3

j

  2

1

k

  2

1

i

 3

i

 3

k

j

k

T i

và

bằng phương pháp định thức

,

x

,

,

, xx i

j

k

TTT j k

i

Có thể giải ra các ẩn là các hệ số 1, 2 và 3 theo như sau. Viết (2.91) ở dạng hệ phương trình:

i

i

(2.92)

x   y 

1

j

j

j

 y  T i T

1   2   2

1

k

k

Ta có

i

i

(2.93)

D





















j

j

 yx i

j

yx j

i

yx k

i

yx i

k

 yx j

k

yx k

j

x x x

y y y

k

k

1   1   1 

    

thấy rằng D = 2A, với A là diện tích của tam giác

x

y

T i

i

i

(2.94)

D

















1

j

j

j

yx j

k

 Tyx j k i

yx k

i

 Tyx k i

j

 yx i

j

 k Tyx i

j

x x

y y

k

k

  T   T  k

    

29

 x y   2 3 x  3  3 T k

i

(2.95)

D













2

j

j

j

 yx k

j

 Tyx i k

j

yx i

k

 Tyx k i

j

  yx j

i

 k Tyx j i

x i x x

y y y

k

T i T T k

k

    

    

i

i

(2.96)

D















3

j

j

j

 yx k

j

 Tyx k

j

i

yx i

k

 Tyx i

k

j

 yx j

i

 k Tyx j i

y y y

x x x

T i T T

k

k

k

    

    

Giải ra

sẽ là







 1

 ; 2

 ; 3

D 1 D

D 2 D

D 3 D

















 1

  yx j

k

 Tyx i j

k

yx k

i

 Tyx k i

j

 yx i

j

 Tyx i

j

k

(2.97)



y





y



y









 



 2

j

 Ty i k

k

 Ty i

j

i

 Ty k j







x









 x

  x

 3

k

 Tx i j

i

 Tx k

j

j

k  Tx i

1 A 2 1 A 2 1 A 2

Với ký hiệu







x-

x-y x a j

b ;y j

k

j

i

(2.98)





k x-

j

i  x a j

k

k i b ; y x- y i j

k

i

k





x c ; y -y k x c ; y - y i j c ; y -y



i x

k x-

b ; y x- y x a j k

k

i

i

j

j

k

i

j

i

(2.97) trở thành









 1

Ta j

j

Ta k

k

 Ta i i

(2.99)









 Tb i i

 2

Tb j

j

Tb k k







 3

Tc j

j

k Tc k

 Tc i i

1 A 2 1 A 2 1 A 2

Thay thế các giá trị của 1, 2 và 3 vào phương trình (2.95) sẽ có

(2.100)

















T





yTc

Ta j

j

Ta k k

 Tb i i

Tb j

j

 xTb k

k

 Tc i i

Tc j

j

k

k

 Ta i i

1 A 2

1 A 2

1 A 2

sắp xếp lại như sau

(2.101)

















T



 a

 a

  a

xb i

 Tyc i i

j

xb j

 Tyc j

j

k

xb k

k  Tyc k

i

1 A 2

1. Hàm nội suy

30

Nhiệt độ tại các vị trí có toạ độ (x,y) bất kỳ trong tam giác được nội suy theo nhiệt độ tại 3 nút của tam giác thông qua hàm nội suy N như sau

T i

(2.102)

,( yxT

)

 TNTN





N

N



 TN 

 N

i

i

j

j

TN k k

j

k

i

j

   T   T  k

    

Từ (2.101) thấy các hàm nội suy là

N







 a



i

i

xb i

yc i

(2.103)

N







 a



j

j

xb j

yc j

N



a







yc

k

k

xb k

k

1 A 2 1 A 2 1 A 2

Viết các hàm nội suy đày đủ theo toạ độ nút:

N









y





 a





 x

y 

i

i

xb i

yc i

  x-yx k

j

j

k

 y -y j

k

 x-x k

j

(2.104)



N











x-

 y x- yx



 x

 x



 

y

 a



j

j

xb j

yc j

i

k

k

i

y - y k

i

i

k



N











y







 a







 x

yx 

k

k

xb k

yc k

  yx i

j

yx j

i

i

 xy j

j

i

1 A 2 1 A 2 1 2 A

1 A 2 1 A 2 1 2 A

(2.105 a)

N

















1

y

Để thấy rõ đặc điểm của các hàm nội suy của phần tử tam giác bậc nhất, chúng ta tính giá trị của chúng tại các nút như sau. - Tính Ni tại nút i có tọa độ là xi và yi 



 x







 yx j i

  yx j

k

 xy k i

yx k

j

k

j

i

i

1 A 2

2 2

A A

- Tính Ni tại nút j có tọa độ

j

(2.105b)



y

N



x

y -

x



y

x-

y





j y x ; 

 y



 x

 0  

j

i

  x-yx k

j

k

j

j

j

k

j

k

j

j

j

1 A 2

k

(2.105c)

N



y





y

x-

y







 x

 0  

- Tính Ni tại nút k có tọa độ là    x-yx k

k

i

j

j

k

j

k

k

k

k

j

k

1 A 2

(2.105d)





x

x y - y



x

y

x-

y

 y x- yx







 

 0  

jN tại nút i có tọa độ là xi và yi 

- Tính  N

k

i

k

i

k

i

i

i

j

i

i

k

i

i

1 A 2

31

x ; k y  xy -xy k

- Tính

j



y

x-

y

k

i

i

j

j

(2.105e)



jN tại nút j có tọa độ  

 N

j

j











 

   1

  y x- yx k   yx i

k i yx j

j y x ;  x x y - y j  yx k

j yx i

i

i

j

k

 x i  yx j

k

k yx k

j

- Tính

jN tại nút k có tọa độ

k

(2.105g)

N





x

x y - y



y

x-

y





 x

   y x- yx

 0  

j

k

i

k

i

k

k

i

k

k

i

k

k

1 A 2

(2.105h)







N













 0  

yx i

j

i

 yx j

i

yx i

i

- Tính Nk tại nút i có tọa độ xi và yi    yx i

yx j

 yx i

k

j

i

i

1 A 2

- Tính Nk tại nút j có tọa độ

j

(2.105h)



N















j y x ; 



 0  

i

k

  yx i

j

yx j

 yx j

i

yx j

j

 yx j

j

yx i

j

j

1 A 2

- Tính

k 





i

i

j

j

k

(2.105k)

N







k

k









  

 

  1

k y kN tại nút k có tọa độ x ;  yx   k  yx k

 yx i  yx i

yx j yx j

i

j

i

yx k yx i

k

 yx j  yx j

k

yx i k yx k

j

Có thể lập bảng giá trị các hàm nội suy tại các nút đối với tam giác bậc nhất như sau Bảng 2.3

Nút i 1

Nút j 0

Nút k 0

0

1

0

0

0

1

iN jN kN

Như vậy ta có thể thấy các hàm nội suy có giá trị bằng 1 ở một nút nhất định và bằng 0 tại tất cả các nút còn lại. Cũng có thể chứng minh được rằng

(2.106)

1

N

N

N





i

j

k

x ; k y

ở tất cả mọi vị trí bên trong phần tử kể cả trên biên giới. 2. Quan hệ giữa biến x,y với các toạ độ nút Từ (2.103) có

32









xN i

i

 xa i

i

xxb i i

yxc i i

(2.107)









xN j

j

 xa j

j

xxb j j

yxc j

j







yxc

xN k

k

 xa k

k

xxb k

k

k

k

1 A 2 1 A 2 1 A 2

Tính tổng





xN i

xN k

k

j

(2.108)























 xN j   xa i

i

xa j

j

xa k

k

 xbx i

i

xb j

j

xb k

k

 xcy i

i

xc j

j

xc k

k 

i 1 A 2

Ký hiệu và tính từng số hạng trong dấu móc đơn của biểu thức trên như sau







a

y





0

 x

y x- yx i

k

k

j

i

j

i

j

k

i

(2.109)









i 

A



b



k

k

i

j

j

i









k 

0



c



 k  xy i  x

   

xa k xb k xc k

k

 x-yx  k j   y -y x j   x x-x k

j

i

 x  y k  x- x i

k

j

  x  y x- yx k i j j   y -y 2 x j i   x x- x j

k

i

xa j xb j xc j

j

 xa i  xb i  xc i

i

Thì sẽ thấy







(

a



bx



cy

)



 20(

Ax



)0

xN i

i

xN j

j

xN k

k

1 A 2

1 A 2

Bởi vậy rút ra

(2.110)

xNx







i

i

xN j

j

xN k

k

Tương tự như vậy, từ (2.107) cũng có











yN i

i

ya i

i

xyb i i

yyc i i

(2.111)









yN j

j

 ya j

j

xyb j j

yyc j

j







yyc

yN k

k

 ya k

k

xyb k k

k

k

1 A 2 1 A 2 1 A 2

Tính tổng







yN i

i

j

yN k

k























yN j   ya i

i

ya j

j

ya k

k

 ybx i

i

yb j

j

yb k

k

 ycy i

i

yc j

j

yc k

k 

1 A 2

Ký hiệu và tính các số hạng trong móc đơn của biểu thức trên như sau

33

a







y

y



y



0

 y

i

j

k

j

k

i

j

k

i

(2.112)

b











 y x- yx j y 0

i 

i

j

k

k

i

j

c











  j y



k 

2

A

 y x- yx i k  yy i  y

   

 ya i  yb i  yc i

i

ya j yb j yc j

j

ya k yb k yc k

k

 x-yx  k j   y y -y j   x-x y k

j

i

  y k  x- x i

 k i j  y -y  i   x- x j

k

j

i

k

Thì sẽ thấy







(

a



bx



cy

)



 00(

x



2

Ay

)

yN i

i

yN j

j

yN k

k

1 A 2

1 A 2

Rút ra

(2.113)

yNy







i

i

yN j

j

yN k

k

Như vậy toạ độ x, y của một điểm bất kỳ luôn thoả mãn mối quan hệ với các toạ độ nút theo hàm nội suy tương tự như quan hệ nhiệt độ tại đó điểm đó với các nhiệt độ nút



xNx





i

(2.114)

i yNy







i

i

xN j j yN j

j

xN k k yN k

k

3. Đạo hàm của hàm nội suy Lấy đạo hàm các hàm nội suy trong (2.103) theo x và y được

 N

j

j



B





 B

 x N 

j

b c

b k c

i

1 A 2

i

j

k

b  i  c 

  

 N  x N   y

 N i  x N   y

 y

 N k  x N  k  y

     

     

     

     

4. Gradient nhiệt độ Gradient nhiệt độ xác định bằng

 N

 N

j

j

i

T i

T

T

T i

j

k









  g

j

 x N 

 x N 

j

j

i

T

T

T i

j

k

  T   T  k

    

 T  x T   y

 N  x N   y

 y

 N k  x N  k  y

 N i  x N  i  y

 y

 N k  x N  k  y

     

     

     

     

     

     

(2.115)

T i

b

j

b k





  TB

j

c

c

1 A 2

i

j

k

b  i  c 

  

k

  T   T 

    

5. Tọa độ khu vực đối với tam giác Hệ tọa độ khu vực hay tự nhiên cũng được sử dụng đối với phần tử tam giác để đơn giản quá trình tính toán. Điểm P ở một vị trí nào đó bên trong tam giác hoàn toàn được xác định bởi ba khoảng cách từ điểm

34

đó đến các cạnh, hình 2.10. Tỷ số giữa các khoảng cách với các đường cao từ đỉnh tương ứng chính là các tọa độ khu vực Li , Lj , Lk.

Hình 2.10

Tọa độ Li được định nghĩa là tỷ số giữa khoảng cách từ điểm P đến cạnh ‘i j’ (tức đoạn PO) và khoảng cách từ điểm i đến cạnh ‘jk’(tức đoạn QR), nghĩa là

(2.116)

Li 

PO QR

Các tọa độ khu vực Lj và Lk cũng được định nghĩa một cách tương tự. Giá trị của Li cũng bằng tỷ số giữa hai diện tích Ai đối diện với điểm ‘i’ và diện tích tam giác toàn phần A, nghĩa là







(2.117)

L i

A i A

.(5,0 .(5,0

PO QR

).( ).(

jk jk

) )

PO QR

Các tọa độ Lj và Lk cũng được tính tương tự có Lj = Aj/A , Lk = Ak/A. Bởi thế tọa độ khu vực còn được gọi là tọa độ diện tích. Do Ai + Aj + Ak = A nên

A

(2.118)







1

A i A

j A

A k A

Nghĩa là Li + Lj + Lk = 1 (2.119) Mối quan hệ giữa tọa độ (x,y) của một điểm bất kỳ trong tam giác với các tọa độ nút được xác định bởi: x = Lixi + Ljxj + Lkxk y = Liyi + Ljyj + Lkyk (2.120) Từ các phương trình (2.118), (2.119) và (2.120) có thể dẫn ra các tọa độ khu vực sau:

35







 a

yc

L i

i

xb i

i

L







 a

yc

j

j

xb j

j

(2.121)

L







 a

yc

k

k

xb k

k

1 A 2 1 2 A 1 A 2

Các hằng số a, b và c trong các phương trình trên cũng được xác định theo phương trình (2.98). Tức là 

x-





b ;y x-y x a k

x c ; y -y k

j

j

i

j

j

k

i

i k x c ; y - y b ; y x- y x a







x-

j

i

k

k

i

i

k



j c ; y -y



i x



k x-

j b ; y x- y x a j k

k

i

i

j

j

i

k

j

i

Nghĩa là tọa độ khu vực cũng chính là hàm nội suy đối với tam giác bậc nhất.

Li = Ni Lj = N Lk = Nk (2.122)

Nói chung các tọa độ khu vực và các hàm nôi suy là như nhau đối với các phần tử tuyến tính, bất kể các phần tử đó là một chiều, hai chiều hay ba chiều. Tích phân hàm nội suy Đối với các phần tử hai chiều tam giác bậc nhất có các tọa độ Li , Lj và Lk chúng ta luôn có công thức đơn giản để tích phân trên toàn tam giác là

(2.123)



2

A

b dALLL j

c k

a i

b dANNN j

a i

c k



 

cba !!!   cba

 !2

A

A

ở đây ‘A’ là diện tích của tam giác. 2.6.4. Phần tử chữ nhật bậc nhất Phần tử tứ giác bậc nhất sẽ có 4 nút như hình 2.11. Phần tử tứ giác đơn giản nhất có dạng hình chữ nhật, trường hợp tổng quát các cạnh hình chữ nhật có thể không song song với các trục toạ độ, hình 2.12.

(x4,y4)

3 (x3,y3)

1

2

(x1,y1)

(x2,y2) x

y Hình 2.12. Phần tử chữ nhật bốn nút

Hình 2.11. Phần tử tứ giác bốn nút

36

Nhiệt độ bên trong tứ giác bậc nhất được đặc trưng bởi phương trình T = 1 + 2x + 3y + 4xy (2.124) Nhiệt độ tại mỗi điểm bên trong tứ giác được nội suy theo nhiệt độ 4 nút T = N1T1 + N2T2 + N3T3 + N4T4 (2.125) trong đó N1, N 2, N 3 và N 4 là các hàm nội suy. 1. Hàm nội suy Để xác định các hàm nội suy, cần viết các giá trị của nhiệt độ tại 4 nút T1, T2, T3 và T4 đối với các toạ độ nút (x1, y1), (x2,y2), (x3,y3), (x4, y4), theo phương trình (2.116), giải ra sẽ nhận được các giá trị 1, 2, 3 và 4. Thay thế các quan hệ này vào phương trình (2.116), sắp xếp lại theo nhiệt độ và so sánh với (2.117) sẽ rút ra các hàm nội suy N1, N 2, N 3 và N 4. Xét trường hợp chữ nhật có gốc toạ độ nằm ở giữa hình và các cạnh song song với hai trục toạ độ, hình 2.13, tức là sau khi đã thực hiện một phép biến đổi chữ nhật trong trường hợp tổng quát ở trên về toạ độ khu vực.

Hình 2.13. Phần tử chữ nhật trong toạ độ khu vực

Khi đó các hàm nội suy N1, N2, N3 và N4 được xác định theo

(

b



ax )(



y

)

N



1

4

(

b



ax )(



y

)

N



2

(2.126)

N



(

b



ax )(



y

)

3

4

N



(

axb

)(





y

)

4

1 ab 1 ab 4 1 ab 1 ab

4

2. Đạo hàm của hàm nội suy Ma trận đạo hàm của hàm nội suy [B] là

37



(

a



y

)

(

a



y

)

(

a



y

)



(

a



y

)

(2.127)





  B

1 ab

4



(

b



x

)



b (



x

)

(

b



x

)

b (



x

)

  

  

N   x N   y

     

     



y

2  

4



x

(2.128)

3  

4

3. Gradient nhiệt độ Từ (2.118), gradient nhiệt độ được viết là  T  x  T  y

Như vậy gradient nhiệt độ trong phần tử thay đổi theo đường thẳng. Do các hàm nội suy là bậc nhất đối với x và y, nên chúng được gọi là có cấu hình song tuyến tính. Các đạo hàm có thể biểu thị như sau









T 1

T 2

T 3

T 4

 N 3  x

 N 4  x

 N 2  x

 T  x

(2.129)

(

a





(

a





(

a





(

a



Ty )



 

Ty ) 1

Ty ) 2

Ty ) 3

4

 N 1 x  1 ab

4

tương tự có



)



(

)



b (





(

 Txb

)

 

 Txb ( 1

 Txb 2

Tx ) 3

4

 T  y

1 ab 4

Ma trận gradient nhiệt độ là

(

a



(2.130)







  g

  TB .

 a ) y   xb ( )

a ( 

 b (

y ) 

x

)

 ) y xb  (

)

 ( ) a y xb  ( )

1 ab 4

  

  

 T  x T   y

     

     

4

T  1  T  2  T  3  T 

      

2.6.5. Các phần tử đẳng tham số 1. Các loại phần tử và hệ tọa độ Phần tử cơ bản, phần tử cong Các phần tử đã khảo sát ở trên có cạnh thẳng gọi là các phần tử thông thường hay cơ bản. Trong thực tế thường gặp các bài toán có hình dạng phức tạp hoặc có biên giới cong. Khi đó cần sử dụng một số lượng

38

rất lớn các phần tử cơ bản có cạnh thẳng dọc theo đường biên giới cong để đạt được đặc tính hình học phù hợp. Trong trường hợp bài toán ba chiều tổng số biến là hết sức lớn và việc giảm tổng số biến là rất quan trọng, đặc biệt khi khối lượng tính toán có liên quan với bộ nhớ máy tính /giá thành. Số phần tử cần thiết trên có thể giảm được đáng kể nếu sử dụng phần tử cong. Có nhiều phương pháp tạo ra phần tử cong, trong đó phương pháp phổ biến nhất được sử dụng là ánh xạ từ các phần tử cơ bản, hình 3.17. Do các hàm nội suy của các phần tử cơ bản đã được biết, có thể viết và xác định trong hệ tọa độ khu vực nào đó, nên các đại lượng đặc trưng của phần tử cong tương ứng cũng sẽ được xác định. Như vậy sẽ có hai hệ thống khái niệm cần được xác định. Một hệ thống là hình dạng các phần tử, hệ thống thứ hai là bậc của các hàm nội suy đối với trường biến. Nói chung không nhất thiết phải sử dụng các hàm nội suy như nhau đối với phép biến đổi tọa độ và phương trình nội suy, và như vậy có hai hệ thống các điểm nút tổng thể khác nhau có thể tồn tại. Hai hệ thống này chỉ đồng nhất trong trường hợp các phần tử là đẳng tham số. Phần tử thực và phần tử quy chiếu Mối quan hệ hàm số của các đại lượng đặc trưng của phần tử (hàm nội suy, đạo hàm,.. tọa độ) trở nên đơn giản, khi sử dụng các phần tử quy chiếu hay phần tử chuẩn hóa. Phần tử ban đầu được rời rạc trong miền khảo sát gọi là phần tử thực. Trong bài toán phẳng, phần tử thực được định vị trong hệ tọa độ gốc (x,y). Phần tử quy chiếu là phần tử đơn giản, định vị trong hệ tọa độ quy chiếu (,) và được dùng để biến đổi thành phần tử thực thông qua phép biến đổi hình học. Để tạo ra phần tử thực từ phần tử quy chiếu, phép biến đổi hình học phải có tính thuận nghịch (song ánh), tức là mỗi điểm trong không gian quy chiếu chỉ ứng với một điểm trong không gian thực và ngược lại. Một phần tử quy chiếu có thể biến thành các phần tử thực cùng loại thông qua các phép biến đổi khác nhau và mỗi phần tử có một phép biến đổi riêng. Bởi vậy phần tử quy chiếu còn được gọi là phần tử “cha-mẹ”.

Hình 2.14. ánh xạ đẳng tham số của tam giác và tứ giác

 TN 

...     TNTNT 11 22 TN mm

Phần tử đẳng tham số Phần tử đẳng tham số là những phần tử có hàm nội suy trường biến đồng nhất với hàm nội suy tọa độ. Trong bài toán nhiệt, trường biến trong phần tử là nhiệt độ được biểu thị bởi hàm số của các nhiệt độ nút

39

N được gọi là hàm nội suy trường biến. Tọa độ trong phần tử được biểu thị bởi hàm số của các tọa độ nút. Hàm số này gọi là hàm nội suy tọa độ

   ....  

   xN  yN 

2

  ...   xNxNx 11 22 yNyNy 11 2 xN mm yN mm

Sự biểu thị nhiệt độ và tọa độ như trên được gọi là biểu thị đẳng tham số, và phần tử như vậy gọi là phần tử đẳng tham số. Nói chung các phần tử đẳng tham số có hàm nội suy nhiệt độ và hàm nội suy tọa độ là đa thức cùng bậc. Các phần tử đã khảo sát ở phần trước đều là đẳng tham số, các phần tử quy chiếu trong hệ tọa độ quy chiếu cũng phải là đẳng tham số. 2. Phần tử đẳng tham số một chiều bậc nhất. Tọa độ quy chiếu đối với phần tử một chiều là tỷ số chiều dài được định nghĩa là: -1   1 , ở đây  là tọa độ quy chiếu. Để chuyển đổi từ tọa độ quy chiếu sang tọa độ gốc, ta thay thế x =  , có gốc là điểm giữa của đoạn thẳng, thay x1 = -1 và x2 =1 vào phương trình (2.75), ta sẽ nhận được







 1



1N

1 2







(2.131)

 1



2N

1   11  )1(   )1(1 

1 2

ở đây i và j là hai nút cuả phần tử một chiều bậc nhất. Vậy hàm nội suy N là

(2.132)





1  

 N





 1

1 2

Nhiệt độ :

,

T





  1

  1

  T 1

2   T

1 2

(2.133)





TNTNT 22 11

, hay



 1

21

N

 1(

N

)



N 1

1

 NN 1

2

2

(2.134)



N



N   11

 22

tức là: Tọa độ: từ (2.131) rút ra:

40

Đạo hàm hàm nội suy [B] theo biến x: Từ (2.132) có



(2.135)

 

11

dN d

1 2

Mặt khác theo đạo hàm hợp thì

J

J 





; với

gọi là Jacobian của phép biến đổi tọa độ. Từ đó suy ra

dN  d

dN dx

dx  d

dN dx

dx d

 1





J



J

  B

 

11

dN dx

dN  d

11  2

3. Phần tử đẳng tham số một chiều bậc hai Phần tử một chiều bậc hai có ba nút, các hàm nội suy theo phương trình (2.85). Tọa độ quy chiếu đối với phần tử một chiều bậc hai là  , với -1   1. Khi chuyển từ tọa độ quy chiếu sang tọa độ gốc, chúng ta thay thế x =  , gốc là điểm giữa của đoạn thẳng, thay x1 = -1, x2 =0 và x3 = 1 vào phương trình (2.85) sẽ được:







 1



iN

 2



  1

2

jN

(2.136)







 1



kN

 2

   0 1        )01   11        )1( 1      10)1(0    )1( 0        01)1(1  

ở đây i, j và k là ba nút của phần tử bậc hai. Để xác định ma trận độ cứng cần tính đạo hàm của hàm nội suy đối với tọa độ gốc x, trong đó tọa độ gốc x là hàm của tọa độ các nút và hàm nội suy

(2.137)







x

Nx i

i

Nx j

j

Nx k

k

Đạo hàm N theo  ở trên viết dưới dạng hàm hợp qua x là

i 

j

(2.138)

j 

k 

dN d  dN d  dN d 

dN i dx dN dx dN k dx

dx  d dx d  dx  d

41



J

với

gọi là Jacobian của phép biến đổi tọa độ. Từ đó suy ra

dx  d

1 J

j

1 J

1 J

(2.139)

dN i dx dN dx dN k dx

dN i d dN j d dN k d

Theo (2.137) thì Jacobian của phép biến đổi phần tử một chiều bậc hai là

dN





x



x



x

  J

i

j

k

dx  d

dN i  d

j  d

dN k  d









2 





 1

 1

  J

 1

 x

j

x k

 2

d  d

d  d

    x i

   2

   

(2.140)





)2( 

x







j

x k

d    d 1   J 2

1 2

  x i 

  

  

  

Đạo hàm của hàm nội suy theo tọa độ gốc được xác định theo phương trình sau:

i

1 

j







  B

  J

dN dx

dN dx dN dx dN k dx

dN i  d dN j  d dN k  d

       

       

        

        

1 

  

(2.141)





)2( 

x







j

x k

1  2

1 2

 x  i 

  

  

     

  





1   2    2  1     2  

         

      

Thí dụ 2.5. Xác định đạo hàm hàm nội suy đối với phần tử một chiều bậc hai với các nút có tọa độ gốc là xi = 2, xj = 4 và xk = 6. Ma trận Jacobian là

42

dN

x



x

x







  J

j

k

i

j  d

dN k  d











  42

1 2

  

 6  

dN dx i  d  d 1     2  2    8 82  



2

1 J

Nghịch đảo Jacobian là  

1 2

Đạo hàm hàm nội suy là

  

  

1 

j

(2.142)









  B

  J

1 2







1  4   1 4

 2    2

1   2     2 1     2  

  

  

         

      

         

      

dN i dx dN dx dN k dx

dN i  d dN j  d dN k  d

       

       

        

        

4. Các phần tử hai chiều Đối với các phần tử là hai chiều, chúng ta có thể biểu diễn các tọa độ x và y là hàm của  và  x = x(,) và y = y(,) (2.143) Để xác định ma trận độ cứng của phần tử, cần biểu thị các đạo hàm hàm nội suy trong tọa độ gốc x,y. Đạo hàm của hàm nội suy viết theo quy tắc hàm hợp như sau

i

, yx









N  i  

N  i  x

x   

N   y

y   

i

(2.144)

, yx









N  i  

N  i  x

x   

N   y

y   

(2.144) có thể viết dạng ma trận

(2.145)

i

i

 

từ đó suy ra đạo hàm hàm nội suy theo tọa độ gốc

43

N  i  x N   y N  i  x N   y N  i   N  i                   J                 x y          x y                   

i

1 

  J

i



Ma trận Jacobian [J]

(2.146)

N   x N   y N  i   N  i                          

  J



Nghịch đảo của ma trận Jacobian [J]-1 được tính theo

      x  y        x  y       

1

(2.147)



  J

 1   det J        y  y         x x        

Các đạo hàm này phải tính được bằng số tại mỗi điểm tích phân, do nghiệm dạng chính xác chưa biết. 5. Phép biến đổi đẳng tham số đối với phần tử tam giác bậc nhất Khi chuyển đổi hệ tọa độ diện tích đối với phần tử tam giác bậc nhất (Li, i = 1, 2, 3) sang tọa độ quy chiếu, biểu thị các hàm nội suy trở nên đơn giản, nếu chọn hệ tọa độ (; ) như trên hình 2.15b, đó là

1

2

   -1 L N 1 L N  1    ;0 2     0; L N 1 (2.148) 3

3

Hình 2.15. Phép biến đổi đẳng tham số của các phần tử tam giác đơn. (a) Tổng thể; (b) Tuyến tính - cục bộ; (c) Bậc hai – cục bộ.

6. Biến đổi đẳng tham số đối với phần tử tam giác bậc hai

44

Đối với tam giác bậc hai có sáu nút, tọa độ quy chiếu được chọn như trên hình 2.15c. Các hàm nội suy tại các góc là N1 = L1(2L1 – 1) = [2(1 -  - ) –1](1 -  - ) N3 = L2(2L2 – 1) = (2 - 1) N5 = L3(2L3 – 1) = (2 - 1) (2.149) Tại các nút giữa cạnh N2 = 4L1L2 = 4(1 -  - ) N4 = 4L2L3 = 4  N6 = 4L3L1 = 4(1 -  - ) (2.150) 7. Phép biến đổi đẳng tham số đối với phần tứ giác bậc hai Đối với phần tử đẳng tham số tám nút, tọa độ quy chiếu được chọn như trên hình 2.16. Nhiệt độ T tại mỗi điểm trong phần tử được xác định bởi

8

(2.151)

T



iTN i



1i 

Hình 2.16. Phần tử đẳng tham số tám nút

Các giá trị tọa độ x và y tại mỗi điểm bên trong phần tử được cho bởi

8

x

 ) ( ,



N

 ). ( , x

i

i



i

 1 8

(2.152)

y

 ) ( ,



N

 ). ( ,

y

i

i



i

 1

ở đây xi và yi là tọa độ của nút ‘i’. Các hàm nội suy của phần tử đẳng tham số quy chiếu là









 1

  1

 1

1N

1 4

45







 1

 2 1  

2N

N









 1

 

 1

1

3



 1

2    1

N 4

N









 1

 

 1

1

5







 1

 2 1  

6N

N







 1

  

 1

1

7

(2.153)



 1

2    1

N 8

1 2 1 4 1 2 1 4 1 2 1 2 1 2

(2.154)

TK

   f 

Các biến  và  là các tọa độ cong và như vậy hướng của chúng sẽ thay đổi theo vị trí. Các nút của phần tử được nhập vào theo trình tự ngược chiều kim đồng hồ bắt đầu từ một góc nào đó. Hướng của  và  được xác định hình 2.16 là  dương theo chiều từ nút 1 đến 3 và  dương theo chiều từ nút 3 đến 5. 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 Phương trình đặc trưng của phần tử là mối quan hệ giữa hàm số cần tìm tại các nút, tức nhiệt độ, và các phụ tải hoặc các lực tương ứng ở dạng ma trận.  Để nhận được phương trình ma trận trên, cần xấp xỉ tích phân phương trình vi phân biểu thị bài toán. Tùy thuộc bài toán mà nhiệt độ biểu thị bởi các hàm số dạng phương trình vi phân khác nhau. - Dạng tổng quát nhất của hàm nhiệt độ trong bài toán dẫn nhiêt là phương trình vi phân dẫn nhiệt

x

y

z

- Bài toán dẫn nhiệt ổn định đối với vật đồng chất đẳng hướng được biểu thi bởi

2

2

2

 k k k q V T      x  T  x  y   T  y  z   T  z                  

- Bài toán dẫn nhiệt ổn định một chiều qua thanh được biểu thi bởi

46

 0     T 2 x T  2 y  T  2 z 

(2.155)

kA





(

)



0

aTThP

2 Td 2 dx

... vv, cùng với các điều kiện biên tùy theo từng trường hợp cụ thể. Các phương trình trên được gọi là phương trình chủ đạo của bài toán. Nghiệm chính xác của bài toán trong nhiều trường hợp không thể giải ra được nên phải tìm nghiệm xấp xỉ. Có nhiều cách tìm nghiệm xấp xỉ của bài toán. Sai phân hữu hạn là phương pháp tìm nghiệm xấp xỉ dạng vi phân, vì nó xấp xỉ vi phân thành số gia, đạo hàm riêng được xấp xỉ thành thương của các số gia và phương trình vi phân được xấp xỉ thành phương trình sai phân. Cuối cùng dẫn tới hệ phương trình bậc nhất của nhiệt độ. Phương pháp phần tử hữu hạn là phương pháp tìm nghiệm xấp xỉ dạng tích phân, vì nghiệm đó là kết quả của việc lấy tích phân phương trình vi phân trong từng phần tử hữu hạn. Nhưng do không thể lấy tích phân trực tiếp phương trình vi phân được, nên phải áp dụng một số lý thuyết trong toán học như lý thuyết biến phân, tích phân hàm trọng số, tích phân từng phần, tích phân số và các phép tính về ma trận... để đưa các phương trình vi phân chủ đạo về dạng xấp xỉ (2.154) đối với mỗi phần tử. Một số phương pháp được áp dụng để xác định nghiệm xấp xỉ tích phân đối với bài toán đã cho là

1. Phương pháp Ritz (tích phân cân bằng nhiệt) 2. Phương pháp Rayleigh Ritz (Biến phân) 3. Phương pháp số dư trọng số

Phương pháp Tích phân cân bằng nhiệt và Biến phân được gọi là phương pháp xấp xỉ tích phân yếu, vì chỉ có thể áp dụng với một số bài toán nhất định. Phương pháp số dư trọng số được gọi là phương pháp xấp xỉ tích phân mạnh, vì có thể áp dụng được với hầu hết các loại bài toán. Phương pháp số dư trọng số gồm có: Phương pháp Collocation (đặt trước giá trị) Phương pháp Sub - domain (miền phụ) Phương pháp Galerkin Phương pháp bình phương nhỏ nhất Trong các phương pháp trên thì Phương pháp Biến phân và Phương pháp Galerkin là hai phương pháp quan trọng nhất, vì chúng có độ chính xác cao nhất và cho kết quả như nhau, nên được ưa chuộng sử dụng trong tính nhiệt. 2.7.1. Phương pháp biến phân Phương pháp biến phân áp dụng lý thuyết quan trọng trong phép tính biến phân phát biểu như sau: “Hàm T(x) sẽ là nghiệm của phương trình vi phân chủ đạo và các điều kiện biên giới, khi nó làm Tích phân biểu thức tương ứng với phương trình vi phân chủ đạo và điều kiện biên của bài toán (gọi là phương trình Euler – Lagrange) đạt cực trị”. Phương trình vi phân chủ đạo trong dẫn nhiệt ổn định là

47

(2.156)

x

y

z

Cùng với điều kiện biên

T = Tb trên mặt S1

k

i



k

j



k

k

q 0

trên mặt S2

x

y

z

T   x

T   y

T   z

k

i



k

j



k

k



(

)



0

trên mặt S3 (2.157)

x

y

z

TTh  a

T   x

T   y

T   z

và k là các pháp tuyến bề mặt, h là hệ số tỏa nhiệt đối lưu, k là hệ số dẫn nhiệt, và q là mật độ

i,

k k k 0 q V  x  T   x  y  T   y  z  T   z                  

j ở đây dòng nhiệt bức xạ. Phương trình Euler- Lagrange Phương trình Euler- Lagrange được phát triển bởi Leonhard Euler và Joseph-Louis Lagrange năm 1750. Đó là công thức cơ bản của phép tính biến phân, được sử dụng để giải các bài toán tối ưu, và kết hợp với nguyên lý hành vi để tính toán đường đi của vật thể trong không gian. Phương trình Euler- Lagrange phát biểu như sau: “Nếu hàm I được cho bởi tích phân dạng

(2.158)

,( yyt ,

dt

)'

I

f



Thì I có giá trị không đổi, nếu phương trình vi phân Euler- Lagrange sau được thỏa mãn

.” (2.160)

Ở đây y’ là đạo hàm của y theo thời gian t:

.

y '

dy dt

Nếu đạo hàm theo thời gian y’ được thay bằng đạo hàm tọa độ yx thì phương trình trên trở thành

(2.161)



f   y

d dx

f  xy 

  

 0 

Phương trình Euler- Lagrange tổng quát đối với ba biến độc lập

(2.162)





f  u 

 x 

f  u 

 y 

f  u 

 z 

f  u 

x

y

z

  

  

  

 0 

   

   

48

 0 f   y d dt f  y  '      

I được gọi là phiếm hàm, đó là biểu thức ở dạng tích phân, chứa các hàm số chưa biết và các đạo hàm của nó. Phiếm hàm I trong bài toán dẫn nhiệt Bằng cách áp dụng phương trình Euler- Lagrange, chúng ta xác đinh được phiếm hàm I tương ứng với phương trình vi phân dẫn nhiệt cùng với các điều kiện biên ở trên là

2

(2.163)



x

y

z

 TTh  a

1 2







2   

2   

2   



S

2

S

2

TI )(  k k k 2  qTds  ds T   x T   x T   x 1 2               dTq  V  

Nguyên tắc và trình tự thiết lập phương trình đặc trưng Chia miền xác định của bài toán  thành ‘n’ phần tử hữu hạn, mỗi phần tử có ‘m’ nút. Nhiệt độ trong mỗi phần tử được biểu thị bằng

m

(2.164)

T



 TN 

TN i i

e  

i

1 

ở đây [N] = [N1, N2 , …, Nm ] là ma trận các hàm hình dạng, và

(2.165)



  T

T  1  T  2  ...   mT 

      

là véc tơ các nhiệt độ nút. Theo nguyên lý biến phân tìm nghiệm xấp xỉ là xác định các giá trị của T để làm I(T) không thay đổi, nghĩa là các giá trị của T thõa mãn biến phân của phiếm hàm I triệt tiêu

n

(2.166)

TI )(



0

 

i

 1

I  iT 

ở đây n là số giá trị rời rạc của T được gán đối với miền của nghiệm. Do Ti là tùy chọn, phương trình (2.166) thỏa mãn chỉ khi

với i = 1, 2,…, n (2.167)

0

I  iT 

Phiếm hàm I(T) có thể viết là tổng của các phiếm hàm riêng lẻ của các phần tử hình thành do việc rời rạc hóa miền nghiệm

n

e

(2.168)



  TI

 e TI





i



1

49

Như vậy thay vì phải xác định phiếm hàm trên toàn miền của nghiệm, ta chỉ cần xác định phiếm hàm của từng phần tử riêng biệt. Từ đó

n

(2.169)

 I

eI 



0

 

i

 1

ở đây Ie được lấy chỉ đối với m giá trị nút liên quan tới phần tử e, nghĩa là

e

e

với i = 1, 2,…, m (2.170)



0

I   T

I  T 

j

  

  

Vì mỗi phần tử có m nút, nên việc giải phương trình (2.170) dẫn tới hệ m phương trình biểu thị đặc tính của mỗi phần tử. Tiếp theo là lắp ghép các phần tử bằng cách cộng toàn bộ các phương trình đặc trưng của tất cả các phần tử theo một nguyên tắc nhất định. Cuối cùng là giải hệ phương trình. Thiết lập phương trình đặc trưng phần tử Tính các số hạng trong phiếm hàm

2

2

2

e

e

e

2

e

e

e

(2.171)

I



k

k



qT

ds





ds

 Th



x

y

z

T a

1 2







T   x

T   x

T   x

1 2



S

2

S

2

  

  

  

  

  

  

  k  

 e  dTq 2 V  

Ma trận gradient nhiệt độ :

e

(2.172)

...

  g

  TB

  ...  T 1 T 2 ...

Ba số hạng đạo hàm đầu:

2

2

2

e

e

e

k

k

k

x

y

z

T   x

T   y

T   z

  

  

  

  

  

  

e

      T  m        ...  N 1  x N  1  y N  1  z  N 2  x N  2  y N  2  z  N m  x N  m  y N  m  z                  T  x e T   y e T   z                  

e

e

e

T

(2.173)

x 0

k 0 0

 gDg   

y 0

z

 0 k  T  y  T  x  T  z       0 k          

T

T



 T  x e T   y e T   z         

   T T B

50

         Theo quy tắc của ma trận chuyển vị của tích   g

T

T

T

(2.174)

gDg



       T

  TBDB

Nên    Nhiệt độ trong phần tử :

e

(2.175)

T



,...,

N



... 

 TN

   

, NN 1

2

m

 TNTN 22

11

TN mm

T  1  T   2  ...   T  m

      

thay thế (2.174) và (2.175) vào phương trình (2.171), ta có

T

T

2

e

(2.176)

I





2







T

ds



   TBDB



   dsTNq

    TNh

a

      T

   dTNq V







1 2

1 2



2 eS

3 eS

0

I e  bây giờ thực hiện cực tiểu hóa tích phân   T 

e

T

T









    dTNq V

      T

    dTBDB





1 2





I    T 

   T



(2.177)

2







ds



0



     TNq ds

    TNh

T a







 1 2



   T

   T

  

   T   

eS 2

eS 3

Vế phải có 4 số hạng dưới dấu tích phân, lần lượt được đánh số (1),(2),(3),(4).





được tính theo (2.175) như

     TN  T 

T e  Đạo hàm của nhiệt độ trong phần tử theo nhiệt độ các nút    T sau

 N 1

2

 N

r

hay

51

 N T e   T 1 T e   T 2 …. e  T  T r

e

(2.178)



N



   T  N

 T    T

N 1 N 2 ... N

m

      

      

Tính riêng từng số hạng của phương trình trên như sau

T

T

T

(2.179)

(1)



      T

    dTBDB

         2 dTBDB





1 2





T

(2.180)

(2)





 

    V2 dTNq

 dNq V





1 2 1 2



   T     T





T

(2.181)



     TNq ds



  Nq

 ds





2 eS

2 eS

 (3)   T 

2

2

2



ds









    TNh

 2

  ds

T a

  TTN a

T a



1 2

1 2

3 eS

3 eS

T

T

(2.182)

 (4)   T 

       dsNTNh

      TNh  ds



    T      TNh a

3 eS

3 eS

Thay các kết quả trên vào (2.177) được

e

T

T

T











        dTBDB

  Nq

 ds

 dNq V







I    T 





eS 2

T

T





0



       dsNTNh

   TNh a

 ds





eS 3

eS 3

T

T

T

T

T

TN









 dNq

 





   ds

  Nq

  Nh

 ds

   TNh

 ds

V

a

        dTBDB

Chuyển vế các số hạng chứa nhiệt độ  T sang vế trái, các số hạng còn lại sang vế phải 









3 eS

3 eS

2 eS

(2.184)

TK

   f 

T

(2.185)

 

  (2.183) Và viết dạng gọn hơn là  trong đó  K

     dBDB

   T dsNNh











3S

52

 

T

T

T

(2.186)

 

  dNq

 Nq



 N



V

a









S

2

S

3

(2.187)

0

f   ds  hT ds

[K] được gọi là ma trận độ cứng của phần tử, nếu phần tử ở bên trong [K] biểu thị nhiệt dẫn. {f} gọi là véc tơ phụ tải nhiệt, chứa các số hạng nguồn trong, bức xạ và đối lưu tại biên. (2.184) là phương trình đặc trưng của phần tử, viết cho một phần tử tổng quát có đủ nguồn trong, bức xạ và đối lưu tại biên giới. Đó là phương trình cốt lõi của phương pháp biến phân trong phương pháp phần tử hữu hạn đối với bài toán dẫn nhiệt. Nếu phần tử không có nguồn sinh nhiệt bên trong (qV = 0), số hạng tương ứng sẽ không có mặt. Tương tự, khi biên giới cách nhiệt (tức q = 0 hoặc h = 0), các số hạng tương ứng cũng không có mặt. 2.7.2. Phương pháp Galerkin Phương pháp Galerkin là một trong các phương pháp số dư trọng số. Phương pháp này yêu cầu biểu thức sau phải thỏa mãn:    dTLk 

là phương trình vi phân chủ

ở đây k là hàm trọng số được chọn là hàm nội suy Nk(x) tại các nút,

 TL

đạo, T là nghiệm xấp xỉ. Điều đó nghĩa là

(2.188)

N

k

k

k

q

d

0

k

x

y

z

V



 x 

T   x

 y 

T   y

 z 

T   z



   

   

   

   

   

   

   

   

Tích phân từng phần được dùng để biến đổi các đạo hàm cấp hai. Khi sử dụng bổ đề Green có thể viết mỗi số hạng đạo hàm cấp hai trong móc vuông thành hai thành phần là

(2.189)

N

k

N

k

ds



k

   dT m

k

x

k

x

x







 x 

T   x

 x 

T   x

N  k  x

N  m  x



S



   

  d  

   

   

ở đây m đặc trưng cho các nút. Với các điều kiện biên giới (2.157), có thể viết (2.189) thành



k



k



k





N

qdS



  dT m

x

x

x

dNq V

k

k







N  k  x

N  m  x

N  m  x

N  k  x

N  m  x

  

  





S

m





dSNhT

0

 TNhN

N  k  x  dS

m

k

a

k





S

S

mT

(2.190) Gộp các hệ số của nhiệt độ nút  lại với nhau sẽ được

(2.191)

K



k



k



k

km

x

y

z

dsNhN k

m





N  k  x

N  m  x

N  k  y

N  m  y

N  k  z

N  m  z



S

  

  d 

Còn lại các số hạng chứa các đại lượng đã cho gồm nguồn trong, dòng nhiệt và nhiệt độ môi trường là

53

(2.192)

k

k

k

a

k









S

S

m

(2.193)

 TK km

  k  f

(2.194)

TK

   f 

f   qN ds  dsNhT dNq V

q



1  2  

sẽ dẫn tới  Tức là  Có thể thấy hai phương trình (2.184) và (2.194) là như nhau, nghĩa là cả hai phương pháp Biến phân và Galerkin cho cùng kết quả như nhau bởi vì có mặt tích phân biến phân kinh điển đối với bài toán dẫn nhiệt. 2.8. Giải bài toán dẫn nhiệt một chiều bằng phương pháp PTHH 1. Vách phẳng một lớp Khảo sát vách phẳng một lớp dày l , hệ số dân nhiệt k , hình 2.17. Phía mặt trái có dòng nhiệt q, mặt phải tiếp xúc với môi trường nhiệt độ Ta , hệ số toả nhiệt tại bề mặt phải là h. Coi nhiệt độ trong vách thay đổi bậc nhất, xác định nhiệt độ hai mặt vách.

Hình 2.17. Vách phẳng và phần tử một chiều tương ứng

Phần tử hữu hạn được chọn là một chiều bậc nhất.Đó là một đoạn thẳng ký hiệu  có hai nút là ‘1’ và ‘2’. a. Ma trận độ cứng và véc tơ phụ tải nhiệt Nhiệt độ hai nút và hàm nội suy tương ứng đã biết là

(2.195)

TN 21

2

(2.196)

và

N



N



1

2

 TNT 21  x  x

x 1 x

 

2

 x x 1

2

x x 1

+ Ma trận độ cứng Ma trận độ cứng của phần tử theo (2.185) đã biết là

54

T

T

(2.197)







     dBDB

dANNh   

e







As

T

,

   dBDB



e1

 

, thì hàm nội suy là :

;0





x

l

x 1

2

 K Trong đó vi phân thể tích d = Adx, diện tích toả nhiệt AS = A diện tích dẫn nhiệt    - Tính số hạng đầu của [K]e :  K trong đó các ma trận [B], [D], [N] xác định như sau Chọn tọa độ

và

(2.198)

N

 1

1

N 2

x l

x l

[D] ma trận hệ số dẫn nhiệt : [D] = k



Đạo hàm của hàm nội suy [B] :   B

 

11

; nên   B T

1 l

 1 l   1  1   

(2.199)

  

  BDB T

 

 11

Vậy số hạng đầu [K]e1

T

lx 

l

1   k   1 1 1 l 1 l k 2l   1  1       1   

T

(2.200)

     dBDB

     dxABDB



x



0

0







.



e

T    dANNh S

2

As

(2.201)

1  1   .  Adx   1 1  1 1 Ak l k 2 l    1       1   







10 



 



22

21

 - Tính số hạng sau của [K]e :  K AS là diện tích toả nhiệt tại mặt phải cũng là A. Mặt khác toả nhiệt xảy ra ở nút 2 nên [N] lấy ở nút 2 , tức là  N Nên

(2.202)

N N  



  dA 10

T    dANNh S

As



A

Vậy ma trận độ cứng của phần tử



Ak l

  

(2.203)



hA



 K e

1  1

00 10

Ak l

  

1    1 

  

  





hA

   Ak l

Ak   l  Ak   l 

  

  

           

     

- Tính véc tơ phụ tải nhiệt {f}.

55

h  hA   0   1     00   10    

T

T

T

Theo (2.186):   f

  dNq

 Nq



 N



V

a









S

2

S

3

Trong đó: - Nguồn nhiệt trong không có nên qV = 0. - Số hạng thứ 2, dòng nhiệt q tại mặt trái, tức nút i của phần tử nên

 [N]

[(N



01 

(N ) i

i

j

[(N



10 

]) i - Số hạng thứ 3, toả nhiệt tại mặt phải, tức nút j của phần tử nên ]) j

i

 [N] (N ) j j Vậy véc tơ phụ tải nhiệt {f} là

T

T

q









ds

ds

dA

qA

ds



  f





 Nq

hT a

 NhT a

AhT a









0

S

S

3

2

A

A

  

1   

  0 

  

0   1 

0   1 

1   

  

  

sẽ là

TK

   f 

  ds  hT ds

qA   AhT  a (2.204) b. Phương trình đặc trưng của phần tử Phương trình đặc trưng 



Ak l

  

(2.205)



j

  

qA   AhT  a

T  i  T 

  





hA

   Ak l

Ak   l  Ak   l 

  

  

           

     

Ví dụ 2.6. Cho bài toán trên với số cụ thể sau: l = 4 cm, k = 0,5 W/m0C, q = 100 W/m2, Ta = 400C, h = 20 W/m2 0C; lấy A = 1m2 , thay số vào được:

12,50

-

12,50



-

12,50

32,50

j

  

  

100   800 

  

T  i  T 

  

Giải ra



T  1  T  2

  

53   

  45 

2. Vách phẳng nhiều lớp Khảo sát vách phẳng có 3 lớp, bề dày và hệ số dẫn nhiệt các lớp tương ứng là l1, l2, l3 và k1, k2, k3. Mặt trái có dòng nhiệt q, mặt phải có môi trường nhiệt độ Ta hệ số toả nhiệt h,hình 2.13. Xác định các nhiệt độ hai mặt ngoài , các chỗ tiếp xúc và dòng nhiệt qua vách. Rời rạc vách thành 3 phần tử và ký hiệu các phần tử và các nút như hình 2.18.

56

1  2  3  4    

q h Ta

l1 l2 l3

L

Hình 2.18. Vách nhiều lớp và cách rời rạc thành các phần tử

a. Phương trình đặc trưng của các phần tử Từ kết quả của một lớp ở trên có thể viết ngay cho từng lớp như sau - Phần tử 1 - (lớp 1) Ma trận nhiệt dẫn và véc tơ phụ tải nhiệt là



(2.206)





 K

;   f

1

1



Ak 1 l 1 Ak 1 l 1

     

     

Ak 1 l 1 Ak 1 l 1 Phương trình đặc trưng của phần tử 1





2

T  1  T 

  

  

qA   0 



Ak 1 l 1 Ak 1 l 1

Ak 1 l 1 Ak 1 l 1

     

     

- Phần tử 2 - (lớp 2) Ma trận nhiệt dẫn và véc tơ phụ tải nhiệt là





 0    qA   

 K



;   f

2

2



Ak 2 l 2 Ak 2 l

Ak 2 l 2 Ak 2 l

2

2

     

     

Phương trình đặc trưng của phần tử 2



(2.207)



0   0 

  

T  2  T  3

  



Ak 2 l 2 Ak 2 l

Ak 2 l 2 Ak 2 l

2

2

     

     

57

 0   0    

- Phần tử 3 - (lớp 3)

3

(2.208)



;   f

 K



3

3

0 hAT a

  

  

 Ak 3 l 

3

3

Phương trình đặc trưng của phần tử 3

  hA Ak 3 l 3 Ak 3 l Ak 3 l                  

3

(2.209)

3

3

 0 Ak 3 l  hAT a T  3  T  4            hA Ak 3 l 3 Ak 3 l Ak 3 l                  

b. Lắp ghép các phần tử Việc lắp ghép các phần tử được thực hiện theo nguyên tắc: cộng các phương trình ở cùng một nút lại với nhau. Để thấy rõ, ta viết các phương trình ma trận của từng phần tử thành các phương trình đại số tại các nút như sau: Từ (2.205)

(nút 1) (2.210)





qA

T 1

T 2

(nút 2) (2.211)





0



T 1

T 2

Ak 1 l 1 Ak 1 l 1

Ak 1 l 1 Ak 1 l 1

Từ (2.206)

(nút 2) (2.212)





0

T 2

T 3

(nút 3) (2.213)







0

T 2

T 3

Ak 2 l 2 Ak 2 l

Ak 2 l 2 Ak 2 l

2

2

Từ (2.207)

(nút 3) (2.214)





0

T 3

T 4

Ak 3 l

Ak 3 l

3

3

(nút 4) (2.215)







T 3

hAT a

Ak 3 l 3

Ak 3 l 3

  

  ThA 4 

Theo nguyên tắc trên thì :

58

- Nút 1: giữ nguyên (2.210):





qA

T 1

T 2

- Nút 2: (2.211) + (2.212):











0

T 1

T 3

Ak 2 l

Ak 2 l

Ak 1 l 1 Ak 1 l 1

Ak 1 l 1 Ak 1 l 1

2

2

- Nút 3: (2.213) + (2.214) :









0

T 2

T 3

T 4

Ak 2 l

Ak 2 l

  T 2  Ak 3 l

Ak 3 l

2

2

3

3

  T 3 

- Nút 4: giữ nguyên (2.215):







T 3

hAT a

Ak 3 l 3

Ak 3 l 3

        

  ThA 4 

Để chuyển trở lại sang dạng ma trận, trong mỗi phương trình trên điền hệ số bằng 0 đối với các nhiệt độ không có mặt, kết quả có dạng ma trận tổng thể như sau

0



0

0





Ak 2 l

Ak 1 l 1 Ak 1 l 1

Ak 1 l 1

2

  

  

(2.216)



0

0





Ak 1 l 1 Ak 2 l Ak 2 l

Ak 2 l

2

2

3

  

  

qA 0 0 hAT a

T  1  T  2  T  3  T  4

      

      

      

0

0





hA

2 Ak 3 l Ak 3 l

Ak 3 l

3

3

  

  

           

           

2.9. Dẫn nhiệt qua vách phẳng có nguồn nhiệt bên trong 1. Giải bằng phần tử bậc nhất Khảo sát vách phẳng dày 2L, hệ số dẫn nhiệt k, nguồn trong qV, nhiệt độ hai mặt như nhau Tm. Trong chương 1 ta đã biết phương trình vi phân đối với bài toán một chiều ổn định có nguồn bên trong là

(2.217)



0



q V k

2 Td 2 dx

59

Nghiệm giải tích của bài toán là hàm bậc 2 của toạ độ, hình 2.19

2

(2.218)

)( xT





x



 2 L



T m

q V 2 k

Để khảo sát bằng phương pháp PTHH, do đối xứng, chúng ta chỉ cần khảo sát một nửa của tấm như trên hình 2.20.

Hình 2.19. Vách phẳng có nguồn trong

eK và  f

a. Rời rạc miền nghiệm Chia nửa tấm thành 4 phần tử, 5 nút, mỗi phần tử dài l = L/8. Coi diện tích mặt cắt ngang truyền nhiệt A = 1 m2. b. Ma trận độ cứng và Véc tơ phụ tải Tính  - Ma trận độ cứng của phần tử theo (2.185), do không có toả nhiệt, nên chỉ còn

Hình 2.20. Rời rạc các phần tử trên nửa tấm phẳng có nguồn trong.

T

(2.219)



  

   dBDB

e

 

 K Hàm nội suy của mỗi phần tử [N] và đạo hàm của hàm nội suy [B] cũng như hai bài toán trước, nên có ngay

T

(2.220)

 K



     dBDB

e

T

(2.221)

f

- Véc tơ phụ tải nhiệt theo (2.186), do không có dòng nhiệt bề mặt và toả nhiệt nên chỉ còn  

 

 dNq V

 



Khi tính véc tơ phụ tải của mỗi phần tử ta lưu ý rằng, do qV phân bố đều trong cả phần tử, nên mỗi hàm nội suy tại mỗi nút lấy giá trị trung bình tại hai vị trí, tức là

1      1 1 kA l    1   

 N



 N



j

j









i

i

i

j

i

j

(2.222)

 N



i

j

 N



Do đó

60

  N N   N  2 2 10  2 01    2          

(2.223)



 N

11 

; 

 TN

1 2

vậy

lx 

T

(2.224)

 1 2 1     1  

  f

  dNq V

x



0







như nhau

eK và  f

  dxA .  q V 1 2 q Al V 2 1     1   1     1  

c. Phương trình đặc trưng của phần tử Phương trình đặc trưng của các phần tử có 



(2.225)







1  1

kA l

q Al V 2

j

j

  

1    1 

1     1  

T  i  T 

  

T  i  T 

  



kA l kA l

Ak i l kA l

q Al V 2 q Al V 2

    

    

    

    

d. Phương trình ma trận tổng thể Lắp ghép các phần tử được ma trận tổng thể của hệ như sau



0

0

0

kA l

q Al V 2







0

0



kA l kA l

kA l

kA l

  

  

(2.226)

0







0





kA l kA l

kA l

kA l

  

  

0

0







T  1  T  2 T  3  T  4  T  5

       

kA l kA l

kA l

kA l

q Al V 2 q Al V 2 q Al V 2

q Al V 2 q Al V 2 q Al V 2

  

  

        

0

0

0



kA l

kA l kA l

                      

q Al V 2

             

             

             

Thí dụ 2.7. Vách phẳng như đầu bài trên với các số liệu cụ thể sau : 2L = 0,06(m); l = L/4 = 0,03/4 (m); k =12 (W/m0C) ; qV = 200000 W/m3 ; nhiệt độ hai mặt ngoài Tm =300C Cho A = 1m2 ; Tính k/l = 1600 m2/W ; q*l/2 = 750 W/m Nếu các phần tử như nhau sẽ có phương trình ma trận đặc trưng tổng thể là



1600 - 1600 0 0

- 1600 3200 1600 - 0

0 1600 - 3200 1600 -

0 0 1600 - 3200

0 0 0 1600

-

0

0

0

-

1600

1600

750

       

T  1  T  2 T  3  T 4   T  5

       

750     1500   1500     1500      

       

61

1600

1600

(2.227)





 750  750

a )3( )3( b

1600

1600



 T 1600 . 4  . T 1600 4

T 1600 . 3  . T 1600 3

Tuy nhiên bài toán cho điều kiện biên tại bề mặt có nhiệt độ Tm =30 = T5 , nên phải áp đặt điều kiện biên tại phần tử 4 như sau : Phương trình ma trận phần tử 3  750   750 

T  3  T  4

  

  

  

  

Phương trình ma trận phần tử 4 cũ

1600



1600

(2.228)





 750  750

a )4( )4( b



1600

1600

1600 T . 4  . T 1600 4

 1600 T . 5  . T 1600 5

  

  

T  4  T  5

  

750   750 

  

để T5 = 30 thì (4b) phải là :

30.

750

a ')4(

30  1600 ')4( b  . T .0 T  5 4  T 1600 4 .

(2.229)

1600 30.  48750 a )'4( 750   

Khi đó ( 4a) sẽ là : Vậy phần tử 4 sẽ có : T 1600 .0 5 T . 5

(4a)’ được cộng với (3b) thuộc nút 4, còn (4b)’ đứng riêng thuộc nút 5. Kết quả được ma trận tổng thể và giải ra nghiệm sau

0

0

750

0

-

1600

1600

0

0

1500

-

-

1600

1600

3200

→ Giải ra

(2.230)





-

- 1600 3200

0 0

0 0

1600 0

3200 - 1600

0

1

30

0

0

       

       

   1500   48750   

T  1  T  2 T  3  T 4   T  5

       

37,0313   36,5625  35,1563   32,8125   30,0000 

       

T  1  T  2 T  3  T 4   T  5

       

       0  So sánh với nghiệm giải tích Bảng 2.4

PTHH

Giải tích

37,0313 36,5625 35,1563 32,8125 30,0000

37,5000 37,0313 35,6250 33,2813 30,0000

Phương pháp T1 T2 T3 T4 T5

b )'4( 30   T . 4 T .0 4

2. Giải bằng phần tử bậc hai Phân bố nhiệt độ của trong vách phẳng có nguồn trong theo giải tích là hàm bậc hai. Vậy có thể khảo sát bài toán lại bài toán trên bằng phần tử bậc hai. Mỗi phần tử bậc hai cần 3 điểm để biểu thị nhiệt độ thay đổi theo hàm bậc hai của toạ độ như đã biết.

62

T = NiTi + NjTj + NkTk (2.231) Các hàm nội suy đã có trong phần trước :



;

;

(2.232)











N k

N j

N i

x l

22 x 2 l

x 4 l

24 x 2 l

x 3 l

22 x 2 l

   

  

  

  

 1  

  

Đạo hàm của hàm nội suy [B] đã xác định được

(2.233)









  B

3 l

4 l

1 l

4 x 2 l

8 x 2 l

4 x 2 l

  

  

  

  

  

     

  

T

, cần phải xác định tích số





     dBDB





 BDB T

a. Ma trận độ cứng Để tính ma trận độ cứng là  K    Trong đó [D] = k; và chú ý các phép nhân ma trận sau:

và

(2.234)

   4.23.1

 43





Nên







k







  

  BDB T

3 l

4 l

1 l

3 l 8 x 2 l

4 x 2 l

8 x 2 l

4 x 2 l

  

  

  

  

  

     

  



1 l

4 x 2 l 4 l 4 x 2 l

        

        

        

        

2











x 2

x 2

x 2

3 l

4 l

8 l

  

  

  

2

(2.235)













x 2

4 l x 2

x 2

4 l

1 l 1 l

x 4 2 l 4 l

x 4 2 l 8 l

  

     

     

           

2











x 2

1 l

3 l 3 l

1 l

4 l

4 l 4 l 1 l

3   l  4 x 2 l 4 x 2 l

4 l 4 x 2 l

   x 4 2 l

8 l

3 l x 8 2 l 4 x 2 l

           

     

  

      8 l   

  

  

  

  

              

        

T



     dBDB



 

+ Tính ma trận độ cứng  K

63

1 43 3  ;11  2 4.23.2 86 4       4.13.1               21     

2

2

2

x

x

x

x





12







9







3

x 2

x 2

x 2

4 x l

16 l

  

  

2

2

12 l 2

   x

x

x

x



12





16





4 



dx



Ak

x 2

x 2



   24 l

8 x l

64 l

1 2 l

x

  

     

     

2

2

2

24 l 32 x 2 l x

x







3



4 





1

x 2

x 2

x 2

12 l

4 x l

16 l    16 l

24 l    8 x l

32 l 8 x l

16 l 16 l

32 l 64 l 32 l

16 l

  

  

  

  

16 l 16 l   

  

           

        

Sau khi lấy tích phân có:

l

3

2

2

3

3

2





9

x





12

x





3

x

40 x 2 l

16 x 2 l

2

3

3

2

3





12

x

x







4

x





Ak

 K



32 x 2 3 l 2 x 64 2 l

32 x 2 l 3

1 2 l

     

     

3

2

2

3

3

2





3

x



4

x







x

24 x 2 l x 32 2 l 3 16 x 2 l

24 x 2 l

8 x 2 l

64 x 2 l 3 32 x 2 l 3

16 x 2 l 3 24 x 2 l 16 x 2 l 3

        

     16    

        

  

  

  16 x   2 l 3     40 x   2 l     16 x   2 l 3  

        

0

Thay cận sẽ được:



12



9

20





12



38

16 3

32 3

16 3





Ak

20





12

16



32



12

4 

 K



Ak 6 l

1 l

32 3

     

14  16 2

16  32  16

2 16  14

    

    



38

12

4 



14

16 3

64 3 32 3

16 3

  

  

        

        

        

   32   3    

              

        

Cuối cùng có ma trận độ cứng của phần tử bậc hai một chiều

(2.236)

14  16

16  32





 K



Ak 6l

2



16

14

    

2   16   

T

, thay [N]T vào sẽ được





  dNq V



b. Véc tơ phụ tải Theo (2.186) :   f

2

1





  

T

(2.237)



Adx

f





q





  dNq

V

V







3 x l 4 x l



2

x  l

2 x 2 l 2 4 x 2 l 2 x 2 l

     

     

           

        

lấy tích phân sẽ được :

64

L

2

3

x







l



l

  

3 2

2 3

x 3 2 l 2

(2.238)

f





l 2









l





Aq V

Aq V

Aq V

l 6

q Al V 6

43



2

3

496    12 8    

    

1   4   1 

    





4 3 l 2 3

l  2

 l       

       

x 4 2 l x 2 l

x 2 2 3 l 3 x 4 2 3 l 2 x 2 3 l

     

     

        

           

0

c. Phương trình ma trận đặc trưng của phần tử

14



16

2

1

T i

(2.239)





j

Ak 6 l

q Al V 6

16 2

32  16

16  14

4 1

    

    

  T   T  k

    

    

    

Thí dụ 2.8. Giải lại với bài toán trên với phần tử bậc hai. Theo đề bài có L = 0,03 m; k =12 W/m0C ; q = 200000 W/m2 ; + Khảo sát bằng 1 phần tử một chiều bậc hai Phần tử có l = L = 0,03 m; A =1 m2. Tính các số hạng : k/6.l = 12/(6.0.015) = 133,3333 ; qVl/6 = 200000.0.015/6 = 500. Thay vào phương trình đặc trưng của phần tử sẽ được

14



16

2

1

14



16

2

1

T i

T i

=









500

j

j

Ak 6 l

q Al V 6

16 2

32  16

16  14

4 1

16 2

32  16

16  14

4 1

    

    

  T   T  k

    

    

    

  33,133   

    

  T   T  k

    

    

    

1866,7

-

2133,3

266,7

500

=



-

2133,3 266,7 -

4266,7 - 2133,3 1866,7

2000 500

    

  2133,3   

T  1  T  2  T  3

    

    

    

Áp đặt điều kiện biên: Do T3 =30, thay vào, hệ trở thành

933,33



1066,7

0

44,1701

-





 giải ra được   T

1066,7 0

2133,4 0

0 1

36001 30

38,9600 30

T 1 T 2 T 3

    

    

    

    

-     

333,33   .   

    

    

+ Khảo sát bằng 2 phần tử một chiều bậc hai Mỗi phần tử có l = L/2 = 0,03/2 = 0,015 m; A =1 m2. Tính k/6.l = 66,6667; qVl/6 = 1000. Nếu các phần tử như nhau, phương trình ma trận đặc trưng của hai phần tử là

65

Phần tử 1

1866,7

-

2133,3

266,7

500



-

2133,3 - 266,7

4266,7 - 1866,7 2133,3

2000 500

    

  2133,3   

T  1  T  2  T  3

    

    

    

Phần tử 2

1866,7

-

2133,3

266,7

500



-

2133,3 - 266,7

4266,7 - 1866,7 2133,3

2000 500

    

  2133,3   

T  3  T  4  T  5

    

    

    

Lắp ghép được

-

1866,7

-



2133,3 4266,7 - 2133,3 0

266,7 2133,3  2133,3

- 1866,7 -

0 0 2133,3 4266,7

-

500 2000  2000

0

0

266,7

-

2133,3

1866,7

500

1866,7   - 2133,3   266,7  0   

0   0   266,7  2133,3   

T  1  T  2 T  3  T 4   T  5

       

   500     

   500     

Áp đặt điều kiện biên T5 =30, hệ trở thành

-

2133,3 4266,7

266,7 2133,3

-

0 0

0 0

500 2000

 T=



-

-



266,7 0 0

2133,3 0 0

3733 4, 2133,3 - 0

2133,3 4266,7 0

266,7 0 1

1000 2133,3.30 30

1866,7   - 2133,3      

       

T  1  T  2 T  3  T 4   T  5

       

     2000   

      65999   

37,4732     37,0070   35,6051     33,2705     30,0000  

So sánh với nghiệm giải tích, phần tử bậc nhất, bậc hai một phần tử và hai phần tử như sau Bảng 2.5

Nghiệm giải tích

PTHH bậc hai

PTHH bậc nhất (5 phần tử) 37,0313 36,5625 35,1563 32,8125 30,0000

1 phần tử 44,1701 38,9600 30,00

2 phần tử 37,4732 37,0070 35,6051 33,2705 30,0000

37,5000 37,0313 35,6250 33,2813 30,0000

T1 T2 T3 T4 T5

Thấy rằng nghiệm PTHH khi dùng hai phần tử bậc hai chính xác hơn 5 phần tử bậc nhất

66

2.10. Dẫn nhiệt qua vách trụ Xét vách trụ đường kính trong d1, ngoài d2 , hệ số dận nhiệt k, mặt trong có nhiệt độ Tm1, mặt ngoài toả nhiệt ra môi trường hệ số toả nhiệt h, nhiệt độ môi trường Ta. Khi sử dụng phương pháp phần tử hữu hạn, có thể coi thay đổi nhiệt độ là tuyến tính. Chọn 1 phần tử một chiều bậc nhất, chiều dài phần tử là bề dầy vách l = r2 – r1, hình 2.21.

Hình 2.21. Vách trụ và chọn phần tử một chiều tương ứng

2- r1

(2.240)





TNTNT 11 22

2)1, vi phân thể tích là d = 2rdr. Như vậy biến số độc lập Thể tích phần tử khảo sát là  = (r2 trong vách trụ là r thay cho x trong vách phẳng. Nhiệt độ trong vách trụ vẫn tuân theo các công thức của phần tử một chiều bậc nhất, được nội suy qua nhiệt độ hai nút 1. Hàm nội suy Khi đặt r1 = 0; r2 = l , các hàm nội suy [N] đối với vách trụ cũng giống như đối với vách phẳng sẽ là

(2.241)





 N

 N

2

1

(2.242)



   N r l               1    

11

 

1 l





rNrNr 11 22

r l 2. Đạo hàm của hàm nội suy Đạo hàm của hàm nội suy [B] cũng như trong vách phẳng   B

Toạ độ r được biểu thị qua hàm nội suy 3. Ma trận độ cứng Ma trận độ cứng của phần tử vách trụ vẫn theo công thức

67

T

(2.243)







     dBDB



   T dANNh s







As

 K

T

+ Tính số hạng thứ nhất :

:

   dBDB

  



đối với phần tử một chiều bậc nhất, ta đã biết là

 BDB T

tích số   

  

  BDB T

 

 11

và với vách trụ ở đây l2 = (r2 – r1)2 ; nên

r

2

2

T

r

2

1



1



rdr



     dBDB

2

ri



 

 2 k 2 l



1

1



1

1

r 2



 2 k  r 1

 r 2

  

1   

  

1   

r

1

Sau khi thay cận có

T

1   k   1 1 1 l 1 l k 2l   1  1       1   



 r 1

(2.244)

     dBDB



- Tính số hạng thứ hai :

:



   T dANNh s

As

Diện tích toả nhiệt mặt ngoài vách trụ AS = 2r21. Toả nhiệt chỉ ở nút 2 đã tính trong (2.202), nên có

1  r 2   2 k l  2  1 1    1   

(2.245)



T    dANNh S

2hr

As



A

Vậy ma trận độ cứng [K] là

0 00 h   2   10   1     dA 10        



 r 1

(2.246)

 K



1  00 r 2   .2 hr 2  1 1 10  2 k l  2    1         

4. Véc tơ phụ tải nhiệt Do chỉ toả nhiệt tại mặt ngoài diện tích AS = 2r21, nên

T

 



 NhT a

As

f  dA S hT a  .2 2r   0   1    

5. Phương trình đặc trưng của phần tử Cuối cùng có phương trình đặc trưng phần tử là

68

1



00



 r 1

r 2

(2.247)

 .2



hT

hr 2

a  .2 r 2

 k 2 l

 2



1

1

10

  

1   

  

  

T  1  T  2

  

0   1 

  

  

  

Thí dụ 2.9. Tính nhiệt độ mặt ngoài và phân bố nhiệt độ trong vách trụ với số liệu sau: r1 = 40 cm, r0 = r2 = 60 cm, k = 10W/m0C, Tm1 =1000C, h = 10W/m2 0C, Ta = 300C. + Khảo sát bằng sơ đồ một phần tử Chiều dài phần tử l = r2 – r1 = 60 – 40 = 20 cm. Ma trận độ cứng và véc tơ tải như sau

(

)

r 1

r 2



 .2

  K e

hr 2

1  1

00 10

k  2 l

 2

  

1    1 

  

  



.

 10.6,0.2

1  1

00 10

50  50

50  62

 10.2 2,0

)4,06,0(  2

  

1    1 

  

  

   

  

 

a

Phương trình ma trận đặc trưng của phần tử là

50



0

(2.248)



50

62

   

50   

T     1       T 360     2

Áp đặt điều kiện biên : T1 = 1000C sẽ có

0

100

 T2 = 86,450C

0

62



100.50

2

1    

  

T    1      T 360   

  

0 f  hT  6,0.2.30.10   r .2 0 0   1     0   1        360    

là khá lớn so với nghiệm chính xác là 86,300C. + Khảo sát bằng sơ đồ hai phần tử bậc nhất Khi coi bề dày vách trụ gồm hai phần tử, sơ đồ sẽ có ba nút:1, 2 và3. Chiều dài mỗi phần tử là: l = (r2 – r1)/2 = (60 – 40)/2 = 10 cm, ba nút tương ứng với các toạ độ là: r1 = 40 cm, r2 = 50 cm và r3 = 60 cm. + Phần tử 1: Phần tử 1 có hai nút 1 và 2, không có đối lưu - Ma trận độ cứng



 r 1

 K



1

69

1  1  90  90 r 2   1 1  1 1  90 90 k .2  l  2  10.2 1.0 5,04,0  2    1       1          

- Véc tơ tải   1f

90



- Phương trình ma trận đặc trưng :





90

90

   

90   

T  1  T  2

  

0   0 

  

+ Phần tử 2: Phần tử 2 có hai nút là 2 và 3, có đối lưu tại nút 3 - Ma trận độ cứng :

(

)

r 2

r 3



 .2

 K



2

hr 0

1  1

00 10

k  2 l

 2

  

1    1 

  

  



.

 10.6,0.2

1  1

110  110

110  122

 10.2 1,0

)6,05,0(  2

  

1    1 

00     10  

   

  

- Véc tơ tải

 0   0    

 

2

110



110

0

- Phương trình ma trận đặc trưng:



110

122

   

  

T    2      T 360    3

  

+ Lắp ráp phương trình đặc trưng tổng thể :

90



90

0

0

(2.249)



110



90 0

90 

 110

     

  110    122 

T  1  T  2  T  3

    

     0     360  

Hay gọn lại là

90



90

0

0

(2.250)





90 0

200  110

110  122

    

    

T  1  T  2  T  3

    

    0     360  

Áp đặt điều kiện biên : do T1 = 1000C, nên

0

0

100

100

 giải ra

(2.251)





100



0 0

200  110

110  122

.90 360

1     

    

T  1  T  2  T  3

    

  0   

    

T  1  T  2  T  3

    

  ,92   ,86 

  4878   3415 

70

0 f   6,0.2.30.10   hT a r .2 3 0   1     0   1        360    

2.11. Dẫn nhiệt qua thanh trụ có nguồn trong Xét thanh trụ đường kính trong d1, ngoài d2 , hệ số dận nhiệt k, mặt trong có nhiệt độ Tm1, mặt ngoài toả nhiệt ra môi trường hệ số toả nhiệt h, nhiệt độ môi trường Ta, bên trong vách có nguồn qV 1. Ma trận độ cứng Khi phần tử có nguồn bên trong, phương trình ma trận độ cứng (2.220) vẫn không thay đổi



 r i

j

(2.252)

 K



0hr

 r 1  00   .2  1 1 10 2  2 k l    1         

2. Véc tơ phụ tải nhiệt

T

Véc tơ phụ tải, ngoài số hạng đối lưu sẽ có thêm số hạng nguồn trong

.2 , nên rdr



 Nq V



r

T

T

.2 (2.253)

rdr





  f





 NhT a

dA S

 Nq V

As





r

T

- Số hạng đối lưu đã biết là



 NhT a

As



f

qV

T

rdr

f

.2 (2.254)

- Tính số hạng nguồn trong ký hiệu  





 Nq V

qV



r

Biến số độc lập r trong tọa độ trụ được biểu thị bằng

(2.255)





r

rN j

j

rN i i

Trong đó Ni và Nj là các hàm nội suy:

(2.256)

N

1 

;

N



i

j

x l

x l

2

thay (2.229) vào (2.228) sẽ được N

i

i

i

i

(2.257)

0  dA S hT a  .2 2r   1    

  f

qV

i

i

j





j

j

i

j

j

r

r

Để tính biểu thức trên, cần áp dụng công thức tích phân :

dlNN



, (2.258)

a i

b j

l



ba !!  ba

(

)!1

71

 ) dr  dr  2 q V rNrN  j  2 q V N rNNrN  j j 2 rNrNN  i     .(          

Với

là hàm nội suy cũng là các toạ độ khu vực; a, b là các số mũ. Các số hạng

i NN ;

j

dlNN









và

(2.259)

1 i

1 j

2 dlN i

l

l



!1!1  )!111(

1 6

!0!2  )!102(

l 3

Thực hiện tích phân số hạng nguồn trong (2.257)

rj

j

j

j

(2.260)

  f

qV

rj ri

j

j

j

ri

Vậy véc tơ lực của phân tố là

2



r

0

j

T

T

V

(2.261)





 .2

rdr





l

  f





 NhT a

dA S

 Nq V

hT a

 .2 r 2





As

r i 

2

r

 q 2 6

r i

j

r

  1 

  

  

  

r

2



r



1





j

T i

j

 r i

(2.262)

l





 .2

hr 0

hT a

 .2 r 2

r i 

2

r

1

2

1



 2 k l

 q 2 V 6

j

r i

j

  

  

  

0   1 

1   

00   10 

  

  T 

  

  

  

 r r i 2  r 2  r   r .  l  q 2 V r i  2 r r i  2 r  q 2 V 6  2 q V 6 r i r i              r r i r 6 r 3 r   3.  r  6      

    1 2 3 4 5     

h = 4000W/m2 Ta = 200C

Tâm thanh trụ Mặt ngoài

3. Phương trình ma trận đặc trưng Phương trình ma trận đặc trưng của phân tố đối vách trụ có nguồn trong là    Thí dụ 2.10. Xác định nhiệt độ trong thanh trụ rất dài bán kính 25 mm, có nguồn nhiệt thể tích 35,3MW/m3, hệ số dẫn nhiệt 21W/m0C. Mặt ngoài tiếp xúc với chất lỏng nhiệt độ 200C, hệ số tỏa nhiệt 4000W/m2 0C. Chúng ta chia nửa miền khảo sát là bán kính thành 4 phần tử, mỗi phần tử dài 6,25 mm như trên hình 2.22.

Hình 2.22. Rời rạc phần tử hữu hạn trong thanh trụ dài vô hạn

Toạ độ các nút r1 = 0; r2 = 0,00625; r3 = 0,0125 ; r4 = 0,01875; r5 = 0,025 + Ma trận độ cứng của các phần tử: - Các phần tử 1, 2, 3 không có đối lưu nên có công thức chung như nhau

72



 r i

j

(2.263)

 K



với



l



r

 r 1    1 1 2    1   

j

1



1



5,10



5,10

,0

r 1

r 2



 K



1



1

1



1

1



5,10

5,10

 21.2 00625 ,0

 2

 21.2 00625 ,0

00625 2

  

1   

  

1   

  2  

  





5,31



1

1

,0

,0

0125

r 2

r 3



 K



2

1

1



5,31

5,31



1



1

 21.2 00625 ,0

 21.2 00625 ,0

00625  2

 2

1   

1   

  2  

5,31   

  

  





5,52



1

1

,0

0125

01875

r 3

r 4



 K



3

1

1



5,52

5,52



1



1

 21.2 00625 ,0

 21.2 00625 ,0

,0  2

 2

1   

1   

  2  

5,52   

  

  

- Phần tử 4 có đối lưu,nên ma trận độ cứng là :

 k 2 l r i



 r i

j

(2.264)

 K



0hr

r 4

r 5



 .2

 K



4

hr 5

1  1

00 10

 21.2 00625 ,0

 2

  

  

1

00



,0

,0

025

 ,0.2

025

.

4000



 21.2 00625

,0

01875  2



1

10

1

  

  

   1   

5,73



5,73

0



0



5,73

5,73

100

1    1     0   2  

  2  

  

  

- Véc tơ lực của các phần tử 1, 2, 3 không có đối lưu, phần tử 4 có đối lưu

j

(2.265)

f



l

 

 r i 2 

r r

.2 q V 6

2 r i

j

  

  

 r 1  00   .2  1 1 10 2  2 k l    1         

+ Véc tơ phụ tải nhiệt của các phần tử - Véc tơ phụ tải của các phần tử 1, 2 và 3 không có thành phần đối lưu là

,00.2 

82,229

 .2

35300000



l



,0



  f

1

r 1 

 2

q  .2 V 6



64,459

6

2 r 1

r 2 r 2

  

  

 00625  ,0.20 

00625   00625 

  2  

  

,0.2

00625



,0

27,919

 .2

35300000



l



,0



  f

2

r 2 

 2

q  .2 V 6

00625



,0.2

09,

6

2 r 2

r 3 r 3

  

  

 00625  ,0 

0125   0125 

  2  1149 

  

73

,0.2

0125



,0

27,

 .2

35300000



l



,0



  f

3

q  .2 V 6

6

0125



,0.2

54,

2 r 3

r r  3 4 2 r  4

  

  

 00625  ,0 

01875   01875 

1608   2  1838 

  

- Véc tơ lực của phần tử 4 có thành phần đối lưu là

0

T

T

(2.266)





hT

ds



l



hT

  f

 N



4

  dNq V

a

 .2. r 5

a





r 4 

 2

 q .2 V 6

2 r 4

r 5 r 5



S

  

  

  1 

  

0 ,0.2 01875  ,0 025  .2 35300000  ,0 00625   .2 4000 025,0.20. 01875  ,0.2 025 6   1          ,0 

5,31

5,52

 5,10 



5,73

 5,52

 5,31 

0 5,31  5,52

5,10  5,31 0

0 0 5,52 

5,10  5,10 0 0

0 0 0 5,73



82,229   

27,919 1608 27, 2298 18,

0

0

0

5,73



4528

00,

5,73



100

       

T  1  T  2 T  3  T 4   T  5

       

  64,459  09, 1149   1838 54,   

       

402,1146 380,2270





(2.267)





0 0 5,52  0,126

5,10  0,42  5,31 0

0 5,31 84 5,52





329,1562 245,9926 130,3081

5,173

5,73



0

0

0

4528

00,

0   0   0  5,73   

T  1  T  2 T  3  T 4   T  5

       

       

82,229   9, 1378  2757 81,   4136 72,   

T  1  T  2 T  3  T 4   T  5

       

       

       

0 2298 18, 2298 18,    2528 00, 2000 0, 4528 00,   2        2        2     

+ Lắp ghép các phần tử         Tức là 5,10    5,10   0  0    So sánh với nghiệm giải tích trong bảng 2.7 như sau Bảng 2.6. So sánh kết quả tính nhiệt độ

PTHH (0C) 402,1146 380,2270 329,1562 245,9926 130,3081

Chính xác (0C) 392,26 376,54 327,29 245,22 130,31

T1 T2 T3 T4 T5

74

2.12. Dẫn nhiệt qua cánh tiết diện thay đổi Khảo sát một phần tử cánh điển hình, tại vị trí i và j có bề dày

jd diện tích

jA , chu vi

id và

iA và

iP và

jP như trên hình 2.23.

Hình 2.23. Cánh mỏng dần và các vị trí i j

Từ hình vẽ chúng ta có diện tích tiết diện cánh và chu vi tại i,j là :

; và

(2.268)

b (2

bd

bd

(2

A













d

)

);

;

Pdb j i

A i

P i

j

j

j

i

Vì A thay đổi theo bậc nhất theo x



A

A i

j



1(



)



A

 AA i

 Ax i

j

L

x L

x L

nên có thể biểu thị theo hàm nội suy

(2.269)





ANANA i

i

j

j

L là chiều dài của phần tử. Bằng cách tương tự chu vi cũng biến đổi được thành

(2.270)





j

T

(2.271)





PNPNP i j i 1. Ma trận độ cứng Từ định nghĩa  K



     dBDB

   T dSNNh









S

lx 

lx 

và

ở đây  và S là thể tích và diện tích bao quanh miền khảo sát, nên:

dS











d  

Adx 0

x



S

Pdx 0

x



- Tính số hạng đầu:

     dBDB T



l



A

1



A i

j

T



x

dx

     dBDB

A i

 



0



1

1

l

k 2 l

  

1   

  

  

75

2



A



A

1



1



A i

j

A i

j





lA i



1

1

l

l 2

k l



1

1

2

k 2 l

  

1   

  

1   

  

  

  

  

Vậy



A

A i

j

(2.272)

     dBDB T





1  1

k l

2

  

1    1 

   

   

- Tính số hạng sau:

   T dANNh



A

l

l

i

j



h



N

h

   T dANNh



dxPNPN  

 N

 Pdx

i

j

i

i

j

j

A







N N

j

2 N i NN i

j

NN i 2 N j

0

0

  

  

   

   

l





2 i

i

j

j

j



h



 

 i 







 PNPNN i j PNPNNN i j

i

i

j

j

PNPNNN i i j PNPNN i j

2 j

i

j

0

   

   dx   

l

(

)

2 i

3 i

i

j

j

i



h



0

(

)

( )  PNNPN j 2 2 PNNPNN  i j i j

j

i

2 2 PNNPNN  i j i j ( PNNPN )  j

3 j

2 j

i

i

   

 dx   

dlNN



áp dụng công thức tích phân (2.113):

, với hai số hạng

a i

b j

l



ba !! ba 

(

)!1





l

l

;

3 dxN i



6 24

l

(2.273)

dxNN



l



l

2 i

j



2 24

!3!0    !130 !2!1    !121

l

sẽ được



P

P i

j



P

P i

j

l



P



P

N

12

P i

j

P i

i

  

(2.274)

   



h

N



hl

 N

 Pdx

i

j



1 12 P





P



3

j P

N

hl 12

j

P i

P i

j

P i

j

j

0

3   

  

  

  

P



j

P i

12

    1 4

1 12

  

  

1 4   

  

         

      

- Vậy ma trận độ cứng:

AA 



1





i

j

(2.275)



 K



1

hl 12

2

k l



1

PP  j i 3 P P  i j

1   

  

3 PP   i j   PP  j i

  

2. Véc tơ phụ tải Từ định nghĩa

76

T

T

T

(2.276)

 



 Nq





 Nq V

 NhT a







l

A

A

qV là mật độ nguồn thể tích, q mật độ dòng nhiệt bề mặt, h hệ số tỏa nhiệt bệ mặt, Ta nhiệt độ môi trường bao quanh; Adx = d , dA = Pdx. Xác định từng số hạng như sau. - Số hạng nguồn trong :

2

f  Adx  dA  dA

i

T

i

i

j

i



 Nq V

j

j

i

i

j 2







j

j

j

i

j

l

l

l

N ANNAN  Adx  ( ANAN  ) dx  dx q V q V N ANANN  i           

2

AA 

j

(2.277)





lq V

j 6 A

i 

2

A

lq V 6

A i

j

  

  

j 3

A i 3 A i 6

     

     

- Số hạng dòng nhiệt bức xạ :

2

   A

i

i

T

i

i

i

 Nq



j

j

i

i

j 2









j

j

j

j

i

A

l

l

l



2

(2.278)





ql

2

ql 6

P i

PP  i j  P j

  

  



P i 3 P i 6

P j 6 P j 3

     

     

- Số hạng đối lưu tại bề mặt xung quanh A:

T

T

(2.279)

dA

Pdx







 NhT a

 NhT a

 



lhT a 6

2 PP  i j P 2  P j

i

l

A

  

  

- Nếu mặt cuối cánh có diện tích An, , phần tử cuối sẽ có toả nhiệt biểu thị bởi

0

T

(2.280)





dA 0

AhT n a



 NhT a 0

A

  1 

  

Vậy véc tơ phụ tải nhiệt của mỗi phân tố sẽ là :

j

(2.281)



  f

AhT n a

 2

A A

PP 2  i j  P P 2

lq V 6

ql 6

lhT a 6

A 2 i A  i

j

2 P i

PP  i j  P 2 j

j

i

0   1 

  

  

  

  

  

  

  

Số hạng cuối chỉ có với phần tử cuối cùng.

77

N N dA  q Pdx  ( dxPNPN  )  q q dx N N PNNPN  j PNPNN  i j                    

3. Phương trình đặc trưng của phần tử Phương trình đặc trưng của phần tử ij đối với cánh có thiết diện thay đổi là



A





A i

j



1  1

2

k l

hl 12

j

PP  j i  P P 3 i

j

j

  

1    1 

PP 3   i j   PP  i

T  i  T 

  

   

(2.282)

2



A

2

j

j

j



AhT n a

A i 

2

A

PP  i 

P

2

   PP  2 i 

P

2

0 1

lq V 6

ql 6

    lhT a 6

A i

j

P i

j

P j

i

  

  

  

  

  

  

  

  

  T1 T2 T3   

10 mm 10 mm

Thí dụ 2.11. Khảo sát cánh phẳng bề dày nhỏ dần từ gốc dày 2 mm đến đỉnh dày 1 mm, như hình 2.19. Đỉnh cũng mất nhiệt ra môi trường, hệ số tỏa nhiệt, h = 120W/m2 0C , nhiệt độ môi trường Ta = 250C. Xác định phân bố nhiệt độ nếu gốc giữ nhiệt độ 1000C. Chiều dài tổng của cánh là L = 20 mm, chiều rộng cánh b = 3 mm. Hệ số dẫn nhiệt của vật liệu bằng 200W/m0C. Chia miền khảo sát thành hai phần tử chiều dài 10 mm như trên hình 2.24. Hình 2.24. Rời rạc phần tử hữu hạn Từ số liệu có

A1 = bd1 = 0,003.0,002 = 6,0.10-6 ; P1 = 2(b+d1) = 2(0,003+0,002) = 0,01 A2 = bd2 = 0,003.0,0015 = 4,5.10-6; P2 = 2(b+d2) = 2(0,003+0,0015) = 0,009 A3 = bd3 = 0,003.0,001 = 3,0.10-6 ; P3 = 2(b+d3) = 2(0,003+0,001) = 0,008 A4 = bd3 = 3,0.10-6

+ Ma trận độ cứng các phần tử - Phần tử 1:

1









 K



1

AA  1 2 2

k l



2

1

hl 12

PP  2 1 P P 3  1 2

  

1   

PP 3   1 2   PP  2 1

  

6 



 105,46 2

1   ,0 009 01,0  009,0  120   2 1 01,0  009,0 01,0  009,0.3 200 01,0 01,0. 12    1    01,0.3      

- Phần tử 2:

78

105,0  ,0 0039 0,019 0,1089 -   105,0 105,0 0,019 0,0037 - 0,1031 0,1086    105,0             0,1031   

1



0,0785

-





A 2

A 3



 K



2

 



2

1

-

0,0733

0,0783

 2

k l

hl 12

P 3 P 3

P 2 P 2

P  3 3 P  3

  

1   

P 3  2  P  2

  

  

0,0733   

+ Véc tơ phụ tải - Phần tử 1:

01,0.2



,0



  f

1

01,0



,0.2

lhT a 6

01.0.25.120 6

PP 2  2 1 P P 2  2

1

  

  

  

009   009 

  

0,1450   0,1400 

- Phần tử 2:

2

0



  f

AhT a 0

2

P

lhT a 6

2

PP  3 2 P 2  3

  1 

  

  

120

01,0.25.



6



10.3.01,0.25.120

008,0  008,0.2

0,1300 0,1340

6

  

   009,0.2   ,0 009  

0   1 

  

  

  

- Lắp ghép các phần tử

0,1089

-

0,1031

0

0,1450

-

(0,1086

0,0785)



0,1340)

0,1031 0

 0,0733

-

 0,1340

    

  - 0,0733   0,0783 

T  1  T  2  T  3

    

  (0,1400   

    

Phương trình đặc trưng tổng thể (chưa kể điều kiện biên)

,0

1089



,0

1031

0

1450





,0 1031 0

0,1871 - 0,0733

,0 270 134,0

    

  - 0,0733   0,0783 

T  1  T  2  T  3

    

,0     

    

- Áp đặt điều kiện biên T1 =1000C , phải thay đổi như sau Dòng 1 : T1 = 100; Dòng 2 : 0T1 + 0,1871T2 – 0,0733.T3 = 0,270 + 0,1031100 - Phương trình đặc trưng tổng thể. Sau khi thay thế, phương trình đặc trưng tổng thể trở thành

0

0

100.0000

 giải ra :





0 0,1871 -0 0,0733

89.2588 85.2703

1     

  - 0,0733   0,0783 

T  1  T  2  T  3

    

0,100     45,10     134,0  

T  1  T  2  T  3

    

    

    

79

2.13. Dân nhiệt ổn định hai chiều dùng phần tử tam giác Khảo sát bài toán dẫn nhiệt hai chiều của một phần tử tam giác 1 2 3, có diện tích tam giác là A, bề dày là  được thể hiện trên hình 2.25. Để bài toán mang tính tiêu biểu, nghĩa là có đủ các thành phần phụ tải nhiệt, ta cho mặt bên dưới ứng với cạnh 12 có dòng nhiệt q, mặt bên phải ứng với cạnh 23 có toả nhiệt với môi trường và trong tam giác có nguồn nhiệt phân bố đều qV. Mặt bên trái ứng với cạnh 31 được cách nhiệt.

Hình 2.25. Phần tử tam giác tiêu biểu.

(2.283)

TK

   f 

T

(2.284)

Phương trình đặc trưng cần xác định  Trong đó ma trận độ cứng phần tử  K

     dBDB



   T dSNNh



2







S

2

T

T

T

(2.285)

 

và véc tơ phụ tải  

 Nq





  dNq V

1

 NhT a

2







1

S

S

2

(2.286)

 f  dS  dS

   TNTNT 11 22 TN 33

 Trong các công thức trên : d = dA với A là diện tích tam giác,  là bề dày ; dS1 = dl12 với S12 và l12 là diện tích và chiều dài cạnh 12 tại mặt bên có dòng nhiệt q ; dS2 = dl23 với S23 và l23 là diện tích và chiều dài cạnh 23 tại mặt bên có dòng nhiệt đối lưu ; 1. Ma trận độ cứng Phân bố nhiệt độ trong phần tử tam giác được viết theo nhiệt đọ 3 nút là Các hàm nội suy

N







yc

1

 a 1

xb 1

1

1 A 2

80

(2.287)

N







 a

yc

2

2

xb 2

2

N







yc

3

xb 3

3

3

1 2 A 1  a 2 A

Với

(2.288)

 x y - x y  y - y  - x ; b 1 ; c 1

3 - x

2

2

3 ; c 2

a 1 a   

2 y x 13 y x 21

3 y 31 y 12

2 ; b 2 ; b 3

2 - y 1 ; c 3

x 3 - x 3    - x y 3 - y 2 y 1 x 1 - x 1 x 2 a 3

đều là các hằng số Đạo hàm hàm nội suy

2

(2.289)





  B

b c

b 3 c

1 A 2

2

3

b  1  c  1

  

N  1  x N  1  y

N  2  x N  2  y

N  3  x N  3  y

     

     

T

T

- Tính



     dBDB

  AdBDB    







A

Trường hợp tổng quát vật liệu không đẳng hướng thì hệ số dẫn nhiệt [D]:

(2.290)

D







0 k

k x 0

y

  

  

c 1

T



  

  BDB

k x 0

0 k

k x 0

0 k

b 2 c

b 3 c

1 A 2

1 A 2

y

y

2

3

b  1  c  1

  

  

  

  

  

c 2 c 3

b  1  b  2  b  3

    















2

2

1 A

4

1 A

4

bk x 1 ck y 1

bk x 2 ck y 2

bk x 3 ck y 3

  

  







c 1 c 2 c 3

ckc y 1 1 ckc y 2 1 ckc y 3 1

bkb x 2 1 bkb x 2 2 bkb x 2 3

ckc y 2 1 ckc y 2 2 ckc y 2 3

bkb x 3 1 bkb x 3 2 bkb x 3 3

ckc y 1 3 ckc y 2 3 ckc y 3 3

b  1  b  2  b  3

    

 bkb x 1 1  bkb  x 2 1  bkb  x 1 3

    

2

(2.291)

  

2 bk x 1 bbk x 21 bbk x 31

3

T

(2.292)



k



k

      dABDB

x

y



 4 A

A

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

2  b 1  bb  21 bb   31

    

    

    

    

    

81

    1 24 A    ck y 1 cck y 21 cck y 31 bbk x 21 2 bk x 2 bbk x 32 cck y 21 2 ck y 2 cck y 32 bbk x 31 bbk x 32 2 bk x 3 cck y 31 cck y 32 2 ck y          

3

- Tính



   T dSNNh



    T  dlNNh

2





2

S

2

N

N

1

2

3

(2.293)



N

N



N



  NN T



 N



1

2

3

2

2

3

N

3

2 1 NN 1 NN 1

3

NN 1 2 N 2 NN 2

3

NN 1 NN 2 2 N 3

    

    

    

    

Tại cạnh 23 có N1 = 0, còn N2 và N3 thay đổi giữa 0 và 1 như đã biết trong phần trước, bảng

Bảng 2.7

Nút 1 1

Nút 2 0

Nút 3 0

0

1

0

0

0

1

1N 2N 3N

0

0

3

T

Bởi vậy

(2.294)



h

   dSNNh

 

2

3

2





S

2

0 0

0 2 N 2 NN 2

3

NN 2 2 N 3

    

   dl .   

dlNN



áp dụng công thức tích phân (2.113):

với các số hạng trên

a i

b j

l



ba !! ba 

(

)!1





l



l

;

2 dxN 2

2 dxN 3





2 6

!2!0    !120

l

l

(2.295)



l



l

dxNN 3

2



1 6

!1!1    !111

l

Nên

(2.296)





2



 lh . 23 6

   T dSNNh 2

S

000     120     210  

Vậy ma trận độ cứng sẽ là

23

(2.297)



k



k



 K



e

x

y

 4 A

 lh . 6

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

2  b 1  bb  21 bb   31

    

    

    

000   120   210 

    

    

    

82

T

T

T

(2.298)



hT

dS

dS





f

Chỉ số e trong phương trình trên biểu thị phần tử đơn. 2. Véc tơ phụ tải nhiệt Công thức tính chung :  

  dNq

 Nq

 N





V

1

a

2







S

1

S

2



N

1

T

- Số hạng nguồn trong là

dA



  dNq V

q V

2







A

    

     

N N 3 Do nguồn trong phân bố đều trong tam giác nên



N



N



N



N



N



   001











 



1

1

11

21

31



N



N



N



N



N













 



2

2

12

22

32

N



N



N



N



N













 



3

3

13

23

33

1 3 1    010 3 1    100 3

1 3 1 3 1 3

1 3 1 3 1 3

Vậy:

T

(2.299)

dA



  dNq V



 

q V 3

 Aq V 3



A

1      1     1  

1     1     1  

N

1

2

T

- Số hạng dòng nhiệt tại cạnh 12 là



q



 Nq



dS 1

2





S

1

1

N N

3

    

   .dl   

Trên cạnh 12 có N3 = 0; còn N1 và N2 thay đổi giữa 0 và 1, nên áp dụng công thức tích phân



đối với N1 và N2 thì đều có

b dlLL j

a i

l



ba !! ba 

(

)!1





dlNN





dlN i

dlN j

1 i

0 j







!0!1  )!101(

1 2

l

l

l

Vậy:

T

(2.300)

dS



 Nq



1

 

 lq . 12 3

S

1

1     1     0  

83

N

1

3

T

- Số hạng toả nhiệt trên cạnh 23 là

dS



hT



 NhT a

2

a

2

2





S

2

N N

3

    

   .dl   

Cũng tương tự như trên, N1 = 0 tại 23, còn N2 và N3 tính theo tích phân số, nên có

T

(2.301)

dS





 NhT a



 Th a 2

S

0   1   1 

    

Véc tơ phụ tải nhiệt

12







  f

(2.302)

e

 GA 3

 lq . 2

 hT l . a 13 2

1     1     1  

1     1     0  

0     1     1  

3. Phương trình đặc trưng của phần tử tam giác

23

000

x

y

2  b 1  bb  21 bb   31

2  c 1  cc  21 cc   31

(2.303)

k  k  120   4 A  lh . 6 210 bb 21 2 b 2 bb 32 bb 31 bb 32 2 b 3 cc 21 2 c 2 cc 32 cc 31 cc 32 2 c 3                     T  1  T  2  T  3                         

12

Nếu nguồn nhiệt không phân bố đều trong miền , mà tập trung tại một điểm có toạ độ (x0, y0) gọi là “nguồn điểm” q* thì số hạng nguồn thứ nhất trong vế phải là

N

1



được thay bằng

N

, với [N] tính tại (x0, y0) (2.304)

2

AqV 3

Aq * 3

N

3

1     1     1  

    

    

x )0,0( y

 hT    Aq V 3  lq . 2  l . a 13 2 1     1     1   1   1   0       0   1   1      

Nguồn điểm q* có đơn vị (W/m) vì phân bố theo bề dày  của tam giác. Thí dụ 2.12. Hình vuông phẳng có bề dày bằng 1m, kích thước 10 cm như trên hình 2.22. Tại cạnh trên cùng ở đỉnh có nhiệt độ 5000C, ba cạnh còn lại đều có nhiệt độ là 1000C. Biết hệ số dẫn nhiệt của vật liệu là không đổi k = 10W/m 0C. Sử dụng các phần tử tam giác bậc nhất để xác định phân bố nhiệt độ trong hình. 1. Rời rạc các phần tử

84

7 8 9  

 

4 5 6    

Miền vuông được chia thành 8 phần tử tam giác bậc nhất có kích thước bằng nhau như trên hình 2.26.

1 2 3 Hình 2.26. Rời rạc các phần tử trong hình vuông

Từ (2.288) có nhận xét rằng, các đại lượng bi, ci trong ma trận là hiệu tọa độ của các nút của tam giác, các hiệu này không phụ thuộc vào vị trí tam giác trong miền. Nói cách khác khi hai tam giác bằng nhau và có cùng hướng thì [K] là như nhau. Như vậy sẽ có hai dạng ma trận độ cứng của phần tử theo sự định hướng của các tam giác. Các phần tử 1,3,5 và 7 cùng có ma trận độ cứng là [K]1 , các phần tử 2,4,6 và 8 cùng có ma trận độ cứng là [K]2. Nghĩa là trong bài toán trên chúng ta chỉ cần xác định hai ma trận [K]1 và [K]2 là đủ. Tuy nhiên ở đây chúng ta vẫn xác định tất cả các ma trận của từng phần tử nhằm trình bày và giải thích cách tính tổng quát, sau này có thể dùng cho các bài toán khác. Để dễ dàng thực hiện tính toán, chúng ta lập bảng thể hiện quan hệ giữa số của phần tử , số các nút, thứ tự của nút cũng như tọa độ của chúng. Lấy gốc tọa độ là nút 1, thì các nút khác sẽ có tọa độ như sau, bảng 2.8 Bảng 2.8. Tọa độ các nút

Nút Tọa độ (m)

x y

1 0 0

2 0,5 0

3 1 0

4 0 0,5

5 0,5 0,5

6 1 0,5

7 0 1

8 0,5 1

9 1 1

2

3

4

5

6

7

8

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

1 2 4 2 5 4 2 3 5 3 6 5 4 5 7 5 8 7 5 6 8 6 9 8

1 0,5

Bảng 2.9 . Số nút của mỗi phần tử, thứ tự nút (còn gọi là bậc tự do) và tọa độ các nút Phần tử 1 Thứ nút Nút số Tọa độ

x 0 0,5 0 0,5 0,5 0 0,5 1 0,5 1 1 0,5 0 0,5 0 0,5 0,5 0 0,5 1 0,5 1 y 0 0 0,5 0 0,5 0,5 0 0 0,5 0 0,5 0,5 0,5 0,5 1 0,5 1 1 0,5 0,5 1 0,5 1 1

Công thức chung (chỉ số theo thứ tự nút trong mỗi phần tử)

b1 = y2 – y3 c1 = x3 – x2

b2 = y3 – y1 c2 = x1 – x3

b3 = y1 – y2 c3 = x2 – x1

Tính các hệ số b1,b1, và c1, c2 trong mỗi phần tử Phần tử 1 b1 = y2 – y4 = 0 – 0,5 = - 0,5 b2 = y4 – y1 = 0,5 – 0 = 0,5 c1 = x4 – y2 = 0 – 0,5 = - 0,5 c2 = x1 – x4 = 0 – 0 = 0

b3 = y1 – y2 = 0 – 0 = 0 c3 = x2 – x1 = 0,5 – 0 = 0,5

Phần tử 3 b1 = y3 – y5 = 0 – 0,5 = - 0,5 b2 = y5 – y2 = 0,5 – 0 = 0,5 c1 = x5 – y3 = 0,5 – 1 = - 0,5 c2 = x2 – x5 = 0,5 – 0,5 = 0

b3 = y2 – y3 = 0 – 0 = 0 c3 = x3 – x2 = 1 - 0,5 = 0,5

85

Phần tử 5 b1 = y5 – y7 = 0,5 – 1 = - 0,5 b2 = y7 – y4 = 1 - 0,5 = 0,5 c1 = x7 – y5 = 0 – 0,5 = - 0,5 c2 = x4 – x7 = 0 – 0 = 0

b3 = y4 – y5 = 0,5 – 0,5 = 0 c3 = x5 – x4 = 0,5 – 0 = 0,5

Phần tử 7 b1 = y6 – y8 = 0,5 – 1 = - 0,5 b2 = y8 - y5 = 1 - 0,5 = 0,5 c1 = x8 – y6 = 0,5 – 1 = - 0,5 c2 = x5 - x8 = 0,5 - 0,5 = 0

b3 = y5 – y6 = 0,5 – 0,5 = 0 c3 = x6 – x5 = 1 - 0,5 = 0,5

Phần tử 2 b1 = y5 – y4 = 0,5 – 0,5 = 0 b2 = y4 – y2 = 0,5 – 0 = 0,5 c1 = x4 – y5 = 0 – 0,5 = - 0,5 c2 = x2 – x4 = 0,5 – 0 = 0,5

b3 = y2 – y5 = 0 – 0,5 = - 0,5 c3 = x5 – x2 = 0,5 – 0,5 = 0

Phần tử 4 b1 = y6 – y5 = 0,5 – 0,5 = 0 b2 = y5 – y3 = 0,5 – 0 = 0,5 c1 = x5 – x6 = 0 – 0,5 = - 0,5 c2 = x3 – x5 = 1 - 0,5 = 0,5

b3 = y3 – y6 = 0 – 0,5 = - 0,5 c3 = x3 – x6 = 1 - 1 = 0

Phần tử 6 b1 = y8 – y7 = 1 – 1 = 0

b2 = y7 – y5 = 1 - 0,5 = 0,5 c1 = x7 – y8 = 0 – 0,5 = - 0,5 c2 = x5 – x7 = 0,5 – 0 = 0,5

b3 = y5 – y8 = 0,5 - 1 = - 0,5 c3 = x8 – x5 = 0,5 – 0,5 = 0

Phần tử 8 b1 = y9 – y8 = 1 - 1 = 0

b2 = y8 – y6 = 1 - 0,5 = 0,5 c1 = x8 – y9 = 0,5 – 1 = - 0,5 c2 = x6 – x8 = 1 - 0,5 = 0,5

b3 = y6 – y9 = 0,5 - 1 = - 0,5 c3 = x9 – x6 = 1 – 1 = 0

Thấy rõ ràng rằng các phần tử 1, 3, 5, 7 có b1,b1, và c1, c2 như nhau và các phần tử 2, 4, 6, 8 cũng vậy, có b1,b2, và c1, c2 như nhau 2. Ma trận độ cứng các phần tử Tính các ma trận độ cứng [K]1 và [K]2 , với  =1, kx = ky có

(2.305)





  K e

k 4 A

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

2  b 1  bb  21 bb   31

    

    

    

    

    

Tính diện tích tam giác A = D/2

1

1

0

0

x 1

y 1

. Vậy A = 0,25/2 2.306)

D





25,0



2

1 1

y y

5,01 1 0

0 5,0

x 2 x 3

3

Tính [K]1 Các số hạng trong móc vuông

2

 







25,0



1  1 0

01  0 1 0 0

  5,0   5.0.5,0    0.5,0 

5.0.5,0   2 5,0  0.5,0



 0.5,0    0.5,0   2 0

2 b 1 bb 21 bb 31

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

    

    

    

    

    

     

86

2

5,0

0

1





 



 0.5,0



25,0

0

0

0

 5,0.5,0



01

1

     



0.5,0   2 0  5,0.0



 5,0.5,0    5,0.0   2 5,0

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

    

    

    

1     

    

     

Vậy

2



1



2



1









5



 K



1

k 4 A

.4

 

1 1

1 0

0 1

 

1 1

1 0

0 1

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

2  b 1  bb  21  bb  31

    

    

    

    

1     

    

1     

    

    

25,0.10 25,0 2

Tính [K]2 Phần tử 2 b1 = y5 – y4 = 0,5 – 0,5 = 0 b2 = y4 – y2 = 0,5 – 0 = 0,5 c1 = x4 – y5 = 0 – 0,5 = - 0,5 c2 = x2 – x4 = 0,5 – 0 = 0,5

b3 = y2 – y5 = 0 – 0,5 = - 0,5 c3 = x5 – x2 = 0,5 – 0,5 = 0

2







0 0 0

0 1  1

 5,0.0   2 5,0   5,0.5,0 

  5,0.0   5,0.5,0   2 5,0 



bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2  b 1  bb  21 bb   31

    

  25,0   

0    1   1 

   0    5,0.0     5,0.0  

     

2









25,0



1  1 0

01  0 1 0 0

  5,0   5,0.5,0     0.5,0

5,0.5,0    2 5,0  0.5,0



 0.5,0    0.5,0   2 0

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

    

    

    

    

    

     

Vậy



25,0





5



 2K

.4

0 0 0

0 1  1

1  1 0

01  0 1 0 0

1  1 0

1  2  1

    

0    1   1 

    

    

    

0    1   1 

    

    

10 25,0 2

với e = 1  8, nghĩa là có thể bỏ hệ số 5 trong

   e TK

Do dẫn nhiệt ổn định, nên các phần tử đều có eK 

2



1



1



1

0

(2.307)



;



;

 K



 K



e

 7,5,3,1

e



8,6,4,2

 

1 1

1 0

0 1

 1 0

2  1

    

1     

    

   1   1 

0e

3. Lắp ghép các phần tử Quá trình lắp ghép các ma trận đặc trưng của từng phần tử thành ma trận tổng thể của cả hệ là thủ tục hết sức quan trọng, đặc biệt các hệ thống lớn gồm hàng trăm, đến hàng ngàn phần tử. Để hiểu rõ quá trình lắp ghép này, dưới đây chúng ta sẽ trình bày và giải thích một cách tỉ mỉ các bước tiến hành: 3.1. Các bước cơ bản , mục đích của mỗi bước

87

a. Tách ma trận của từng phần tử thành các phương trình đại số theo nhiệt độ tại các nút, đánh số phương trình : Phần tử 1: gồm các nút 1,2 và 4

Phần tử 2: gồm các nút 2,5 và 4

Nút 1: (2T1 – T2 – T4) = 0 (1) Nút 2: (-T1 + T2 + 0T4) = 0 (2) Nút 4: (-T1 + 0T2 + T4) = 0 (3) Phần tử 3: gồm các nút 2,3 và 5

Nút 2: ( T2 – T5 + 0T4) = 0 (4) Nút 5: (-T2 + 2T5 - T4) = 0 (5) Nút 4: ( 0T2 – T5 + T4) = 0 (6) Phần tử 4: gồm các nút 3,6 và 5

Nút 2: (2T2 – T3 – T5) = 0 (7) Nút 3: (-T2 + T3 + 0T5) = 0 (8) Nút 5: (-T2 + 0T3 + T5) = 0 (9) Phần tử 5: gồm các nút 4,5 và 7

Nút 3: (T3 – T6 + 0T5) = 0 (10) Nút 6: (-T3 + 2T6 – T5) = 0 (11) Nút 5: (0T3 – T6 + T5) = 0 (12) Phần tử 6: gồm các nút 5,8 và 7

Nút 4: (2T4 – T5 – T7) = 0 (13) Nút 5: (-T4 + T5 + 0T7) = 0 (14) Nút 7: (-T4 + 0T5 + T7) = 0 (15) Phần tử 7: gồm các nút 5,6 và 8

Nút 5: (T5 – T8 + 0T7) = 0 (16) Nút 8: (-T5 + 2T8 – T7) = 0 (17) Nút 7: (0T5 – T8 + T7) = 0 (18) Phần tử 8: gồm các nút 6,9 và 8 Nút 6: (T6 – T9 + 0T8) = 0 (22) Nút 9: (-T6 + 2T9 – T8) = 0 (23) Nút 8: (0T6 – T9 + T8) = 0 (24)

Nút 5: (2T5 – T6 – T8) = 0 (19) Nút 6: (-T5 + T6 + 0T8) = 0 (20) Nút 8: (-T5 + 0T6 + T8) = 0 (21) b. Lập bảng Thống kê các nút và số của phương trình trong mỗi phần tử thành bảng để dễ dàng nhận biết , rồi chuyển sang bảng số phương trình theo các nút. Bảng 2.10. Bảng thể hiện số nút, số của phương trình theo mỗi phần tử

2

3

5

6

8 6 8 6 9 8

7 4 1 1 2 4 2 5 4 2 3 5 3 6 5 4 5 7 5 8 7 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Phần tử Nút số Phương trình số

Bảng 2.11. Bảng thể hiện số của phương trình tại mỗi nút

7

9

1 1

2 2,4,7

3 8,10

4 3,6,13,

6 11,20,22 15,18

8 17,21,24 23

Nút số Phương trình số

5 5,9,12, 14,16,19

Số nút của toàn miền là 9, nghĩa là cần có 9 phương trình biểu thị nhiệt độ của hệ tại 9 nút. Muốn vậy phải lắp ghép các phần tử rời rạc biểu thị bởi 24 phương trình trên lại với nhau để tạo thành 9 phương trình. Do có những nút là chung của các phần tử xung quanh nên phải thỏa mãn mối quan hệ nhiệt độ với các phần tử xung quanh. Bởi vậy nguyên tắc lắp ghép là cộng toàn bộ các phương trình có cùng số nút lại để mỗi nút chỉ có 1 phương trình duy nhất. Muốn vậy phải lập các bảng sau c. Cộng các phương trình có trong một nút Nút 1: chỉ có phương trình 1: 2T1 – T2 – T4 = 0 (25) Nút 2: có các phương trình 2, 4 và 7: ( -T1 + T2 + 0T4 ) = 0 (2)

( T2 – T5 + 0T4 ) = 0 (4)  -T1 + 4T2 - T3 – 2T5 = 0 (26) ( 2T2 – T3 – T5 ) = 0 (7)

88

Nút 3: có các phương trình 8 và 10

(-T2 + T3 + 0T5 ) = 0 (8)

(T3 – T6 + 0T5 ) = 0 (10)  -T2 + 2T3 - T6 = 0 (27) Nút 4: có các phương trình 3,6 và 13 ( -T1 + 0T2 + T4 ) = 0 (3) ( 0T2 – T5 + T4 ) = 0 (6)  -T1 + 4T4 – 2T5 – T7 = 0 (28) ( 2T4 – T5 – T7 ) = 0 (13) Nút 5: có các phương trình 5,9,12,14,16 và 19 (-T2 + 2T5 - T4 ) = 0 (5) ( -T2 + 0T3 + T5 ) = 0 (9) ( 0T3 – T6 + T5 ) = 0 (12) ( -T4 + T5 + 0T7) = 0 (14)  - 2T2 – 2T4 + 8T5 – 2T6 – 2T8 = 0 (29) ( T5 – T8 + 0T7 ) = 0 (16) ( 2T5 – T6 – T8 ) = 0 (19)

Nút 6: có các phương trình 11, 20 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

2  1

1  4

0  1

1  0

0  2

0 0

0 0

0 0

0 0

0



1

2

0

0



1

0

0

0



1

0

0

4



2

0



1

0

0





2

0

0



2

8



2

0



2

0

0

0



1

0



2

4

0

0



0 0

0 0

0 0

1  0

0  2

0 0

2  1

1  4

0

0

0

0

0



1

0



1

2

             

        1   0   1   

0   0   0  0   0   0  0   0   0 

              

T  1  T  2  T 3  T  4  T  5  T 6  T  7  T  8  T  9

              

89

1

2

0

1









0 1

1 1

1 0

 

 1 0

  4   5  6

  1   2   3

    

    

    

    

    

1     

0   0   0 

T  1  T  2  T  4

T  2  T  5  T  4

   1   1 

0

1

1

2









1 1

0 1

1 0

 

 1 0

  7   8   9

  10   11  12

    

    

    

    

    

1     

0   0   0 

T  3  T  6  T  5

T  2  T  3  T  5

   1   1 

2

1





1

0

1





 1 0

 

0 1

1 0

1 1

  16   17 18 

  13   14   15

   1   1 

T  5  T  8  T  7

T  4  T  5  T  7

0   0   0 

1     

    

    

    

    

    





1

1

2

0

1





 1 0

2  1

 

0 1

1 0

1 1

  19   20   21

  22   23 24 

   1   1 

T  6  T  9  T  8

T  5  T  6  T  8

0   0   0 

1     

    

    

    

    

    

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:  0 1     0 2     0  1   Phần tử 3: Phần tử 4:  0 1     0 2     0  1   Phần tử 5: Phần tử 6:  0     0 2     0  1   Phần tử 7: Phần tử 8:  0     0     0   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ử

3

2

4

8

5

7

6

1 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

Thứ tự nút trong phần tử 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:

90

Để 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ố

Hệ số của nhiệt độ có mặt trong phương trình T5

T9

T2

T1

T3

T4

T6

T7

T8

1 2

3

4

5

6

7

8

Số của phương trình tại mỗi nút 1 2 4 7 8 10 3 6 13 5 9 12 14 16 19 11 20 22 15 18 17 21 24 23

9

2 -1 -1

-1 +1 +1 +2 -1 -1 -1

-1 +1 +1 -1

-1 +1 +1 +2 -1 -1 -1

-1 -1 -1 -1 +2 +1 +1 +1 +1 +2 -1 -1 -1 -1

-1 -1 -1 +2 +1 +1 -1

-1 +1 +1 -1

-1 -1 -1 +2 +1 +1 -1

-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 độ.

91

2  1

1  4

0  1

1  0

0  2

0 0

0 0

0 0

0 0

0



1

2

0

0



1

0

0

0



1

0

0

4



2

0



1

0

0



0



2

0



2

8



2

0



2

0

(2.308)

0

0



1

0



2

4

0

0



0 0

0 0

0 0

1  0

0  2

0 0

2  1

1  4

0

0

0

0

0



1

0



1

2

             

        1   0   1   

0   0   0  0   0   0  0   0   0 

              

T  1  T  2  T 3  T  4  T  5  T 6  T  7  T  8  T  9

              

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:

sinh

y

n

 1



 )1(



1

(2.310)

yxT ,(

)





)

sin

x



T ( 2

T 1

T 1



2 

n

 n b

n

 1

  

  

sinh

H

 n b  n b

     

     

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

0

(2.311)













100



100



500



200

C



 100



T 5

T 4

T 6

T 8

 T 2

1 4

1 4

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)

92

3(0;50) 2(50;50)  

4(40,40)

1(50;0)

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

; i = 1,2,3 (2.312)

N







 a

yc

i

i

xb i

i

1 A 2

Toạ độ các nút:

Nút Tọa độ (m)

x y

1 50 0

2 50 50

3 0 50

Xác định các hệ số ai, bi , ci :

50.0

2

 

 

50.50   0

 50.50

 

y y

 

50 50

 

50  ;0 c  1 c  0 ;50

0   50

50   0

50 50

a 1 a

2

;  2500 b 1 b  2500 ; 2 

3 

y 3 y 1  0

x  3 x  1 x

2 c

yx 3 2 yx 13 yx 1

2

yx 2 3 yx 31 yx 12

2

x  2 x  3 x 3

2

3

3

Tính D = 2A:

2

A



2500

50 50 0

0 50 50

1   1   1 

    





N







.0

x



50

y





 2500



 50

y



1

yc 1

xb 1

N









2500



50

x



50

y



50

 x

 

 

y

1 50 

 a



yc 2

xb 2

2

2

1 50

N











50

x



50

y



50(



x

)



 a



yc 3

xb 3

3

3

Tính các hàm nội suy N1, N2 và N3 1  a 1 2500 1 2500 1  2500 2500

1 A 2 1 2 A 1 A 2

1 50

93

   50.50  0 2500 ; y 50  ;50    50  50  0 b 3 y 1 a

Tại điểm 4 có x = 40; y = 40

(2.313)

N



;

N



;

N



1

2

3

1 5

3 5

1 5

Tính nhiệt độ tại điểm 4:









100



200



100



20



40.3



20



160

C0

TN T 11

4

TN 2 2

TN 33

1 5

3 5

1 5

Mật độ dòng nhiệt:



k



k



k













     TN

qx

 TNTN 22

11

TN 33

 Tb 11

Tb 22

33 Tb

k 2 A

 x 

 x 

2





200.50



.50

100

cmW

20

/

 100.0





k



k



k















     TN

q y

TNTNTN 33 22

11

 Tc 11

Tc 22

33 Tc

k 2 A

 x 

 y 

2



.50

100



200.50



100.0

cmW

20

/

 



 T  x 10 2500  T  y 10 2500

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 = 300C, 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 = 1,2 W/cm2 Ta = 300C

qV =1,2W/cm2 3 4

q = 2W/cm2



T=1000C

5 cm

q* (5W/cm)

 (1,1)

1 2

Cách nhiệt

Hình 2.28

5 cm

Bảng 2.14

Nút

x y

Tọa độ (m) Phần tử 1 Phần tử 2

1 0 0 1

2 5 0 2 1

3 0 5 3 3

4 5 5 2

94

T

(2.314)

     dBDB

   T dSNNh



Công thức chung  K







T

T

T

(2.315)

 









 S  Nq

  dNq V

 . NTh 









S

S

f   dS  dS

Phần tử 1: Tam giác 1 2 3 -Tính A, các hệ số ai,bi, ci :

2

A



25





2

2

x 1 x x

y 1 y y

    



5

3 yx 2 3 

a 1 

     3 yx 2 3 

001   051   501  0.05.5      ;05.00.0

 y

y 2 

 50 y 3 c  ;505

 x 

50  x 2  000

2

a

 y 1 y

;5 c 1 x  1 x

2 c

1   1   1   yx 13 yx 1

2

 yx 31 yx 12

;25 b 1 b  2 b 3

3 y 1

2

2

3 x 3 x 3

3

a1 = 25,0 b1=-5,0 b2 = 5,0 a2 = 0,0 b3 = 0,0 a3 = 0,0

c1 = -5,0 c2 = 0,0 c3 = 5,0

+ Tính [K]1 , do không có đối lưu nên



k



k

 K



1

x

y

1 A 4

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

2  b 1  bb  21 bb   31

    

    

    

    

    

25



25

0

25

0



2



1





25

0



0

0

0



1

1

0

25

2 25.2

0

0

25

0



1

0

1

0

     

    

     

25     25 

      

1     

    

    

T

T



   ;00.50.0      ;000    505 a 3

  dNq

 Nq



+ Tính [f]1 ,do có nguồn điểm và bức xạ cạnh bên:   f 1







S

N

i

- Số hạng nguồn điểm trong

, tại điểm (1;1) có

N

j

N

k

   q    

    

)1;1(

N











5

x



5

y







 25



1

 a 1

xb 1

yc 1

1 A 2

1 A 2

15 25

3 5

95

  dS

N









 50(

x



y )0





 a



2

2

xb 2

yc 2

N









x



5

y





 a



  00



3

3

xb 3

yc 3

1 2 A 1 A 2

5 25 5 25

1 5 1 5

1 2 A 1 A 2

Vậy nguồn điểm trong  =1

N

i

N



5



j

1 5

N

k

  q  *   

    

3     1     1  

3   1   1 

    

)1;1(

3

- Số hạng bức xạ cạnh bên:

1

1

T



dS



q

 Nq



2

3





1 2

q  

S

N N N

3

    

   dl   

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

dlNN



đối với N1 và N3 thì đều có

a i

b j

l



ba !!  ba

(

)!1





dlNN





dlN 3

1 i

0 j

dlN 1





!0!1  )!101(

1 2

l

l

 l Vậy:

T

dS







 Nq



 

31ql 2

5.2 2

S

1     0     1  

1   0   1 

    

5   0   5 

    

Véc tơ tải



T

T





dS







  f

 dNq .*



 Nq



1







S

3     1     1  

5   0   5 

    

    

2   1    4 

- Phương trình đặc trưng của phần tử 1

2



1





2

T 1

(2.316)



 

1 1

1 0

0 1

1  4

    

1     

  T  2  T  3

    

    

    

Phần tử 2:

96

T

     dBDB

   T dSNNh



+ Tính [K]2 : 

 K 2





S



Phần tử 2 Nút Tọa độ (m)

x y

1 2 5 0

2 4 5 5

3 3 0 5

Diện tích





25



25



25



25

2



A

2

2

x 1 x x

y 1 y y

3

3

    

051   551   501 

    



3

2

1   1   1  Các hệ số a 1 a

 

5.05.5      5.50.0

 2 

y y

 

x 

 

2

5

 ;25 b y 1 b  ;25 2 

;055 c  x 3 1 y c  05 ;5 1  50

2 

 50 5 x x   05 1 3  

3 y

yx 3 2 yx 13 yx 1

2

 yx 2 3 yx  31 yx 12

2

2

Tính [K]2

23





k



000 120

 K



2

x

y

1 A 4

 lh . 6

210

2 b 1 bb 21 bb 31

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

2 c 1 cc 21 cc 31

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

    

    

    

    

    

    

  k   

    

- Số hạng đầu

0

25



25

0

0

1



1

0











k 4 A

2 25.2

25  25

25 0

25 0

0 0

1  0

2  1

bb 21 2 b 2 bb 32

bb 31 bb 32 2 b 3

cc 21 2 c 2 cc 32

cc 31 cc 32 2 c 3

    

2  c 1  cc  21 cc   31

    

0   0   0 

  25    25 

    

    

    

  1    1 

 2  b 1   bb   21  bb    31

    

    

    

000

000





. 23lh 6

5.1.2,1 6

120 210

- Số hạng sau     

  120   210 

    

    

000     120     210  

Vậy ma trận độ cứng

97

    0.55.5  ;25  ;5 x 055 a 3 b 3 y 1 c 3 x 1









 2K

1  1 0

1  2  1

1  1 0

    

0    1   1 

    

000   120   210 

    

01    0 4   0 3 

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



12

23







  f

e

Aq V 3

 lq . 2

 l . hT a 2

1     1     1  

1   1   0 

    

  1   1 

0     

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à

Ta = 30



43

3

4





  f

2

Aq V 3

 l . hT a 2

qV

1     1     1  

  1   1 

0     

T = 100

2

0

5











  f

2

1.5.30.2,1 2

1.25.2,1 2.3

Thay số với qV =1,2; A=25/2; =1; h=1,2; Ta=30; l34=5     

5     5     5  

  90   90 

1     1     1  

0   1   1 

    

    95     95  

- Phương trình đặc trưng của phần tử 2

1



01

5

(2.317)



1  0

4 0

    

  0   3 

T  2  T  4  T  3

    

    95     95  

Hình 2.31. Phầ n tử 2

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ử

1



01

5

)4(

2



1





2

)1(

; ghép với





 1 0

4 0

)5( )6(

1 0

0 1

1  4

)2( )3(

    

  0   3 

T  2  T  4  T  3

    

    95     95  

    

1     

T  1  T  2  T  3

    

    

    

 1  1 - Bảng lắp ghép Bảng 2.15

98

Phụ tải

Nút số 1 2

3

4

Phươn g trình 1 2 4 3 6 5

T1 2 -1 -1

Hệ số nhiệt độ T3 T2 -1 -1 1 1 1 3 -1

-2 1 5 -4 95 95

T4 -1 4

- Lập phương trình ma trận đặc trưng tổng

2



1



1

0



1

2

0







1

0

4

0

91

0



1

0

4

95

     

  1    

T  1  T  2  T  3  T  4

      

  6     

2       

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

 giải ra

(2.318)





2 0 0 1  01

0

0

1

0

Nên dạng ma trận là 01    0 0   4 0  

     

T  1  T  2  T  3  T  4

      

98   100   91   100 

      

T  1  T  2  T  3  T  4

      

69   100   40   100 

      

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

99

(2.319)

Lấy gốc tại điểm 3 (k), các hàm nội suy sẽ là

N



1



1

y 2 a

x 2 b

  

  

1



N



2

(2.320)

  1     y 2 a

  

N



3

N



1



4

x 2 b

x   2 b  xy 4 ab y   2 a 

  

Ma trận đạo hàm của hàm nội suy là



2(

a



y

)

2(

a



y

)



y

y

(2.321)





  B

1 ab

4



b 2(



y

)



x

b 2(



y

)

x

  

  

N  3  x N  3  y

N  2  x N  2  y

N  4  x N  4  y

     

     

T

    TNTNT 11 22 TN 44 TN 33

N  1  x N  1  y Ma trận dộ cứng sẽ là  K



   dVBDB

 

   T dSNNh









S

, và [B]T[D][B] là

D



trong đó 



k x 0

0 k

y

  

  

 

 2( a  y )  b 2(  y )

T

  

  BDB

x 0

y

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à

2

b

2

a

2

a

2

b

2

a

2

b

k

2

2

(2.322)

)

y

a



2(

dxdy



2(

 xb

)

dxdy



dxdy

2

2

0

0

0

0

0

0



 

 

xy 4 ab

k x 2 ba

16

y 2 ba

16

 Sau khi tích phân sẽ được

100

k 0 2( a  y )  x  2( a  y ) 2( a  y )  y y  k 1 ab 4 y x 1 ab 4  2( b  y )  x b 2(  y ) x              y 2( b  y )            

2  2  1 1 2  2  1 1 0000

(2.323)

 K



2 1  1 2 1  1  2 0000  2    1 2  2 1 2  2  1  1 ak x 6 b bk y a 6 hl 12

Véc tơ tải là

1

2

b

2

a

T

2

(2.324)

f



dA



G



 

 NG



0

0







GA 4

3



N N N N

4

      

   dxdy    

1     1     1     1  

1  1  2 2  1  2 2 1                             2400  4200       

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)

2  2  1 1 2  2  1 1 0000

 K



2 1  1 2 1  1  2 0000  2    1 2  2 1 2  2  1  1 ak x 6 b bk y a 6 hl 12

Thay số vào được

101

1  1  2 2  1  2 2 1                             2400  4200       

2  2  1 1 2  2  1 1

 K



 2 2 1  1  2 2 1  1     1 1 2  2  1 1 2  2 5 15 5 15

1  1  2 2 1  1  2 2                         0000   0000   2400  4200       

sau khi biến đổi được

8  2  4  2

(2.325)

 K



 2 8  2  4   4  2 20 4 1 6

1

hT

2

31







f



q

 

 lq 14 2

a l . 2

 6 4

3

N N N N

4

1   0   0   1 

      

0   0   1   1 

      

 2  4 4 20            

Véc tơ phụ tải 1     1     1     1  

   *     

      

thay các giá trị vào được

7,5

(2.326)



  f

3,93

  3,8   7,97   

      

Phương trình ma trận đặc trưng cuả phần tử

2

(2.327)



8  2 4 

2  8  2

4  2  20

2  4  4

1 6

3



2



4

4

20

3,93

4

     

     

T  1  T   T   T 

      

7,5   3,8   7,97   

      

Áp điều kiện biên, với T2 = T3 = 1000C, sẽ được phương trình đặc trưng và giải ra nghiệm

8

0

0

- 2

, 2 634

T 1

88 , 4846

0

1

0

0

100 , 0

(2.328)

 giải ra   T



0

0

1 0

100 , 0

100 , 0000  100 , 0000

-

2

0

0

, 8 559

4

     20 

  T  2  T  3  T 

      

      

      

     

102

36 , 8385              

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

103