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