Bài giảng Lớp cao học cơ khí phương pháp sai phân hữu hạn giải bài toán dẫn nhiệt - Trịnh Văn Quang

Nội dung text: Bài giảng Lớp cao học cơ khí phương pháp sai phân hữu hạn giải bài toán dẫn nhiệt - Trịnh Văn Quang

1. Bài giảng Lớp cao học Cơ khí PHƯƠNG PHÁP SAI PHÂN HỮU HẠN giải bài toán dẫn nhiệt PGS.TS. Trịnh Văn Quang- ĐHGT 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). 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 (1) x 2 y 2 1
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) x 2 ( x) 2 ( x) 2  2T ( T ) (T T ) (T T ) i, j 1 i, j i, j i, j 1 (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 (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 ) (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 . 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 (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
3. 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 (7)  x2 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 (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 (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 (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 (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 ) (12) vậy : p 1 p 1 p 1 p -FoTi 1 (1 2Fo)Ti Fo.Ti 1 Ti (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 3
4. 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  (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)  (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 ) (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 ) (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 (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 (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 (20) an1T1 an2T2 annTn Cn trong đó: 4
5. 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  (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): 1 Ti  aij  Ci  (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 (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 (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 (25)  ( x) 2 Để giải (2.25) cần biến đổi như sau : 5
6. a.  T p 1 T p (T p 2T p T p ) (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 (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. Để các nghiệm hội tụ cần điều kiện : (1- 2Fo) 0 (28) tức là : 1 Fo 2 hay phải chọn bước thời gian đủ nhỏ : ( x) 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  (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 )  (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 ) (32) 1 1 2 1 1 - Độ tăng nội năng dU bằng tổng hai dòng nhiệt trên : 6
7. k x h T p T p  (T p T p )  c (T p 1 T p ) (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 (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 (35) 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 (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 (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() (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 (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 : 7
8. 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  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 (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 (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 (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 (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 (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 (45) i, j 2 2 i, j ( x) ( y) c. 8
9. (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 (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 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 q ( ) h(TK Tm ) .I h TK Tm h(TK Tm ) (47) 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. 9
10. + 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 (49) 2 i, j i, j + 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 (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 (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 (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 (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 (54) 10
11. 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 a11T1 a12T2 a1nTn C1 a T a T a T C 21 1 22 2 21 n 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   (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 11
12. D D D Nghiệm sẽ là T 1 ;T 2 ; ;T n ; (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 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 T 1 1 a21 / a21 a22 / a21 a2n / a21 2 C2 / a21   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  12
13. 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) 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 a nn 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 ; 6 6 6 6 hàng chứa T3 có : T3 a3nTn C3 T3 a3nTn C3 ,. 13
14. 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  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 C1 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ó. 14
15. 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  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  (59) k k 0 0 0 1 0 T3 C3 T3 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 15
16. - 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ó : Hình 2.4. Chia mạng tiết diện 1 T T T T T , ngang dầm bêtông i, j 4 i1 i2 i3 i4 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à : 16
17. 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   (60) an1 an2 an3 ann Tn  Cn  Hay ở dạng gọn sau : [aij]. [Ti] = [Ci] (61) Từ đó sẽ rút ra được : - 1 [Ti] = [Ci] [aij] (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 (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 (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 . 17