【编程训练】二维导热 - 差分模拟

请注意,文章中部分理解可能已被笔者废弃,本文最后更新于:1 个月前

显示迭代法

问题

模拟绝热环境下有初始温度差的二维平面的温度随时间的变化情况。

(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

模型说明

二维热传导方程为

Tt=k(2Tx2+2Ty2)\frac{\partial T}{\partial t}=k\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}\right)

显式有限差分格式

Ti,jn+1Ti,jnΔt=(Ti+1,jn2Ti,jn+Ti1,jnΔx2+Ti,j+1n2Ti,jn+Ti,j1nΔy2)\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)

显式迭代格式免去了矩阵求解的麻烦,但却是条件稳定的。即只有保证时间步长和空间步长满足一定关系时才能获得收敛解,对于热传导问题需要满足:

Δtmin(Δx2,Δy2)12\frac{\Delta t}{\min(\Delta x^2,\Delta y^2)}\le \frac{1}{2}

隐式有限差分格式(BTCS)

Ti,jn+1Ti,jnΔt=(Ti+1,jn+12Ti,jn+1+Ti1,jn+1Δx2+Ti,j+1n+12Ti,jn+1+Ti,j1n+1Δy2)\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)

Ti,jn+1Ti,jnΔt=0.5(Ti+1,jn+12Ti,jn+1+Ti1,jn+1Δx2+Ti,j+1n+12Ti,jn+1+Ti,j1n+1Δy2)+0.5(Ti+1,jn2Ti,jn+Ti1,jnΔx2+Ti,j+1n2Ti,jn+Ti,j1nΔy2)\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)

  1. x方向隐式更新

Ti,jn+1Ti,jnΔt=k(Ti+1,jn+12Ti,jn+1+Ti1,jn+1Δx2+Ti,j+1n2Ti,jn+Ti,j1nΔy2)\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)

  1. y方向隐式更新

Ti,jn+1Ti,jnΔt/2=k(Ti+1,jn2Ti,jn+Ti1,jnΔx2+Ti,j+1n+12Ti,jn+1+Ti,j1n+1Δy2)\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)

交替隐式迭代法是无条件稳定的,且时间和空间离散精度均为二阶,同时矩阵求解计算量优于上述其他隐式算法,具有独特的优势。

参考