【编程训练】二维导热 - 差分模拟
显示迭代法
问题
模拟绝热环境下有初始温度差的二维平面的温度随时间的变化情况。
（gif见知乎文章）
代码
%二维热传导过程显式差分模拟
clear;close all;clc;
t = 0.03; %时间范围，计算到0.03h
x = 1;y = 1; %空间范围，0-1m
m = 320; %时间t方向分320个格子
n = 32; %空间x方向分32个格子
k = 32; %空间y方向分32个格子
ht = t/(m-1); %时间步长dt
hx = x/(n-1); %空间步长dx
hy = y/(k-1); %空间步长dy
u = zeros(m,n,k);
%设置边界
[x,y] = meshgrid(0:hx:1,0:hy:1);
u(1,:,:) = sin(4*pi*x)+cos(4*pi*y);
%按照公式进行差分
for ii=1:m-1
for jj=2:n-1
for kk=2:k-1
u(ii+1,jj,kk) = ht*(u(ii,jj+1,kk)+u(ii,jj-1,kk)-2*u(ii,jj,kk))/hx^2 + ...
ht*(u(ii,jj,kk+1)+u(ii,jj,kk-1)-2*u(ii,jj,kk))/hy^2 + u(ii,jj,kk);
end
end
end
for i=1:200
figure(1);
mesh(x,y,reshape(u(i,:,:),[n k]));
axis([0 1 0 1 -2 2]);
F=getframe(gcf); %捕获显示在屏幕上的当前坐标区作为影片帧
I=frame2im(F); %将捕获的影片帧转换为图像数据
[I,map]=rgb2ind(I,256);
if i == 1
imwrite(I,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.05);
else
imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.05);
end
end交替隐式迭代法
问题
模拟外界温度周期变化的影响下，二维地板的温度随时间的变化情况。
（gif见知乎文章）
代码
%二维热传导过程差分交替隐式迭代法(Alternating-direction implicit method)
%c过大会导致交替迭代非常不稳定
clear;close all;clc;
%% 参数定义
t = 72; %时间范围
xl = 20;yl = 20; %空间范围
ht = 0.5; %时间步长dt
nt = t/ht; %时间t方向分格子
nx = 50; %空间x方向分格子
ny = 50; %空间y方向分3格子
hx = xl/(nx-1); %空间步长dx
hy = yl/(ny-1); %空间步长dy
k = 1.93; %二维材料导热系数
c = 1e-1; %与空气换热参数
u = zeros(nt,ny,nx);
sqt = @(x) -12.5*cos(3.1415926535*x/12)+72.5; %温度计算
%设置边界与初始温度
[x,y] = meshgrid(-hx:hx:xl+hx,-hy:hy:yl+hy);
u = 72.5*ones(size(x));
u(1,:)=sqt(ht)*ones(size(u(1,:)));
u(end,:)=sqt(ht)*ones(size(u(end,:)));
u(:,1)=sqt(ht)*ones(size(u(:,1)));
u(:,end)=sqt(ht)*ones(size(u(:,end)));
log(1,:,:)=u;
%% 构造迭代矩阵 A,B
A = eye(ny)*(-2);
A(2:end,1:end-1) = A(2:end,1:end-1)+eye(ny-1);
A(1:end-1,2:end) = A(1:end-1,2:end)+eye(ny-1);
A = A*(k*ht/hy^2);
%A = inv(eye(ny)-A);
tmp = zeros(size(A)+[2,2]);
tmp(2:end-1,2:end-1) = A;
tmp(2,1:3) = [ht*c/2,-ht*c/2,0];
tmp(end-1,end-2:end) = [0,-ht*c/2,ht*c/2];
%tmp(1,1)=1;tmp(end,end)=1;
A = tmp;
B = eye(nx)*(-2);
B(2:end,1:end-1) = B(2:end,1:end-1)+eye(nx-1);
B(1:end-1,2:end) = B(1:end-1,2:end)+eye(nx-1);
B = B*(k*ht/hx^2);
%B = eye(nx)+B;
tmp = zeros(size(B)+[2,2]);
tmp(2:end-1,2:end-1)=B;
tmp(1:3,2) = [ht*c/2,-ht*c/2,0]';;
tmp(end-2:end,end-1) = [0,-ht*c/2,ht*c/2]';
%tmp(1,1)=1;tmp(end,end)=1;
B = tmp;
%% 差分模拟
for i = 2:nt
%维护空气温度
u(1,:)=sqt(i*ht)*ones(size(u(1,:)));
u(end,:)=sqt(i*ht)*ones(size(u(end,:)));
u(:,1)=sqt(i*ht)*ones(size(u(:,1)));
u(:,end)=sqt(i*ht)*ones(size(u(:,end)));
%交替隐式迭代法
if mod(i,2)==0
u=(eye(size(A))-A)\u*(eye(size(B))+B);
else
u=(eye(size(A))+A)*u/(eye(size(B))-B);
end
log(end+1,:,:)=u;
end
%计算平均温度
for i=1:nt
logMeanT(i)=mean(mean(mean(log(i,2:end-1,2:end-1))));
end
figure(1);
plot(ht:ht:t,logMeanT);
fprintf('极差:%d max:%d min:%d\n',max(logMeanT(48/ht:72/ht))-min(logMeanT(48/ht:72/ht)),max(logMeanT(48/ht:72/ht)),min(logMeanT(48/ht:72/ht)));
%% 可视化过程
for i=1:nt
figure(2);
mesh(x,y,reshape(log(i,:,:),size(u)));
xlabel('x');
title(num2str(i*ht));
axis([0-hx,xl+hx,0-hy,yl+hy 60 85]);
pause(0.1);
F=getframe(gcf); %捕获显示在屏幕上的当前坐标区作为影片帧
I=frame2im(F); %将捕获的影片帧转换为图像数据
[I,map]=rgb2ind(I,256);
if i == 1
imwrite(I,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.05);
else
imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.05);
end
end模型说明
二维热传导方程为
\[\frac{\partial T}{\partial t}=k\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}\right)\]
显式有限差分格式
\[\frac{T_{i, j}^{n+1}-T_{i, j}^{n}}{\Delta t}=\left(\frac{T_{i+1, j}^{n}-2 T_{i, j}^{n}+T_{i-1, j}^{n}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n}-2 T_{i, j}^{n}+T_{i, j-1}^{n}}{\Delta y^{2}}\right)\]
显式迭代格式免去了矩阵求解的麻烦，但却是条件稳定的。即只有保证时间步长和空间步长满足一定关系时才能获得收敛解，对于热传导问题需要满足：
\[\frac{\Delta t}{\min(\Delta x^2,\Delta y^2)}\le \frac{1}{2}\]
隐式有限差分格式(BTCS)
\[\frac{T_{i, j}^{n+1}-T_{i, j}^{n}}{\Delta t}=\left(\frac{T_{i+1, j}^{n+1}-2 T_{i, j}^{n+1}+T_{i-1, j}^{n+1}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n+1}-2 T_{i, j}^{n+1}+T_{i, j-1}^{n+1}}{\Delta y^{2}}\right)\]
需要注意BTCS法是无条件稳定的，即时间步长可以任意取都不会计算发散（不发散不代表算得准）。由于BTCS时间项采用向后差分，因此该方法在时间层精度为一阶，而空间层精度为二阶。
隐式有限差分格式(Crank–Nicolson)
\[\frac{T_{i, j}^{n+1}-T_{i, j}^{n}}{\Delta t}=0.5\left(\frac{T_{i+1, j}^{n+1}-2 T_{i, j}^{n+1}+T_{i-1, j}^{n+1}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n+1}-2 T_{i, j}^{n+1}+T_{i, j-1}^{n+1}}{\Delta y^{2}}\right) \\ +0.5\left(\frac{T_{i+1, j}^{n}-2 T_{i, j}^{n}+T_{i-1, j}^{n}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n}-2 T_{i, j}^{n}+T_{i, j-1}^{n}}{\Delta y^{2}}\right)\]
作为隐式方法，Crank–Nicolson格式是无条件稳定的，且其时间离散精度为二阶。
交替隐式迭代法(Alternating-direction implicit method)
- x方向隐式更新
\[\frac{T_{i, j}^{n+1}-T_{i, j}^{n}}{\Delta t}=k\left(\frac{T_{i+1, j}^{n+1}-2 T_{i, j}^{n+1}+T_{i-1, j}^{n+1}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n}-2 T_{i, j}^{n}+T_{i, j-1}^{n}}{\Delta y^{2}}\right)\]
- y方向隐式更新
\[\frac{T_{i, j}^{n+1}-T_{i, j}^{n}}{\Delta t / 2}=k\left(\frac{T_{i+1, j}^{n}-2 T_{i, j}^{n}+T_{i-1, j}^{n}}{\Delta x^{2}}+\frac{T_{i, j+1}^{n+1}-2 T_{i, j}^{n+1}+T_{i, j-1}^{n+1}}{\Delta y^{2}}\right)\]
交替隐式迭代法是无条件稳定的，且时间和空间离散精度均为二阶，同时矩阵求解计算量优于上述其他隐式算法，具有独特的优势。
参考
- matlab练习程序（差分法解二维热传导方程）. http://t.zoukankan.com/tiandsp-p-14406223.html
- 一文看懂交替隐式迭代法. https://zhuanlan.zhihu.com/p/522899015