Xem mẫu
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
GIẢI SỐ PHƢƠNG TRÌNH TRUYỀN NHIỆT 2D
Nguyễn Thị Trà My, Lớp K60E, Khoa Toán – Tin
GVHD: TS. Nguyễn Hùng Chính
Tóm tắt: Mô phỏng toán học là một ngành đã và đang phát triển hết sức mạnh mẽ trên thế giới và có
vai trò quan trọng trong hầu hết các lĩnh vực của đời sống xã hội. Cùng với sự phát triển của công
nghệ thông tin nói chung và công nghệ tính toán nói riêng, những mô hình toán học phức tạp (xuất
phát từ các khoa học và thực tiễn) đã được số hóa thành công nghệ và đem lại hiệu quả kinh tế cao.
Báo cáo này trình bày nghiên cứu về phương pháp giải số để giải phương trình truyền nhiệt hai chiều,
một trong những phương trình toán học có nhiều ứng dụng, và đề xuất kĩ thuật đưa thuật toán vào
máy tính để xây dựng chương trình mô phỏng quá trình truyền nhiệt theo thời gian.
Từ khóa: Phương pháp xấp xỉ sai phân, lược đồ tường minh, sự ổn định của lược đồ, phương trình
truyền nhiệt, điều kiện biên Dirichlet.
I. MỞ ĐẦU
Tại sao phải giải số (xấp xỉ nghiệm) của một phƣơng trình toán học? Nhƣ chúng ta
đã biết, phần lớn các mô hình toán học trong thực tế đều không giải đƣợc nghiệm đúng, vì
vậy cần xấp xỉ nghiệm và điều khiển đƣợc sai số của nghiệm gần đúng. Hơn nữa, việc giải
số sẽ đƣa đến thuật toán, tức là ta có thể ra lệnh cho máy tính thực hiện các phép tính để
tìm ra kết quả.
Phƣơng trình truyền nhiệt là một phƣơng trình đạo hàm riêng quan trọng xuất phát từ
mô hình vật lí và có giá trị thực tiễn nhất định. Việc giải phƣơng trình truyền nhiệt cho ta
khảo sát sự phân bố nhiệt lƣợng theo thời gian của một vùng chất điểm, nhƣ một thanh kim
loại (với trƣờng hợp 1 chiều) và mảnh kim loại (với trƣờng hợp 2 chiều).
Phƣơng pháp phổ biến để giải đúng phƣơng trình truyền nhiệt vẫn đƣợc biết đến là
phƣơng pháp Fourier (đƣợc phát triển từ năm 1822 bởi nhà toán học Joseph Fourier). Tuy
nhiên, trong thực hành thí nghiệm, việc cần thay đổi các dữ kiện bài toán, cũng nhƣ việc số
liệu trong thực hành phải tính toán lớn đều làm cho việc giải đúng gặp khó khăn, vì vậy mà
ngƣời ta cần đến phƣơng pháp số để giải phƣơng trình truyền nhiệt.
Phƣơng pháp số là một lĩnh vực của toán học chuyên nghiên cứu các phƣơng pháp
giải gần đúng các bài toán dựa trên những số liệu cụ thể và cho kết quả dƣới dạng số. Với
sự hỗ trợ của máy tính, phƣơng pháp số là công cụ không thể thiếu cho phép thực hiện tính
toán với tốc độ tính toán nhanh và khối lƣợng tính toán lớn.
Báo cáo tập trung vào việc xây dựng phƣơng pháp giải số, thuật toán để giải bài toán
truyền nhiệt, sự ổn định, điều kiện ổn định và đề xuất kĩ thuật số hóa lƣợc đồ, lập chƣơng
trình máy tính và mô phỏng số trong MATLAB.
II. NỘI DUNG
1. Phƣơng pháp số và thuật toán
1.1. Bài toán
Xét phƣơng trình truyền nhiệt trên miền Ω = [a,b]×[c,d] với nguồn nhiệt f(x,y,t) thay
đổi theo thời gian:
7
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
với là một hằng số dƣơng đặc trƣng cho vận tốc truyền nhiệt.
Giả sử trạng thái ban đầu là:
và điều kiện trên biên xác định bởi:
,
1.2. Giải số (xấp xỉ nghiệm) phương trình
Ta phân hoạch miền Ω bởi lƣới hình chữ nhật (xem Hình 1) có các đỉnh xác định bởi
điểm với i = 1,2,...,Nx và j = 1,2,...,Ny trong đó:
Đồng thời, ta xét phân hoạch thời gian t bởi các điểm , và kí
hiệu:
Để giải gần đúng phƣơng trình (1), trƣớc tiên ta xấp xỉ các đạo hàm [1] bậc nhất theo
biến thời gian và đạo hàm bậc hai theo biến không gian bởi:
Áp dụng vào phƣơng trình (1), ta đƣợc:
Biểu diễn theo các phần tử còn lại, ta có lƣợc đồ:
trong đó:
8
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
1.3. Điều kiện ổn định của lược đồ
Trong phần này, chúng ta nghiên cứu điều kiện ổn định của lƣợc đồ (2) trong trƣờng
hợp không có nguồn nhiệt tác động vào hệ và điều kiện biên đồng nhất bằng không.
Nghiệm gần đúng của (1) là bộ các giá trị rời rạc có dạng . Ta
thực hiện phép biến đổi Fourier hai chiều cho nghiệm xấp xỉ trên nhƣ sau:
Ta có tính chất sau đây: Giả sử với , ta có:
Nhƣ vậy:
Nếu với , thì ,
Tƣơng tự, ta có:
- Nếu với , thì ,
- Nếu với , thì ,
- Nếu với , thì .
Áp dụng tính chất trên đây của phép biến đổi Fourier vào lƣợc đồ (2) với , ta
đƣợc:
Vì nên ta suy ra:
với
Lƣợc đồ (2) đƣợc gọi là ổn định khi và chỉ khi
tức là:
hay:
9
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
Từ đó suy ra :
Nhƣ vậy, ràng buộc (*) chính là điều kiện ổn định của lƣợc đồ. Trong quá trình mô
phỏng, các dữ liệu đầu vào cần phải thỏa mãn điều kiện (*) để đảm bảo thuật toán trong
máy tính hội tụ và cho kết quả tin cậy.
1.4. Dạng ma trận của lược đồ
Điều đầu tiên, ta nhận thấy với cách biểu diễn lƣợc đồ ở trên, ta đã phân hoạch miền
Ω thành lƣới chữ nhật tƣơng ứng với ma trận vuông cấp , và mỗi điểm lƣới ui,j trong
lƣợc đồ đều phải đƣợc tính từ 4 số hạng lân cận là ui−1,j (bên trái), ui+1,j (bên phải), ui,j−1
(bên trên), ui,j+1 (bên dưới):
Hình 1. Phân hoạch, sơ đồ các điểm trên lưới ảnh hướng đến việc tính
Xét lƣợc đồ (2) với , ta có:
Cho chạy, tức là viết phƣơng trình trên lần lƣợt với và biểu diễn kết
quả thu đƣợc dƣới dạng ma trận, ta có:
10
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
(3)
Mặt khác, quan sát ma trận , ta nhận thấy có 2 thành phần ,
không đƣợc lặp lại trong quá trình tính toán, đó cũng chính là giá trị hàm nhiệt độ xác định
ở thời điểm tm, tức là α(tm), β(tm), ở điểm biên đầu x0 và điểm biên cuối trên đƣờng
thẳng phân hoạch tại . Vì vậy, khi chuyển hệ trên thành dạng ma trận ta cần đƣa hai
thành phần này ra ngoài thành số hạng tự do (vectơ b).
Tiếp theo, ta đặt:
Khi đó, (3) đƣợc viết gọn nhƣ sau:
trong đó, I là ma trận vuông đơn vị cấp , B là ma trận vuông cấp xác định bởi:
Thực hiện tƣơng tự nhƣ trên với , ta thu đƣợc một hệ gồm
phƣơng trình tuyến tính với các biến là các vectơ chiều:
Tƣơng tự nhƣ đã xử lí với ma trận B, ta chuyển 2 phần tử mang giá trị ở biên trên và
dƣới ra một vectơ tự do. Nhƣ vậy dạng ma trận của lƣợc đồ (3) là:
11
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
.
trong đó, ma trận A là một ma trận vuông đối xứng có cấp là .
Ta xác định cụ thể tại những vị trí 1, 2, 3,.., Nx chứa giá trị của biên dƣới (Γ3), vị trí
Nx, 2Nx, 3Nx,..., NxNy chứa giá trị của biên phải (Γ2) 1, Nx + 1, 2Nx+1 ,..., (Ny−1)Nx+1 chứa
giá trị của biên trái (Γ1), NxNy−Nx, ..., NxNy chứa giá trị của biên trên (Γ4).
Việc xác định này sẽ rất thuận lợi trong việc gán các giá trị biên vào đúng vị trí trong
các vectơ điều kiện biên tùy theo yêu cầu của bài toán.
2. Mã lệnh trong MATLAB và mô phỏng nghiệm của bài toán
2.1. Sơ đồ hóa lược đồ ta cần thực hiện các bước
Bảng 1. Đoạn code của ma trận A và vectơ điều kiện biên trong bài toán
% Dieu kien bien % Xay dung ma tran A
Left = 0.0d0 * ones(Ny,1); % Duong cheo chinh
Right = 80.d0 * ones(Ny,1); e = ones(Nx,1);
Top = 20.d0 * ones(1,Nx); V0 = -2*mu*(LambdaX + LambdaY)*e;
Bottom = 20.d0 * ones(1,Nx); V1 = mu*LambdaY * e;
% Dieu kien bien trai + phai T1 = spdiags([V1 V0 V1],-1:1,Nx,Nx);
G = zeros(Nx*Ny,1); I1 = eye(Nx);
G(1:Nx:Nx*Ny) = mu* LambdaY * T1 = T1 + I1;
Left; A1 = kron(eye(Ny),T1);
G(Nx:Nx:Nx*Ny) = mu * LambdaY * % Duong cheo tren, duoi
Right; I2 = mu * LambdaX * eye(Nx);
% Dieu kien bien tren + duoi T2 = spdiags([ones(Ny,1)],1,Ny,Ny);
H = zeros(Nx*Ny,1); A2 = kron(T2,I2);
H(1:Nx) = mu * LambdaX * Bottom; A3 = kron(T2',I2);
H(Nx*Ny - Nx + 1:Nx*Ny) = mu * A = A1 + A2 + A3;
LambdaX * Top;
- Mở một tệp tin "heatequation2D.m"
- Khai báo các dữ kiện của bài toán.
12
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
- Khai báo các bƣớc phân hoạch theo thời gian và không gian.
- Xây dựng vectơ điều kiện biên b.
- Xây dựng vectơ hàm F
- Xây dựng ma trận A.
- Thực hiện thuật toán để giải.
- Vẽ đồ thị nghiệm của bài toán tại thời điểm muốn khảo sát.
Điều khó khăn nhất của ta trong việc sơ đồ hóa lƣợc đồ để thực hiện thuật toán trong
MATLAB chính là viết thuật toán xây dựng ma trận A và vectơ điều kiện biên b. Bảng 1
mô tả đoạn code của ma trận A và vectơ điều kiện biên trong bài toán.
2.2. Sử dụng chương trình máy tính để mô phỏng quá trình truyền nhiệt
Dựa vào kết quả mô phỏng, ta có thể đánh giá về tính đúng đắn của lƣợc đồ mà
chúng ta xây dựng. Hơn nữa, bằng mô phỏng và chuẩn hóa chƣơng trình tính toán trên
MATLAB, một phƣơng trình đạo hàm riêng thuần túy toán học có thể mang lại hiệu quả
thực tiễn nhất định khi giải quyết đƣợc trong nhiều trƣờng hợp, nhanh chóng và dễ thực
hiện hơn nhiều so với việc giải đúng.
Sau đây, tôi xin giới thiệu một số kết quả mô phỏng từ các bài toán có điều kiện ban
biên khác nhau:
a) Ví dụ 1: Bài toán truyền nhiệt 2 chiều với điều kiện ban đầu:
u0= 20 và điều kiện biên: α0 = 0, α1 = 80, β0 = 20, β1 = 20.
Hình 2. Sự truyền nhiệt của mảnh kim loại (5x5) có hệ số C=1, tại thời điểm t = 5
b) Ví dụ 2: Bài toán truyền nhiệt 2 chiều với điều kiện ban đầu: u0= 20 và điều kiện
biên: α0 = 20, α1 = 20, β0 = 200, β1 = 200.
13
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
Hình 3. Sự truyền nhiệt của mảnh kim loại (5x5) có hệ số C=1 qua các thời điểm t’
c) Ví dụ 3: Bài toán truyền nhiệt 2 chiều với điều kiện ban đầu: u0= 20 và điều kiện
biên: α0 = 200, α1 = 200, β0 = 200, β1 = 200.
Hình 4. Sự truyền nhiệt của mảnh kim loại (5x5) có hệ số C=1 qua các thời điểm t”
III. KẾT LUẬN
Đề tài đã tập trung vào những công việc sau:
1. Đề xuất một phƣơng án giải số phƣơng trình truyền nhiệt 2 chiều.
2. Nghiên cứu lí thuyết về sự ổn định của lƣợc đồ và kiểm chứng, so sánh kết quả
thông qua mô phỏng.
14
- KỈ YẾU HỘI NGHỊ SINH VIÊN NGHIÊN CỨU KHOA HỌC NĂM HỌC 2013-2014
3. Xây dựng đƣợc chƣơng trình tính toán, cho phép mô phỏng nhiều bài toán cụ thể
từ thực tế, có thể tái sử dụng dễ dàng vào các nghiên cứu khác.
Một số hƣớng nghiên cứu tiếp theo:
1. Giải số phƣơng trình truyền nhiệt với vận tốc truyền nhiệt khác hằng số (mô hình
sát thực tế là tổng quát hơn).
2. Xây dựng lƣợc đồ giải gần đúng cho bài toán trên một miền xác định không phải
là hình chữ nhật. Đây là bài toán khó hơn rất nhiều, có thể đòi hỏi những kiến thức toán
học khác.
3. Song song hóa thuật toán phục vụ cho những ứng dụng đòi hỏi khối lƣợng tính
toán lớn.
TÀI LIỆU THAM KHẢO
[1] Nguyễn Minh Chƣơng, Nguyễn Văn Khải, Khuất Văn Ninh, Nguyễn Văn Tuấn,
Nguyễn Tƣờng, Giải tích số, NXB Giáo dục, 2007.
[2] Laurent Di Menza, Analyse numérique des équations aux derivée partielles,
Published by Cassini, 2009.
[3] Peter J. Olver, Numerical Methods, Lecture Notes, 2011.
15
nguon tai.lieu . vn
