Phương pháp số giải phương trình truyền nhiệt

Ôn Ngũ Minh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br /> 102(02): 155 - 159<br /> <br /> Ngày nhận bài: 03/11/2012, ngày phản biện: 07/3/2013, ngày duyệt đăng: 26/3/2013<br /> <br /> PHƯƠNG PHÁP SỐ GIẢI PHƯƠNG TRÌNH TRUYỀN NHIỆT<br /> Ôn Ngũ Minh*<br /> Trường Đại học Kỹ thuật Công nghiệp – ĐH Thái Nguyên<br /> <br /> TÓM TẮT<br /> Kỹ thuật số giải phương trình đạo hàm riêng gồm hai phương pháp chủ yếu: sai phân hữu hạn và<br /> phần tử hữu hạn. Trong bài viết này, chúng tôi quan tâm tới phương pháp sai phân hữu hạn để giải<br /> phương trình dạng parabolic, sử dụng phương trình truyền nhiệt làm ví dụ. Chúng tôi cũng sử<br /> dụng ngôn ngữ MATLAB để viết mã cho các giải thuật.<br /> Từ khóa: sai phân hữu hạn, phương trình truyền nhiệt, ba đường chéo, phân tích LU.<br /> <br /> Bài toán truyền nhiệt*<br /> Nhiệt độ trong thanh mỏng có thể được mô tả<br /> bởi phương trình truyền nhiệt một chiều:<br /> ut = auxx, 0 < x < X, 0 < t.<br /> Nhiệt độ ban đầu tại mỗi điểm trên thanh<br /> được cho bởi: u(x, 0) = f(x), 0 ≤ x ≤ X.<br /> Ngoài ra, cần có thêm thông tin tại các điểm<br /> mút của thanh. Nếu các mút của thanh được<br /> giữ ở nhiệt độ xác định (ví dụ, ngâm mỗi mút<br /> trong một môi trường chất lỏng với nhiệt độ<br /> dao động theo thời gian), thì các điều kiện<br /> biên là: u(0, t) = g1(t), u(X, t) = g2(t), 0 < t.<br /> Các trạng thái vật lý cũng sẽ xác định dạng<br /> điều kiện biên khác nhau. Ví dụ, giữ một mút<br /> của thanh tương ứng với đạo hàm riêng ux<br /> bằng 0 sẽ là ux(0,t) = 0 hoặc ux(X, t) = 0.<br /> Khả năng thứ ba là một mút của thanh phải<br /> được làm mát (ví dụ, tiếp xúc với không khí<br /> lạnh). Điều kiện biên tương ứng liên quan đến<br /> sự kết hợp của u và ux tại mút thích hợp.<br /> Sự kết hợp của điều kiện ban đầu với các điều<br /> kiện biên nói trên là ví dụ điển hình cho<br /> phương trình đạo hàm riêng dạng parabolic.<br /> <br /> T<br /> *<br /> <br /> Nhiệt độ tại thời điểm t<br /> <br /> Tel: 0913 351286<br /> <br /> Giải phương trình truyền nhiệt<br /> Xét phương trình truyền nhiệt một chiều:<br /> ut = auxx, 0 ≤ x ≤ X, 0 < t ≤ T,<br /> với điều kiện đầu:<br /> u(x, 0) = f(x), 0 ≤ x ≤ X,<br /> và điều kiện biên:<br /> u(0, t) = g1(t), u(X, t) = g2(t), 0 < t ≤ T.<br /> Chia đoạn [0, X] thành n đoạn bằng nhau có<br /> độ dài h = X/n bởi các điểm chia<br /> xi = kh, i = 0, 1, 2, ..., n.<br /> Chia đoạn [0, T] thành m phần bằng nhau có<br /> độ dài k = T/m bởi các điểm chia<br /> tj = jk, j = 0, 1, 2, ..., m.<br /> Với mọi i = 0, 1, 2, ..., n và j = 0, 1, 2, ..., m ta<br /> ký hiệu:<br /> ui,j = u(xi, tj), fi = f(xi),<br /> 155<br /> <br /> x<br /> x=0<br /> <br /> x=X<br /> <br /> Ôn Ngũ Minh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br /> g1j = g1(tj), g2j = g2(tj).<br /> Phương pháp thứ nhất (Explicit)<br /> Tại bước thứ j theo thời gian, ta xấp xỉ ut bởi<br /> công thức sai phân tiến và uxx bởi công thức<br /> sai phân trung tâm:<br /> ut =<br /> uxx =<br /> <br /> 1<br /> ( ui, j +1 − ui, j ) ,<br /> k<br /> <br /> 1<br /> ( ui−1, j − 2ui, j + ui+1, j ) .<br /> h2<br /> <br /> Khi đó<br /> <br /> 1<br /> a<br /> ui , j +1 − ui , j ) = 2 ( ui −1, j − 2ui , j + ui +1, j ) .<br /> (<br /> k<br /> h<br /> ak<br /> Đặt r = 2 ta nhận được:<br /> h<br /> ui,j+1 = rui-1,j + (1 – 2r)ui,j + rui+1,j,<br /> Phương pháp này đơn giản vì cho phép tính<br /> ui,j+1 thông qua ba giá trị của bước trước theo<br /> thời gian. Tuy nhiên, để đảm bảo sự ổn định<br /> của nghiệm, ta cần chọn h và k sao cho<br /> <br /> tj+<br /> t<br /> <br /> xi-<br /> <br /> x<br /> <br /> xi+<br /> <br /> 0 ≤ r ≤ 0.5.<br /> Để cài đặt trên MATLAB, ta tổ chức dữ liệu<br /> vào một ma trận u gồm m+1 hàng và n+1 cột.<br /> Trong đó các giá trị của g1(t) và g2(t) được lưu<br /> tương ứng vào cột thứ nhất và cột cuối cùng<br /> (không kể vị trí đầu tiên trong cột).<br /> Các giá trị của f(x) được lưu tại hàng đầu tiên.<br /> Các giá trị ui,j được lưu ở hàng j+1 cột i+1,<br /> với i = 1, 2, ..., n – 1, j = 1, 2, ..., m.<br /> Như vậy, mỗi giá trị ui,j sẽ được tính dựa vào<br /> ba giá trị của hàng trước mà kề với nó.<br /> <br /> 102(02): 155 - 159<br /> <br /> g1(t) và g2(t). Hàm trả lại ma trận nghiệm u,<br /> véc tơ hàng x và véc tơ hàng t.<br /> function [u,x,t] = heat_explicit(f, g1, g2, X, T, n, m, a)<br /> %solve heat equation u_t = au_xx<br /> % for 0