›
›

Bài giảng Truyền nhiệt cho các lớp cao học cơ khí - Phương pháp sai phân hữu hạn và phần tử hữu hạn trong truyền nhiệt - Trịnh Văn Quang

pdf 103 trang cucquyet12 6510
Download
Bạn đang xem 20 trang mẫu của tài liệu "Bài giảng Truyền nhiệt cho các lớp cao học cơ khí - Phương pháp sai phân hữu hạn và phần tử hữu hạn trong truyền nhiệt - Trịnh Văn Quang", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfbai_giang_truyen_nhiet_cho_cac_lop_cao_hoc_co_khi_phuong_pha.pdf

Nội dung text: Bài giảng Truyền nhiệt cho các lớp cao học cơ khí - Phương pháp sai phân hữu hạn và phần tử hữu hạn trong truyền nhiệt - Trịnh Văn Quang

  1. 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
  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 2
  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ạ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
  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 :  2T  2T 0 (2.1) x 2 y 2 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 T T i, j i 1, j x x x T T T T i , j i, j 1 y y y Hình 2.1. Mạng các điểm nút  2T ( T ) (T T ) (T T ) i 1 j i, j i, j i 1, j (2.2) x 2 ( x) 2 ( x) 2  2T ( T ) (T T ) (T T ) i, j 1 i, j i, j i, j 1 (2.3) y 2 ( y) 2 ( y) 2 Thay (2.2) và (2.3) vào phương trình vi phân (2.1) sẽ được : (T T ) (T T ) (T T ) (T T ) i 1 j i, j i, j i 1, j i, j 1 i, j i, j i, j 1 0 (2.4) ( x)2 ( y)2 (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 : 1 T (T T T T ) (2.5) i, j 4 i 1, j i 1, j i, j 1 i, j 1 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
  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 : a11T1 a12T2 a1nTn C1 a T a T a T C 21 1 22 2 21 n 2 (2.6) an1T1 an2T2 annTn Cn 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 : T  2T a (2.7)  x 2 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 : T T T p 1 T p i i (2.8)    Vế phải của (2.8) viết cho thời điểm sau (p+1)  2T ( T) (T p 1 T p 1 ) (T p 1 T p 1 ) i 1 i i i 1 (2.9) x 2 ( x)2 ( x) 2 thay (2.8) và (2.9) vào (2.7): T p 1 T p (T p 1 T p 1 ) (T p 1 T p 1 ) i i a i 1 i i i 1 (2.10)  ( x) 2 (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: a.  (T p 1 2T p 1 T p 1 ) T p 1 T p (2.11) ( x) 2 i 1 i i 1 i i a.  Đặt Fo sẽ được ( x)2 p 1 p p 1 p 1 p 1 Ti Ti Fo.(Ti 1 2Ti Ti 1 ) (2.12) vậy : p 1 p 1 p 1 p -FoTi 1 (1 2Fo)Ti Fo.Ti 1 Ti (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 5
  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  : p 1 p 1 qh h TK T1  (2.14) - Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian : k q (T p 1 T p 1)  (2.15) k x 2 1 Độ tăng nội năng dU phân tố sau thời gian  : x dU c . V (T p 1 T p ) c (T p 1 T p ) (2.16) 1 1 2 1 1 Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : k x h T p 1 T p 1  (T p 1 T p 1)  c (T p 1 T p ) (2.17) K 1 x 2 1 2 1 1 k h x  k  2 T p 1 T p 1 2 (T p 1 T p 1) T p 1 T p (2.18) c k ( x)2 K 1 c ( x)2 2 1 1 1 a.  h. x k Đặt Fo = , Bi , a ; Fo là tiêu chuẩn Phuriê, Bi là tiêu chuẩn Biô, a là hệ số khuyếch ( x)2 k c tán nhiệt độ sẽ được : p 1 p 1 p 1 p 1 p 1 p 2Bi.Fo TK T1 2Fo(T2 T1 ) T1 T1 Chuyển nhiệt độ và các đại lượng đã biết sang vế phải p 1 p 1 p 1 p 2Bi.Fo 2Fo 1 .T1 2Fo.T2 2Bi.Fo.TK T1 (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 : a11T1 a12T2 a1nTn C1 a T a T a T C 21 1 22 2 21 n 2 (2.20) an1T1 an2T2 annTn Cn trong đó: aij là các hệ số của nhiệt độ phải tìm, P 1 Ti là nhiệt độ cần tìm ở thời điểm (p+1), viết gọn của Ti 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 : aij  Ti  Ci  (2.21) trong đó: aij  là ma trận vuông gồm các hệ số của nhiệt độ phải tìm, Ti  là ma trận cột gồm nhiệt độ cần tìm ở thời điểm (p+1) Ci  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
  7. 1 Ti  aij  Ci  (2.22) 1 aij  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 : T T T p 1 T p i i (2.23)    Vế phải của (2.7) viết cho thời điểm trước p :  2T ( T) (T p T p ) (T p T p ) i 1 i i i 1 (2.24) x 2 ( x)2 ( x) 2 thay (2.23)và (2.24)vào (2.7) : T p 1 T p (T p T p ) (T p T p ) i i a i 1 i i i 1 (2.25)  ( x) 2 Để giải (2.25) cần biến đổi như sau : a.  T p 1 T p (T p 2T p T p ) ( 2.26) i i ( x) 2 i 1 i i 1 a.  Đặt Fo sẽ được : ( x)2 p 1 p p p p Ti Ti Fo.(Ti 1 2Ti Ti 1 ) Vậy : p 1 p p p Ti Fo.Ti 1 (1 2Fo.)Ti Fo.Ti 1 (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. 7
  8. Để các nghiệm hội tụ cần điều kiện : (1- 2Fo) 0 (2.28) tức là : 1 Fo 2 hay phải chọn bước thời gian đủ nhỏ : ( x) 2  (2.29) 2a 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  : p p qh h TK T1  (2.30) - Dòng nhiệt dẫn từ phân tố bên trong tới sau thời gian : k q (T p T p )  (2.31) k x 2 1 - Độ tăng nội năng dU phân tố dày x/2 sau thời gian  : x dU c . V (T p 1 T p ) c (T p 1 T p ) (2.32) 1 1 2 1 1 - Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : k x h T p T p  (T p T p )  c (T p 1 T p ) (2.33) K 1 x 2 1 2 1 1 Hay k h x  k  2 T p T p 2 (T p T p ) T p 1 T p (2.34) c k ( x)2 K 1 c ( x)2 2 1 1 1 a.  h. x k Với Fo = , Bi = , a = sẽ được : ( x)2 k c p p p p p 1 p 2Bi.Fo TK T1 2Fo(T2 T1 ) T1 T1 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 p 1 p p p T1 2Bi.Fo.TK (1 2Fo 2Bi.Fo)T1 2FoT2 (2.35) 8
  9. p 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 t1 theo nhiệt độ các điểm thời điểm trước. p 1 Điều kiện để xác định T1 , 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: T  2T  2T a. 2T a (2.37) 2 2  x 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() - Điều kiện biên loại 3 : T h T h 1 T ; 2 T x x 0 k x x a k T h T h 3 T ; 4 T (2.39) x y 0 k x y b 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
  10.  2T ( T ) (T T ) (T T ) T 2.T T i 1 j i, j i, j i 1, j i 1, j i, j i 1, j (2.40) x 2 ( x) 2 ( x) 2 ( x) 2  2T ( T ) (T T ) (T T ) T 2.T T i, j 1 i, j i, j i, j 1 i, j 1 i, j i, j 1 (2.41) y 2 ( y) 2 ( y) 2 ( y) 2 Riêng đạo hàm theo thời gian luôn có T T T p 1 T p i, j i, j (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 : T p 1 T p k T p 2.T p T p T p 2.T p T p i, j i, j i 1, j i, j i 1, j i, j 1 i, j i, j 1 (2.43) 2 2  c. ( x) ( y) 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 : T p 1 T p k T p 1 2.T p 1 T p 1 T p 1 2.T p 1 T p 1 i, j i, j i 1, j i, j i 1, j i, j 1 i, j i, j 1 (2.44) 2 2  c. ( x) ( y) (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 p 2.T p T p T p 2.T p T p k T p 1 i 1, j i, j i 1, j i, j 1 i, j i, j 1 .  T p (2.45) i, j 2 2 i, j ( x) ( y) 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ó: t p 1 2.t p 1 t p 1 t p 1 2.t p 1 t p 1 k i 1, j i, j i 1, j i, j 1 i, j i, j 1 . .  t p 1 t p (2.46) 2 2 i, j i, j ( x) ( y) 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
  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 : P P Dòng bức xạ là q R ( ) .I , với  là hệ số hấp thụ của vật, I là năng suất bức xạ chiếu tới - Điều kiên loại 3 : P P Dòng đối lưu từ không khí là qK ( ) h(TK Tm ) - Dòng nhiệt tổng : .I P P P P P P P P (2.47) q ( ) h(TK Tm ) .I h TK Tm h(TK Tm ) h trong đó : P P TK ,Tm 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 .I P là nhiệt độ quy đổi của bức xạ h .I P T P T P là nhiệt độ tương đương của không khí có kể đến bức xạ K K 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. + 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 : p 1 p 1 k p 1 p 1 k p 1 p 1 k p 1 p 1 k Ti 1, j Ti, j y T1 1, j Ti, j y Ti, j 1 Ti, j x Ti, j 1 Ti, j x  x x y y p 1 p c . x. y Ti, j Ti, j (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ên: p 1 p 1 k y p 1 p 1 k y p 1 p 1 k p 1 p 1 Ti 1, j Ti, j Ti 1, j Ti, j Ti, j 1 Ti, j x h TK Ti, j x  x 2 x 2 y y p 1 p c . x T T (2.49) 2 i, j i, j 11
  12. + 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 : p 1 p 1 k y p 1 p 1 y p 1 p 1 k x p 1 p 1 x Ti 1, j Ti, j h TK Ti, j Ti, j 1 Ti, j h TK Ti, j  x 2 2 y 2 2 x y p 1 p c . T T (2.50) 2 2 i, j i, j + 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 : p 1 p 1 k y p 1 p 1 k p 1 p 1 y Ti 1, j Ti, j . Ti 1, j Ti, j . y h T K Ti, j  x 2 x 2 p 1 p 1 k x p 1 p 1 x p 1 p 1 k 3 p 1 p Ti, j 1 Ti, j . h T K Ti, j Ti, j 1 Ti, j . x  c x y Ti, j Ti, j y 2 2 y 4 (2.51) k  h. x Sau khi lấy x = y , và đặt Fo , Bi , thay vào các phương trình trên sẽ được : c x 2 k Phương trình tại các phần tử thuộc nút bên trong : p 1 p 1 p 1 p 1 p 1 p Fo(Ti 1, j Ti 1, j Ti, j 1 Ti, j 1 ) (1 4)FoTi, j Ti, j (2.52) Phương trình tại các phần tử thuộc nút trên biên : p 1 p 1 p 1 p 1 p p 1 Ti 1, j Ti 1, j 2Ti, j 1 Fo 1 4Fo 2Bi.Fo Ti, j Ti, j 2Bi.Fo.TK (2.53) Phương trình tại các phần tử thuộc nút ở góc lồi : p 1 p 1 p 1 p p 1 2Fo(Ti 1, j Ti, j 1) 4Fo Bi 1 Ti, j Ti, j 4Bi.Fo.TK (2.54) Phương trình tại các phần tử thuộc nút ở góc lõm : 2 p 1 p 1 p 1 p 1 4 p 1 p 4 p 1 Fo(Ti 1, j 2Ti 1, j Ti, j 1 2Ti, j 1 ) 4Fo Bi.Fo 1 Ti, j Ti, j Bi.Fo.TK 3 3 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
  13. a11T1 a12T2 a1nTn C1 a T a T a T C 21 1 22 2 21 n 2 (2.56) an1T1 an2T2 annTn Cn Viết dạng ma trận : a11 a12 a13 a1n T1  C1  a21 a22 a23 a21 T2 C2   (2.57) an1 an2 an3 ann Tn  Cn  Các phương pháp giải thông dụng 1. Phương pháp định thức a11 a12 a13 a1n C1 a12 a13 a1n a a a a C a a a D 21 22 23 2n ;D 2 22 23 2n ; 1 an1 an2 an3 ann Cn an2 an3 ann a11 C1 a13 a1n a11 a12 a13 C1 a C a a a a a C D 21 2 23 2n ; ,D 21 22 23 2 ; 2 n an1 Cn an3 ann an1 an2 an3 Cn D D D Nghiệm sẽ là T 1 ;T 2 ; ;T n ; (2.58) 1 D 2 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 áp dụng tính chất 1 với (a1) áp dụng tính chất 2 với (b1) 2x + 2y = 4 (a1) (a1)/2 x + y = 2 (a2) hệ 2  hệ 1 x + y = 2 (a2) hệ 3  hệ 2 13
  14. x + 4y = 3 (b1) x + 4y = 3 (b1) (b1)-(a2) 0 + 3y = 1 (b2) (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 1 1 1 1 a11 a12 a1n T1  C1  1 1 1 T 1 a21 a22 a2n 2 C2   (1) a1 a1 a1 a1 31 32 33 3n 1 1 1 1 1 an1 an2 an3 ann Tn  Cn  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 đó: 1 1 1 1 1 1 1 1 a11 / a11 a12 / a11 a1n / a11 T1  C1 / a11  1 1 1 1 1 1 1 1 a / a a / a a / a T2 C / a 21 21 22 21 2n 21  2 21  a1 / a1 a1 / a1 a1 / a1 a1 / a1 / a1 31 31 32 31 33 31 3n 31 31 1 1 1 1 1 1 1 1 1 1 an1 / an1 an2 / an1 an3 / an1 ann / an1 Tn  Cn / an1  2 2 2 1 a12 a1n T1  C 1  2 2 T 2 1 a 22 a 2 n 2 C 2   (2) 1 a 2 a 2 a 2 32 33 3 n 2 2 2 2 1 a n 2 a n 3 a nn T n  C 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 2 2 2 2 2 1 a12 a1n T1  C1  1 a12 a1n T1  C1  2 2 2 2 T 2 2 3 3 T 3 1 1 (a22 a12) (a2n a1n ) 2 (C2 C1 ) 0 a22 a2n 2 C2     (3) 1 1 (a2 a2 ) (a2 a2 ) (a2 a2 ) 0 a 3 a 3 a 3 32 12 33 13 3n 1n 32 33 3n 2 2 2 2 2 2 2 2 3 3 3 3 1 1 (an2 a12) (an3 a13) (ann a1n ) Tn  (Cn C1 ) 0 an2 an3 ann Tn  Cn  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
  15. 1 a2 a2 T  C 2  2 2 2 12 1n 1 1 1 a12 a1n T1  C1  3 3 3 3 3 3 0 a / a a / a T C / a 4 4 4 22 22 2n 22 2 2 22 0 1 a23 a2n T2 C2   (4) 3 3 3 3 3 3 4 4   0 a32 / a32 a33 / a32 a3n / a32 0 1 a33 a3n 3 3 3 3 3 3 3 3 4 4 4 0 an2 / an2 an3 / an2 ann / an2 Tn  Cn / an2  0 1 an3 ann Tn  Cn  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) 2 2 2 2 2 2 1 a12 a1n T1  C1  1 a12 a1n T1  C1  4 4 4 4 4 4 0 1 a a T2 C 0 1 a a T2 C 23 2n  2  23 2n  2  (5) 0 1 1 a 4 a 4 a 4 a 4 0 0 a 5 a 5 33 23 3n 2n 33 3n 4 4 4 4 4 4 5 5 5 0 1 1 an3 a23 ann a2n Tn  Cn C2  0 0 an3 ann Tn  Cn  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 2 2 2 2 2 2 1 a12 a1n T1  C1  1 a12 a1n T1  C1  4 4 4 4 4 4 0 1 a a T2 C 0 1 a a T2 C 23 2n  2  23 2n  2  (6) 0 0 a 5 / a 5 a 5 / a 5 C 4 / a 5 0 0 1 a 6 C 6 33 33 3n 33 3 33 3n 3 5 5 5 5 5 5 6 6 0 0 an3 / an3 ann / an3 Tn  Cn / an3  0 0 1 ann Tn  C n  k g. Tiếp tục như vậy cho đến khi số hạng ann 1 , thì sẽ được tam giác sau 2 2 2 1 a12 a1n T1  C1  4 4 4 0 1 a a T2 C 23 2n  2  (7) 0 0 1 a 6 T C 6 3n 3 3 k 0 0 0 1 Tn  Cn  k h. Giải ra tính ngược từ dưới lên: hàng dưới cùng : Tn C n ; hàng chứa T có : T a6 T C 6 T a6 T C 6 ,. 3 3 3n n 3 3 3n n 3 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 1 1 1 a12 a13 a1n T1  C1  1 1 1 0 1 a a T2 C2 23 2n   0 0 1 a1 T C1 3n 3 3 1 0 0 0 1 Tn  Cn  15
  16. 1 a. Lấy hàng 2 làm gốc, nhân hàng 2 với a12 sẽ được: 1 2 2 2 0 a12 a23 a2n T2  C2  Lấy hàng 1 trừ đi hàng vừa có ở trên 1 1 1 2 1 2 1 2 2 2 2 1 0 a12 a12 a13 a23 a1n a2n T1  C1 C2  1 0 a13 a1n T1  C1  1 1 1 1 1 1 0 1 a a T2 C 0 1 a a T2 C 23 2n  2  23 2n  2  0 0 1 a1 T C1 0 0 1 a1 T C 1 3n 3 3 3n 3 3 1 1 0 0 0 1 Tn  Cn  0 0 0 1 Tn  Cn  (1) 1 1 2 2 b. Lấy hàng 3 làm gốc, nhân hàng 3 với a23 sẽ được: 0 0 a23 a3n T3 C3 ; Lấy hàng 2 trừ đi hàng vừa có 2 2 2 2 2 2 1 0 a13 a1n T1  C1  1 0 a13 a1n T1  C1  1 1 1 2 1 2 2 2 0 1 a a . a a T2 C C 0 1 0. a T2 C 23 23 2n 3n  2 3  2n  2  (2) 0 0 1 a1 T C1 0 0 1 a1 T C1 3n 3 3 3n 3 3 1 1 0 0 0 1 Tn  Cn  0 0 0 1 Tn  Cn  c. Tiếp tục như vậy sẽ được 2 2 2 1 0 a13 a1n T1  C1  2 2 0 1 0 a T2 C 2n  2  (4) 0 0 1 0 T C 2 3 3 1 0 0 0 1 Tn  Cn  2 2 d. Để triệt tiêu a13 của hàng 1, lấy hàng 3 làm gốc, nhân hàng 3 với a13 , rồi lấy hàng 1 trừ đi kết quả mới có. 3 3 3 1 0 0 a14 a1n T1  C1  2 2 0 1 0 a T2 C 2n  2  (5) 0 0 1 0 T C 2 3 3 1 0 0 0 1 Tn  Cn  3 3 e. Để triệt tiêu a14 của hàng 1, lấy hàng 4 làm gốc, nhân hàng 4 với a14 rồi lấy hàng 1 trừ đi kết quả mới có. 4 4 1 0 0 0 a1n T1  C1  2 2 2 0 1 0 a a T2 C 24 2n  2 (6) 0 0 1 0 T C 2 3 3 1 0 0 0 1 Tn  Cn  16
  17. 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 k k 1 0 0 0 0 T1  C1  T1  C1  k k 0 1 0 0 0 T2 C T2 C  2   2  (2.59) 0 0 0 1 0 T k T k 3 C3 3 C3 k k 0 0 0 1 Tn  Cn  Tn  Cn  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 T1 a21T2 a31T3 an1Tn ;(1) T a T a T a T : (2) 2 12 1 32 3 n2 n Tn a1nT1 a2nT2 an 1.nTn 1;(n) 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 - 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
  18. 1 Hình 2.4. Chia mạng tiết diện T T T T T , i, j 4 i1 i 2 i3 i4 ngang dầm bêtông 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) T4 = (T3 +100 + 80 +70 )/ 4 (4) T2 = (T1 + T3 + T5 + 100)/4 (2) T5 = (T2 + T6 + 50 + 40 )/ 4 (5) T3 = (T2 + T4 + T6 + 100)/4 (3) 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 : a11 a12 a13 a1n T1  C1  a21 a22 a23 a21 T2 C2   (2.60) an1 an2 an3 ann Tn  Cn  18
  19. Hay ở dạng gọn sau : [aij]. [Ti] = [Ci] (2.61) Từ đó sẽ rút ra được : - 1 [Ti] = [Ci] [aij] (2.62) - 1 trong đó [aij] là ma trận nghịch đảo của [aij] có dạng : b11 b12 b13 b1n b a a b a 1 21 22 23 21 (2.63) ij bn1 bn2 bn3 bnn 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à : T1 = b11C1 + b12C2+ b13C3 + b1nCn T2=b21C1+ b22C2+ b23C3 + b2nCn T3 = b31C1+ b32C2+ b33C3 + b3n Cn (2.64) Tn = bn1C1 bn2C2+ bn3C3 + bnnCn 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
  20. 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
  21. 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 T N1T1 N 2T2 N k Tk  N iTi (2.65) i 1 Hoặc ở dạng ma trận T1  T2 T N1 N2 Nk   N T (2.66) Tk  ở đâ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
  22. (2.67) K e Te f e ở đây: e là chỉ số biểu thị cho phần tử Te là nhiệt độ phải tìm tại các nút. K e 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 e là véc tơ phụ tải nhiệt hoặc nhiệt độ cho trước tại nút biên nào đó. 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ố 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ử Te dạng ma trận (2.67) ở trên thành hệ (n phần tử) cũng dưới dạng: K T f  (2.70) [ K ] là ma trận các hệ số của cả hệ T là véctơ ẩn của cả hệ f  là tải nhiệt tại các nút của cả hệ 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
  23. 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 Trong phần tử bậc nhất, nhiệt độ là hàm bậc nhất của toạ độ: T 1 2 x (2.71) 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à xi và x j thì nhiệt độ tương ứng tại đó là : Ti 1 2 xi ; T j 1 2 x j (2.72) 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: Ti  T N iTi N jT j N i N j   N  T  (2.73) T j  ở đây Ni và N j gọi là các hàm nội suy tại hai nút. 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 x x T j T i T (2.74) i j x j xi x j xi Theo (2.74) các hàm nội suy N i , N j là x j x x xi N N i N j  (2.75) x j xi x j xi Nếu lấy xi 0; x j l , thì x x N 1 (2.76) l l Lấy tổng hàm nội suy N i N j 1 (2.77) 23
  24. Từ (2.76) , (2.77) có thể thấy hàm nội suy Ni và N j là hàm bậc nhất theo x, biến đổi ngược chiều nhau có giá trị tại các vị trí khác nhau như sau, bảng 2.1 Bảng 2.1 Hàm Nút i Nút j x Ni 1 0 Giữa 0 và 1 Nj 0 1 Giữa 0 và 1 Ni +Nj 1 1 1 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 Ni rồi thay N j 1 N i sẽ được như sau: xi  x N i xi N j x j Ni N j   (2.78) 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 dN d 1 N  1 1 B (2.79) dx dx 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 Ti ,T j là ẩn số chưa biết phải tìm, nhưng trong một phân tố Ti ,T j có giá trị không đổi, nên nhiệt độ T trong phân tố chỉ phụ thuộc vào x , vậy gradient nhiệt độ ký hiệu g sẽ là dT dNi dN j dNi dN j Ti  g Ti T j  B T (2.80) dx dx dx dx dx T j  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
  25. 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) T N iTi N jT j 0 0 Với : Ti 100 C;T j 160 C; xi 0; x j 12cm; x 8cm. x j x 12 8 4 Và : Ni x j xi 12 0 12 x xi 8 0 8 N j x j xi 12 0 12 thay các giá trị trên vào (2.45) có : 4 8 T N T N T 100 160 156,6660 C = 156,6660C i i j j 12 12 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. 2 T(x) = 1 + 2x + 3x (2.81) 25
  26. 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 l phần tử. Trong mỗi phần tử có độ dài l x x , nếu lấy x 0 thì x ; x l . k i i j 2 k Nhiệt độ tại ba nút ứng với các toạ độ là 2 l l 2 Ti 1 ; Tj 1 2 3 ; Tk 1 2l 3l (2.82) 2 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 : 3x 2x2 x x2 x x2 T (x) 1 2 T1 4 4 2 Tj 2 2 Tk (2.83) l l l l l l 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: Ti  T (x) NiTi N jT j N kTk Ni N j N k  T j  N  T (2.84) Tk  Trong đó N i , N j và Nk 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 suy của phần tử một chiều bậc hai là 3x 2x 2 x x 2 x x 2 (2.85) N N i N j N k  1 2 4 4 2 2 2 l l l l l l 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 Giá trị của hàm nội suy Nút i j k Hàm Ni 1 0 0 nội suy Nj 0 1 0 Nk 0 0 1 Tổng: Ni+ Nj + Nk 1 1 1 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
  27. 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 dN dN j dN 3 4x 4 8x 1 4x B i k (2.86)   2 2 2 dx dx dx dx l l l l l 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 độ dT Ký hiệu gradient nhiệt độ: g , dx Ti  dT dN dN j dN dN i k (2.87) T j  T B T dx dx dx dx dx Tk  Gradient nhiệt độ là dT dN dN dN i T j T k T (2.88) dx dx i dx j dx k đạo hàm các biểu thức trong (2.70) thay vào (2.73) sẽ có dT 3 4x 4 x 1 x T 8 T 4 T (2.89) 2 i 2 j 2 k dx l l l l l l 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
  28. 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à 3x 2x2 3(l / 4) 2(l / 4)2 Ni 1 2 1 2 = 1-3/4 + 1/8 = 0,3750 l l l l x x2 (l / 4) (l / 4)2 N j 4 4 2 4 4 2 1 1/ 4 = 0,75 l l k l l x x2 l / 4 (l / 4)2 Nk 2 2 2 2 1/ 4 1/8 = -0,1250 l l l 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à : 3x 2x2 3(l /3) 2(l /3)2 Ni 1 2 1 2 1 1 2/9 0,2222 l l l l x x2 (l /3) (l /3)2 N j 4 4 2 4 4 2 4/3 4/9 = 0,8889 l l k l l x x2 l /3 (l /3)2 Nk 2 2 2 2 1/3 2/9 = -0,1111 l l l l Cũng thấy ngay rằng Ni + Nj + Nk = 0,2222 + 0,8889 – 0,1111 = 1 28
  29. 2.6.3. Phần tử hai chiều tam giác bậc nhất 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ố. y Tk k (xk,yk)  Ti  i (x ,y )  Tj j (xj,yj) 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ệ Ti 1 2 xi 3 yi ; T j 1 2 x j 3 y j ; Tk 1 2 xk 3 yk (2.91) Có thể giải ra các ẩn là các hệ số 1, 2 và 3 theo xi , x j , xk và Ti ,T j ,Tk bằng phương pháp định thức như sau. Viết (2.91) ở dạng hệ phương trình: 1 2 xi 3 yi Ti 1 2 x j 3 y j T j (2.92) 1 2 xk 3 yk Tk Ta có 1 x i y i D 1 x y x y x y x y x y x y x y (2.93) j j i j j i k i i k j k k j 1 x k y k thấy rằng D = 2A, với A là diện tích của tam giác Ti xi yi D T x y x y x y T x y x y T x y x y T (2.94) 1 j j j j k k j i k i i k j i j j i k Tk xk yk 29
  30. xi Ti yi D x T y x y x y T x y x y T x y x y T (2.95) 2 j j j k j j k i i k k i j j i i j k xk Tk yk yi xi Ti D y x T x y x y T x y x y T x y x y T (2.96) 3 j j j k j j k i i k k i j j i i j k yk xk Tk D D D Giải ra 1 ; 2 ; 3 sẽ là 1 D 2 D 3 D 1  x y x y T x y x y T x y x y T  1 2A j k k j i k i i k j i j j i k 1  y y T y y T y y T  (2.97) 2 2A j k i k i j i j k 1  x x T x x T x x T  3 2A k j i i k j j i k Với ký hiệu a i x j y k - x k y j ; b i y j - y k ; c i x k - x j a j x k y i - x i y k ; b j y k - yi ; c j x i - x k (2.98) a k x i y j - x j y i ; b k y i - y j ; c k x j - x i (2.97) trở thành 1 a T a T a T  1 2A i i j j k k 1 b T b T b T  (2.99) 2 2A i i j j k k 1 c T c T c T  3 2A i i j j k k Thay thế các giá trị của 1, 2 và 3 vào phương trình (2.95) sẽ có 1 1 1 T a T a T a T  b T b T b T x c T c T c T y (2.100) 2A i i j j k k 2A i i j j k k 2A i i j j k k sắp xếp lại như sau 1 T  a b x c y T a b x c y T a b x c y T  (2.101) 2A i i i i j j j j k k k k 1. Hàm nội suy 30
  31. 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 Ti  T (x, y) N iTi N jT j N k Tk Ni N j N k  T j  N  T (2.102) Tk  Từ (2.101) thấy các hàm nội suy là 1 N a b x c y i 2A i i i 1 N a b x c y (2.103) j 2A j j j 1 N a b x c y k 2A k k k Viết các hàm nội suy đày đủ theo toạ độ nút: 1 1 N a b x c y  x y - x y y - y x x - x y i 2A i i i 2A j k k j j k k j 1 1 N a b x c y  x y - x y y - y x x - x y (2.104) j 2A j j j 2A k i i k k i i k 1 1 N a b x c y  x y x y y y x x x y k 2A k k k 2A i j j i i j j i Để 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 1 2A N  x y x y y y x x x y  1 (2.105 a) i i 2A j k k j j k i k j i 2A - Tính Ni tại nút j có tọa độ x j ; y j 1 N  x y - x y y x - y x x y - x y  0 (2.105b) i j 2A j k k j j j k j k j j j - Tính Ni tại nút k có tọa độ là xk ; yk 1 N  x y - x y y x - y x x y - x y  0 (2.105c) i k 2A j k k j j k k k k k j k - Tính N j tại nút i có tọa độ là xi và yi 1 N  x y - x y x y - x y x y - x y  0 (2.105d) j i 2A k i i k i k i i i i k i 31
  32. - Tính N j tại nút j có tọa độ x j ; y j x y - x y x y - x y x y - x y N k i i k j k j i i j k j 1 (2.105e) j j xi y j x j yi xk yi xi yk x j yk xk y j - Tính N j tại nút k có tọa độ xk ; yk 1 N  x y - x y x y - x y x y - x y  0 (2.105g) j 2A k i i k k k k i i k k k - Tính Nk tại nút i có tọa độ xi và yi 1 N  x y x y x y x y x y x y  0 (2.105h) k i 2A i j j i i i i j j i i i - Tính Nk tại nút j có tọa độ x j ; y j 1 N  x y x y x y x y x y x y  0 (2.105h) k j 2A i j j i j i j j j j i j - Tính N k tại nút k có tọa độ xk ; yk xi y j x j yi xk yi xk y j x j yk xi yk N k k 1 (2.105k) xi y j x j yi xk yi xi yk x j yk xk y 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 Nút j Nút k 1 0 0 Ni 0 1 0 N j 0 0 1 N k 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 N i N j N k 1 (2.106) ở 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
  33. 1 N x a x b x x c x y i i 2A i i i i i i 1 N x a x b x x c x y (2.107) j j 2A j j j j j j 1 N x a x b x x c x y k k 2A k k k k k k Tính tổng N i xi N j x j N k xk 1  a x a x a x x b x b x b x y c x c x c x  (2.108) 2A i i j j k k i i j j k k i i j j k k 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 ai xi a j x j ak xk x j y k - x k y j xi x k y i - x i y k x j x i y j - x j y i xk 0 b bi xi b j x j bk xk y j - y k xi yk yi x j yi - y j xk 2A (2.109) c ci xi c j x j ck xk x k - x j xi x i - x k x j x j - x i xk 0 Thì sẽ thấy 1 1 N x N x N x (a bx cy) (0 2Ax 0) i i j j k k 2A 2A Bởi vậy rút ra x N i xi N j x j N k xk (2.110) Tương tự như vậy, từ (2.107) cũng có 1 N y a y b y x c y y i i 2A i i i i i i 1 N y a y b y x c y y (2.111) j j 2A j j j j j j 1 N y a y b y x c y y k k 2A k k k k k k Tính tổng N i yi N j y j N k yk 1  a y a y a y x b y b y b y y c y c y c y  2A i i j j k k i i j j k k i i j j k k 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
  34. a ai yi a j y j ak yk x j y k - x k y j yi x k y i - x i y k y j x i y j - x j yi yk 0 b bi yi b j y j bk yk y j - y k yi yk yi y j yi - y j yk 0 (2.112) c ci yi c j y j ck yk x k - x j yi x i - x k y j x j - x i yk 2A Thì sẽ thấy 1 1 N y N y N y (a bx cy) (0 0x 2Ay) i i j j k k 2A 2A Rút ra y N i yi N j y j N k yk (2.113) 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 x N i xi N j x j N k xk y N i yi N j y j N k yk (2.114) 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 Ni N j N k x x x x 1 bi b j bk B B N N N N i j k 2A ci c j ck y y y y 4. Gradient nhiệt độ Gradient nhiệt độ xác định bằng T N i N j N k Ni N j N k T T T Ti  i j k g x x x x x x x T T N N N N N N j  i T j T k T i j k i j k Tk  y y y y y y y (2.115) Ti  1 bi b j bk T j  B T c c c 2A i j k Tk  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
  35. đó đế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à PO L (2.116) i 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à A 0,5.(PO).( jk) PO L i (2.117) i A 0,5.(QR).( jk) 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 A A i j k 1 (2.118) A A 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
  36. 1 L a b x c y i 2A i i i 1 L a b x c y j 2A j j j 1 L a b x c y (2.121) k 2A k k k 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à ai x jyk - x k y j; bi y j - yk ; ci x k - x j a j x k yi - x i yk ; b j yk - yi ; c j x i - x k a k x i y j - x jyi ; bk yi - y j ; ck x j - x 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à a!b!c! La Lb Lc dA N a N b N c dA 2A (2.123) i j k i j k A A a b c 2 ! ở đâ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) y 3 (x3,y3) 1 2 (x ,y ) 1 1 (x2,y2) x 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
  37. 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 1 N (b x)(a y) 1 4ab 1 N (b x)(a y) 2 4ab 1 N (b x)(a y) (2.126) 3 4ab 1 N (b x)(a y) 4 4ab 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
  38. N x 1 (a y) (a y) (a y) (a y) B N (2.127) 4ab (b x) (b x) (b x) (b x) y 3. Gradient nhiệt độ Từ (2.118), gradient nhiệt độ được viết là T y x 2 4 T x (2.128) y 3 4 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 N N N N 1 T 2 T 3 T 4 T x x 1 x 2 x 3 x 4 1  (a y)T (a y)T (a y)T (a y)T  (2.129) 4ab 1 2 3 4 tương tự có T 1  (b x)T (b x)T (b x)T (b x)T  y 4ab 1 2 3 4 Ma trận gradient nhiệt độ là T T  1  x 1 a y) (a y) (a y) (a y) T2 g T   B.T (2.130) 4ab (b x) (b x) (b x) (b x) T3 y  T4  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
  39. 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 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 T N1T1 N2T2 NmTm N T 39
  40. 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 độ x N x N x N x N x 1 1 2 2 m m y N1 y1 N2 y2 Nm ym N y 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 1 N 1  1 1 1 2  ( 1) 1 N 1  (2.131) 2 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à 1 N  1  1   (2.132) 2 Nhiệt độ : 1 T  1  T 1  T , 2 1 2 tức là: T N1T1 N2T2 (2.133) Tọa độ: từ (2.131) rút ra:  1 2N1 1 N1 (1 N2 ) N1 N2 , hay  N11 N2 2 (2.134) 40
  41. Đạo hàm hàm nội suy [B] theo biến x: Từ (2.132) có dN 1  1 1 (2.135) d 2 Mặt khác theo đạo hàm hợp thì dN dN dx dN dx J ; với J gọi là Jacobian của phép biến đổi tọa độ. Từ đó suy ra d dx d dx d dN dN 1 B J 1 J 1  1 1 dx d 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:  0  1  N 1  i 1 0) 1 1 2  ( 1)  1 N 1  2 j 0 ( 1) 0 1  ( 1)  0  N 1  (2.136) k 1 ( 1) 1 0 2 ở đâ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 x xi Ni x j N j xk Nk (2.137) Đạo hàm N theo  ở trên viết dưới dạng hàm hợp qua x là dN dN dx i i d dx d dN dN dx j j (2.138) d dx d dN dN dx k k d dx d 41
  42. dx với J gọi là Jacobian của phép biến đổi tọa độ. Từ đó suy ra d dN dN i J 1 i dx d dN dN j J 1 j dx d dN dN k J 1 k (2.139) dx d Theo (2.137) thì Jacobian của phép biến đổi phần tử một chiều bậc hai là dx dN dN j dN J  i x x k x d d i d j d k d  d 2 d  J  1  xi 1  x j 1  xk d 2 d d 2 1 1 J   xi ( 2 )x j  xk (2.140) 2 2 Đạ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: dN  dNi  i dx d dN dN j 1 dN j B  J   dx dx d dN k dNk dx  d  1   1 2 1 1  xi ( 2 )x j  xk 2  (2.141) 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
  43. dx dN dN j dN J  i x x k x d d i d j d k 1 1  2 2 4  6 2 2 2 8 8 2 Nghịch đảo Jacobian là 1 J  1 2 Đạo hàm hàm nội suy là dN  dNi  i 1  1   dx d  2 4 2 dN j 1 dN j 1 B  J   2    (2.142) dx d 2 dN 1 1  k dNk  2  4 2  dx  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 N N x N y i x, y i i  x  y  N N x N y i x, y i i (2.144)  x  y  (2.144) có thể viết dạng ma trận N i  x y Ni  N i     x x J  (2.145) N  x y N  N  i i i     y  y  từ đó suy ra đạo hàm hàm nội suy theo tọa độ gốc 43
  44. N i  N i  x 1  J  N  N  i i y    Ma trận Jacobian [J] x y   J  (2.146) x y   Nghịch đảo của ma trận Jacobian [J]-1 được tính theo y y 1 1   J  (2.147) x x detJ    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à N1 L1 1-  N2 L2  ;0  1 N3 L3 ;0  1 (2.148) 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
  45. Đố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 T  N iTi (2.151) i 1 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( ,)  Ni ( ,).xi i 1 8 y( ,)  N i ( ,).yi (2.152) 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 N 1  1  1   1 4 45
  46. 1 N 1  2 1  2 2 1 N 1  1    1 3 4 1 N 1  1  2 4 2 1 N 1  1    1 5 4 1 N 1  2 1  6 2 1 N 1  1    1 7 2 1 N 1  1  2 (2.153) 8 2 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. K T f  (2.154) Để 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 T  T  T  T k x k y k z qV  x x y y z 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  2T  2T  2T 0 x 2 y 2 z 2 - Bài toán dẫn nhiệt ổn định một chiều qua thanh được biểu thi bởi 46
  47. d 2T kA hP(T T ) 0 (2.155) dx 2 a 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
  48.  T  T  T (2.156) k x k y k z qV 0 x x y y z z Cùng với điều kiện biên T = Tb trên mặt S1 T T T k i k j k k q 0 trên mặt S2 x x y y z z T T T k i k j k k h(T T ) 0 trên mặt S3 (2.157) x x y y z z a ở đây i, j 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 độ 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 I f (t, y, y')dt (2.158) 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 f d f 0 .” (2.160) y dt y' dy Ở đây y’ là đạo hàm của y theo thời gian t: y' . 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 f d f 0 (2.161) y dx yx Phương trình Euler- Lagrange tổng quát đối với ba biến độc lập f  f  f  f 0 (2.162) u x ux y u y z uz 48
  49. 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 2 T T T 1 2 I(T) 1 k k k 2q T d qTds h T T ds (2.163) 2 x x y x z x V 2 a  S 2 S 2 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 e T  N iTi N T  (2.164) i 1 ở đây [N] = [N1, N2 , , Nm ] là ma trận các hàm hình dạng, và T1  T2 T   (2.165) Tm  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 I I(T )  0 (2.166) i 1 Ti ở đâ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 I 0 với i = 1, 2, , n (2.167) Ti 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 I T  I e T e (2.168) i 1 49
  50. 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 I I e 0 (2.169) 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à I e  I e  0 với i = 1, 2, , m (2.170) T  Tj 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 e 2 e 2 e 2 T T T 1 2 I e 1 k k k 2q T e d qT e ds h T e T ds (2.171) 2 x x y x z x V 2 a  S 2 S 2 Ma trận gradient nhiệt độ : e T  N1 N 2 N m T  x x x x 1 e T N N N T2 g  1 2 m  B T (2.172) y y y y e T N1 N 2 N m T m  z  z z z Ba số hạng đạo hàm đầu: 2 2 2 T e T e T e kx k y kz x y z T e  k 0 0 x e e e x e T T T  T T 0 k 0 g D g (2.173)  y      x y z  y 0 0 k e z T z  Theo quy tắc của ma trận chuyển vị của tích gT TT BT 50
  51. Nên gT D g TT BT DB T (2.174) Nhiệt độ trong phần tử : T1  e T2 T N  T N1, N2 , ,Nm   N1T1 N2T2 NmTm (2.175) Tm  thay thế (2.174) và (2.175) vào phương trình (2.171), ta có 1 1 I e T T B T D B T 2q N T d q N T ds h N T T 2 ds (2.176)        V          a 2  S 2e S 3e 2 I e bây giờ thực hiện cực tiểu hóa tích phân 0  T I e 1   T T B T D B T d q N T d        V     T  2  T    T   (2.177)   1 2 q N T ds h N T T ds 0       a  T  S 2e  T S 3e 2 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). T e  N T Đạo hàm của nhiệt độ trong phần tử theo nhiệt độ các nút được tính theo (2.175) như  T  T sau T e N1 T1 T e N2 T2 . T e Nr Tr hay 51
  52. N1  e T N2 T  N N (2.178)  T Nm  Tính riêng từng số hạng của phương trình trên như sau 1  1 (1) TT BT DB T d 2BT DB T d (2.179) 2  T 2  1  (2) 2q N T d q N T d (2.180)  V    V   2  T    (3) qN T ds qNT ds (2.181)  TS 2e S 2e  1  1 (4) h N T T 2 ds h N T 2 2 N T T T 2 ds    a       a a  T S 3e 2  T S3e 2 h N T T N ds h N T T ds (2.182)       a S3e S3e Thay các kết quả trên vào (2.177) được I e BT DB T  d q N T d qN T ds  T  V   S 2e h N T T N ds h N T T ds 0       a S 3e S 3e 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 B T D B T d h N T N T ds q N T d q N T ds h N T T ds            V       a  S 3e  S 2e S 3e (2.183) Và viết dạng gọn hơn là K T f  (2.184) trong đó K  BT DBd hN T N ds (2.185)  S 3 52
  53. f q N T d q N T ds hT N T ds (2.186)  V     a    S 2 S 3 [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:  L T d 0 (2.187) k  ở đâ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, L T là phương trình vi phân chủ đạo, T là nghiệm xấp xỉ. Điều đó nghĩa là  T  T  T N k k k q d 0 (2.188) k x x x y y y z z z V  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à  T  T Nk Nm N k d N k ds k T m d (2.189) k x k x x  x x S x x  x x ở đâ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 N k N m N k N m N k N m k k k T m d q N d N qdS x x x x x x x x x V k k   S hN N T m dS hT N dS 0 k m a k S S (2.190) Gộp các hệ số của nhiệt độ nút T m  lại với nhau sẽ được N N N N N N K k k m k k m k k m d hN N ds (2.191) km x y z k m  x x y y z z S 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
  54. f q N d qN ds hT N ds (2.192) k V k k a k  S S sẽ dẫn tới Kkm  T m fk  (2.193) Tức là K T f  (2.194) 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. 1  2 q 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à T N1T2 N1T2 (2.195) x2 x x1 x N1 và N 2 (2.196) x2 x1 x2 x1 + Ma trận độ cứng Ma trận độ cứng của phần tử theo (2.185) đã biết là 54
  55. T T K e B DBd hN NdA (2.197)  As 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 - Tính số hạng đầu của [K]e : K e1 B DBd ,  trong đó các ma trận [B], [D], [N] xác định như sau Chọn tọa độ x1 0; x2 l , thì hàm nội suy là : x x N 1 và N (2.198) 1 l 2 l [D] ma trận hệ số dẫn nhiệt : [D] = k 1 T 1 1 Đạo hàm của hàm nội suy [B] : B  1 1; nên B l l 1 1 1 1 k 1 1 T (2.199) B DB k  1 1 2 l 1 l l 1 1 Vậy số hạng đầu [K]e1 T x l l k 1 1 Ak 1 1 B DBd BT DBA.dx Adx (2.200)  x 0 0 2 l 1 1 l 1 1 T - Tính số hạng sau của [K]e : K e2 hN N dAS . As 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  N1 2 N 2 2  0 1 (2.201) Nên T 0 0 0 hN NdAS h 0 1dA hA (2.202) As A 1 0 1 Vậy ma trận độ cứng của phần tử Ak Ak Ak 1 1 0 0 l l K  hA (2.203) e l 1 1 0 1 Ak Ak hA l l - Tính véc tơ phụ tải nhiệt {f}. 55
  56. Theo (2.186): f q N T d q N T ds hT N T ds  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 i )i (N j )i ] 1 0 - 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 [N] [(N i ) j (N j ) j ] 0 1 Vậy véc tơ phụ tải nhiệt {f} là 1 0 1 0 qA f q N T ds hT N T ds q dA hT ds qA hT A    a   a a S 2 S3 A 0 A 1 0 1 hTa A (2.204) b. Phương trình đặc trưng của phần tử Phương trình đặc trưng K T f  sẽ là Ak Ak l l Ti  qA  (2.205) Ak Ak T hT A hA j  a l l 0 2 0 Ví dụ 2.6. Cho bài toán trên với số cụ thể sau: l = 4 cm, k = 0,5 W/m C, q = 100 W/m , Ta = 40 C, h = 20 W/m2 0C; lấy A = 1m2 , thay số vào được: 12,50 -12,50 Ti  100   -12,50 32,50 T j  800 Giải ra T1  53   T2  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
  57. 1 2  3  4     q h Ta l l l 1 2 3 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à k A k A 1 1 l l qA K 1 1 ; f  (2.206) 1 k A k A 1  1 1 0  l1 l1 Phương trình đặc trưng của phần tử 1 k1 A k1 A l l T1  qA 1 1   k1 A k1 A T 0 2   l1 l1 - Phần tử 2 - (lớp 2) Ma trận nhiệt dẫn và véc tơ phụ tải nhiệt là k A k A 2 2 l l 0 K  2 2 ; f  2 k A k A 2  2 2 0 l2 l2 Phương trình đặc trưng của phần tử 2 k2 A k2 A l l T2  0 2 2   (2.207) k2 A k2 A T 0 3   l2 l2 57
  58. - Phần tử 3 - (lớp 3) k A k A 3 3 l l 0  K  3 3 ; f  (2.208) 3 k A k A 3  3 3 hATa  hA l3 l3 Phương trình đặc trưng của phần tử 3 k A k A 3 3 l l T  0  3 3 3 (2.209) k A k A   3 3 T4  hATa  hA l3 l3 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) k1A k1 A T1 T2 qA (nút 1) (2.210) l1 l1 k1 A k1 A T1 T2 0 (nút 2) (2.211) l1 l1 Từ (2.206) k2 A k2 A T2 T3 0 (nút 2) (2.212) l2 l2 k2 A k2 A T2 T3 0 (nút 3) (2.213) l2 l2 Từ (2.207) k3 A k3 A T3 T4 0 (nút 3) (2.214) l3 l3 k A k A 3 3 (nút 4) (2.215) T3 hA T4 hATa l3 l3 Theo nguyên tắc trên thì : 58
  59. k1A k1 A - Nút 1: giữ nguyên (2.210): T1 T2 qA l1 l1 k A k A k A k A - Nút 2: (2.211) + (2.212): 1 1 2 2 T1 T2 T3 0 l1 l1 l2 l2 k A k A k A k A - Nút 3: (2.213) + (2.214) : 2 2 3 3 T2 T3 T3 T4 0 l2 l2 l3 l3 k A k A - Nút 4: giữ nguyên (2.215): 3 3 T3 hA T4 hATa l3 l3 Để 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 k A k A 1 1 0 0 l l 1 1 k1 A k1 A k2 A k2 A T1  qA  0 l1 l1 l2 l2 T2 0 (2.216) k A k A k A   2 2 3 T3 0 0 0 l2 l2 l3 hAT T4  a  k3 A k3 A 0 0 hA l3 l3 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à d 2T q V 0 (2.217) dx 2 k 59
  60. 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 q T (x) V L2 x 2 T (2.218) 2k m Để 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 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 Ke và f  - 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 Ke B DBd (2.219)  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 kA 1 1 K B D B d (2.220)  e       l 1 1 - 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 f q N T d (2.221)  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à N N N i i N i j j i j j 1 0 0 1 (2.222) N N i N j  2 2 2 2 Do đó 60
  61. 1 T 1 1 N 1 1 ; N (2.223) 2 2 1 vậy x l T 1 1 qV Al 1 f  qV N d qV A.dx  (2.224) x 0  2 1 2 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ó Ke và f  như nhau kA k A q Al  i V kA 1 1 Ti  q Al 1 Ti   V  l l  2  (2.225) T kA kA T q Al l 1 1 j  2 1 j  V l l 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 kA kA q Al  0 0 0 V l l 2 kA kA kA kA T qV Al qV Al 0 0 1  l l l l T 2 2 2 kA kA kA kA qV Al qV Al (2.226) 0 0 T3   l l l l 2 2 T4 kA kA kA kA qV Al qV Al 0 0 T l l l l 5  2 2 kA kA q Al 0 0 0 V l l 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); 0 3 0 k =12 (W/m C) ; qV = 200000 W/m ; nhiệt độ hai mặt ngoài Tm =30 C 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 0 T1  750  -1600 3200 - 1600 0 0 T2 1500 0 -1600 3200 -1600 0 T3  1500 0 0 -1600 3200 -1600 T 1500 4 0 0 0 -1600 1600 T5  750  61
  62. 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 1600 1600 T3  750 1600.T3 1600.T4 750 (3a)   (2.227) 1600 1600 T4  750 1600.T3 1600.T4 750 (3b) Phương trình ma trận phần tử 4 cũ 1600 1600 T4  750 1600.T4 1600.T5 750 (4a)   (2.228) 1600 1600 T5  750 1600.T4 1600.T5 750 (4b) để T5 = 30 thì (4b) phải là : 0.T4 .T5 30 (4b)' Khi đó ( 4a) sẽ là : 1600.T4 1600.30 750 (4a)' Vậy phần tử 4 sẽ có : 1600.T 0.T 750 1600.30 48750 (4a)' 4 5 (2.229) 0.T4 .T5 30 (4b)' (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 T1  37,0313 1600 -1600 0 0 0 T1  750  -1600 3200 - 1600 0 0 T 1500 T2 36,5625 2 → Giải ra T 35,1563 (2.230) 0 -1600 3200 -1600 0 T3  1500  3   0 0 -1600 3200 0 T 48750 T 32,8125 4 4 0 0 0 0 1 T5  30  T5  30,0000 So sánh với nghiệm giải tích Bảng 2.4 Phương PTHH Giải tích pháp T1 37,0313 37,5000 T2 36,5625 37,0313 T3 35,1563 35,6250 T4 32,8125 33,2813 T5 30,0000 30,0000 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
  63. T = NiTi + NjTj + NkTk (2.231) Các hàm nội suy đã có trong phần trước : 2 2 4x 4x 2 x 2 x 3x 2 x ; N ; N (2.232) N i 1 2 j 2 k 2 l l l l l l Đạo hàm của hàm nội suy [B] đã xác định được 4x 3 4 8x 4x 1 B (2.233)   2 2 2 l l l l l l a. Ma trận độ cứng Để tính ma trận độ cứng là K BT DBd , cần phải xác định tích số  BT DB Trong đó [D] = k; và chú ý các phép nhân ma trận sau: 3 1 1.3 1.4 3 4 1 2 1.3 2.4 11; và 3 4 (2.234) 4 2 2.3 2.4 6 8 Nên 4x 3 2 l l T 4 8x 4x 3 4 8x 4x 1 B DB k 2 2 2 2 l l l l l l l l 4x 1 2 l l 2 4x 3 4x 3 4 8x 4x 3 4x 1 2 2 2 2 2 l l l l l l l l l l 2 4 8x 4x 3 4 8x 4 8x 4x 1 (2.235) 2 2 2 2 2 l l l l l l l l l l 2 4x 1 4x 3 4x 1 4 8x 4x 1 2 2 2 2 2 l l l l l l l l l l + Tính ma trận độ cứng K BT DBd  63
  64. 16x 2 24x 16x 32x 2 24x 16x 2 4x 12x 9 12 3 2 2 2 l l l l l l l l 1 16x 32x 2 24x 64x 64x 2 16x 32x 2 8x Ak 12 16 4 dx 2 2 2 2 x l l l l l l l l l 16x 2 12x 4x 16x 32x 2 8x 16x 2 8x 3 4 1 2 2 2 l l l l l l l l Sau khi lấy tích phân có: l 16x 3 24x 2 40x 2 32x 3 16x 3 16x 2 9x 12x 3x 2 2 2 3l 2l 2l 3l 3l 2l 1 40x 2 32x 3 64x 2 64x 3 24x 2 32x 3 K Ak 12x 16x 4x   2 2 2 2 l 2l 3l 2l 3l 2l 3l 16x 3 16x 2 24x 2 32x 3 16x 3 8x 2 3x 4x x 3l 2 2l 2l 3l 2 3l 2 2l 0 Thay cận sẽ được: 16 32 16 12 9 20 12 8 3 3 3 3 14 16 2 1 32 64 32 Ak K Ak 20 12 16 32 12 4 16 32 16 l 3 3 3 6l 2 16 14 16 32 16 8 3 12 4 4 1 3 3 3 Cuối cùng có ma trận độ cứng của phần tử bậc hai một chiều 14 16 2 Ak K 16 32 16 (2.236) 6l 2 16 14 b. Véc tơ phụ tải T T Theo (2.186) :  f  qV N d , thay [N] vào sẽ được  3x 2x 2 1 2 l l 2 T 4x 4x (2.237)  f  qV N d qV Adx  l 2  l x 2x 2 2 l l lấy tích phân sẽ được : 64
  65. L 3x2 2x3 x 3 2 2l 3l 2 l l l 2 3 6 9 4 1 4x2 4x3 4 l q Al f q A q A 2l l q A 12 8 V 4 (2.238)   V 2 V V 2l 3l 3 6 6 2 3 l 2l 3 4 1 x 2x 2 2 3 2l 3l 0 c. Phương trình ma trận đặc trưng của phần tử 14 16 2 Ti  1 Ak qV Al 16 32 16 T 4 (2.239) j  6l 6 2 16 14 Tk  1 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 Ti  1 14 16 2 Ti  1 Ak qV Al 16 32 16 T 4 =133,33 16 32 16 T 500 4 j  j  6l 6 2 16 14 Tk  1 2 16 14 Tk  1 1866,7 - 2133,3 266,7 T1  500  = - 2133,3 4266,7 - 2133,3 T 2000 2   266,7 - 2133,3 1866,7 T3  500  Áp đặt điều kiện biên: Do T3 =30, thay vào, hệ trở thành 933,33 1066,7 0 T1  - 333,33 44,1701  -1066,7 2133,4 0 T 36001 . giải ra được T 38,9600 2     0 0 1 T3  30  30  + 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
  66. Phần tử 1 1866,7 - 2133,3 266,7 T1  500  - 2133,3 4266,7 - 2133,3 T 2000 2   266,7 - 2133,3 1866,7 T3  500  Phần tử 2 1866,7 - 2133,3 266,7 T3  500  - 2133,3 4266,7 - 2133,3 T 2000 4   266,7 - 2133,3 1866,7 T5  500  Lắp ghép được 1866,7 - 2133,3 266,7 0 0 T1  500  - 2133,3 4266,7 - 2133,3 0 0 T2 2000 266,7 - 2133,3 1866,7 1866,7 - 2133,3 266,7 T3  500 500 0 0 - 2133,3 4266,7 - 2133,3 T 2000 4 0 0 266,7 - 2133,3 1866,7 T5  500  Áp đặt điều kiện biên T5 =30, hệ trở thành 37,4732 1866,7 - 2133,3 266,7 0 0 T1  500   - 2133,3 4266,7 - 2133,3 0 0 T 2000 37,0070 2 T= 266,7 - 2133,3 3733,4 - 2133,3 266,7 T3  1000  35,6051 0 0 - 2133,3 4266,7 0 T 2000 2133,3.30 65999 33,2705 4 0 0 0 0 1 T5  30  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 nhất PTHH bậc hai (5 phần tử) 1 phần tử 2 phần tử T1 37,5000 37,0313 44,1701 37,4732 T2 37,0313 36,5625 37,0070 T3 35,6250 35,1563 38,9600 35,6051 T4 33,2813 32,8125 33,2705 T5 30,0000 30,0000 30,00 30,0000 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
  67. 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 2 Thể tích phần tử khảo sát là  = (r2 - r1 ) 1, vi phân thể tích là d = 2 rdr. Như vậy biến số độc lập 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 T N1T1 N 2T2 (2.240) 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à r r N  N1 N 2  1 (2.241) l 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 1 B  1 1 (2.242) l Toạ độ r được biểu thị qua hàm nội suy r N1r1 N 2 r2 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
  68. T T K  B DBd hN NdAs (2.243)  As T + Tính số hạng thứ nhất : B DBd :  tích số BT DB đối với phần tử một chiều bậc nhất, ta đã biết là 1 1 1 k 1 1 T B DB k  1 1 2 l 1 l l 1 1 2 2 và với vách trụ ở đây l = (r2 – r1) ; nên 2 r 2 T r2 2 k 1 1 2 k 1 1 r B DBd 2 rdr 2  ri l 1 1 1 1 2 r2 r1 r1 Sau khi thay cận có T 2 k r r 1 1 B DBd 1 2 (2.244)  l 2 1 1 T - Tính số hạng thứ hai : hN NdAs : 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ó T 0 0 0 hN N dAS h 0 1dA 2 r2 h (2.245) As A 1 0 1 Vậy ma trận độ cứng [K] là 2 k r1 r2 1 1 0 0 K  2 .r2 h (2.246) l 2 1 1 0 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 0 f  hTa N  dAS hTa 2 .r2  As 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
  69. 2 k r1 r2 1 1 0 0 T1  0 2 .r2 h  hTa 2 .r2  (2.247) l 2 1 1 0 1 T2  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: 0 0 2 0 0 r1 = 40 cm, r0 = r2 = 60 cm, k = 10W/m C, Tm1 =100 C, h = 10W/m C, Ta = 30 C. + 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 2 k (r1 r2 ) 1 1 0 0 K e 2 .r2 h l 2 1 1 0 1 2 .10 (0,6 0,4) 1 1 0 0 50 50 . 2 .0,6.10 0,2 2 1 1 0 1 50 62 0 0 0  f  hTa 2 .r0  10.30.2 .0,6   1 1 360 Phương trình ma trận đặc trưng của phần tử là 50 50 T1  0    (2.248) 50 62 T2  360 0 Áp đặt điều kiện biên : T1 = 100 C sẽ có 1 0 T1  100  0   T2 = 86,45 C 0 62 T2  360 50.100 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 2 .k r1 r2 1 1 2 .10 0,4 0,5 1 1 90 90 K 1 l 2 1 1 0.1 2 1 1 90 90 69
  70. 0 - Véc tơ tải f 1  0 90 90 T1  0 - Phương trình ma trận đặc trưng :   90 90 T2  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 : 2 k (r2 r3 ) 1 1 0 0 K 2 2 .r0 h l 2 1 1 0 1 2 .10 (0,5 0,6) 1 1 0 0 110 110 . 2 .0,6.10 0,1 2 1 1 0 1 110 122 - Véc tơ tải 0 0 0  f 2 hTa 2 .r3  10.30.2 .0,6   1 1 360 110 110 T2  0  - Phương trình ma trận đặc trưng:   110 122 T3  360 + Lắp ráp phương trình đặc trưng tổng thể : 90 90 0 T1  0  90 90 110 110 T 0 (2.249) 2   0 110 122 T3  360 Hay gọn lại là 90 90 0 T1  0  90 200 110 T 0 (2.250) 2   0 110 122 T3  360 0 Áp đặt điều kiện biên : do T1 = 100 C, nên 1 0 0 T1  100  T1  100  0 200 110 T 0 90.100 giải ra T 92,4878 (2.251) 2   2   0 110 122 T3  360  T3  86,3415 70
  71. 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 2 k ri rj 1 1 0 0 K  2 .r0 h (2.252) l 2 1 1 0 1 2. Véc tơ phụ tải nhiệ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 q N T 2 .rdr , nên V   r T T f  hTa N dAS qV N 2 .rdr (2.253) As r T 0 - Số hạng đối lưu đã biết là hTa N dAS hTa 2 .r2  As 1 - Tính số hạng nguồn trong ký hiệu  f qV f q N T 2 .rdr (2.254)  qV V   r Biến số độc lập r trong tọa độ trụ được biểu thị bằng r N i ri N j r j (2.255) Trong đó Ni và Nj là các hàm nội suy: x x N 1 ; N (2.256) i l j l thay (2.229) vào (2.228) sẽ được 2 Ni Ni ri Ni N j rj   f qV 2 qV .(Ni ri N j rj )dr 2 qV dr (2.257) N 2 r j r N j Ni ri N j rj  Để tính biểu thức trên, cần áp dụng công thức tích phân : a b a!b! N i N j dl , (2.258) l (a b 1)! 71
  72. Với N i ; N j 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 1 1 1!1! 1 2 2!0! l N i N j dl và Ni dl (2.259) l (1 1 1)! 6 l (2 0 1)! 3 Thực hiện tích phân số hạng nguồn trong (2.257) rj r r ri rj 2 q 2ri rj  2 q 2ri rj   f  2 q . 3 6 V r .rj V l (2.260) qV V r r  ri  6 ri 2rj  6 ri 2rj  ri rj 6 3 ri Vậy véc tơ lực của phân tố là 2r r T T 0 2 qV i j  f  hTa N  dAS qV N 2 .rdr hTa 2 .r2  l  (2.261) As r 2r r 1 6 i j  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à r r T 2r r 2 k i j 1 1 0 0  i  0 2 qV i j  2 .r0h   hTa 2 .r2  l  (2.262) l 2 1 1 0 1  T j  1 6 ri 2rj  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. 2 1 2  3  4  5 h = 4000W/m 0 Ta = 20 C Tâm thanh trụ Mặt ngoài 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