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

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