用matlab写一个二维温度场模拟程序

% 二维温度场模拟程序
clear all; clc;

%% 参数设置
% 空间网格大小
dx = 0.1; % x方向步长
dy = 0.1; % y方向步长
Lx = 10; % x方向长度
Ly = 10; % y方向长度
nx = Lx/dx; % x方向网格数
ny = Ly/dy; % y方向网格数

% 时间参数
dt = 0.01; % 时间步长
tmax = 10; % 总模拟时间
nt = tmax/dt; % 时间步数

% 物理参数
k = 0.1; % 热传导系数
rho = 1; % 密度
cp = 1; % 热容
alpha = k/(rho*cp); % 热扩散系数

% 初始条件
T0 = 0; % 初始温度

% 边界条件
Tleft = 1; % 左侧温度
Tright = 0; % 右侧温度
Ttop = 1; % 上侧温度
Tbottom = 0; % 下侧温度

%% 初始化
% 初始化温度场
T = ones(ny,nx)*T0; % 全场初始温度
T(1,:) = Tbottom; % 下边界
T(ny,:) = Ttop; % 上边界
T(:,1) = Tleft; % 左边界
T(:,nx) = Tright; % 右边界

% 绘制初始温度场
figure(1)
imagesc(T);
colorbar;
title('Initial Temperature Distribution');

% 初始化温度场数组
Tn = T;

%% 时间迭代
for n = 1:nt
% 内部网格点的温度变化
for i = 2:nx-1
for j = 2:ny-1
Tn(j,i) = T(j,i) + alpha*dt*((T(j,i+1)-2*T(j,i)+T(j,i-1))/dx^2 ...
+ (T(j+1,i)-2*T(j,i)+T(j-1,i))/dy^2);
end
end

% 处理边界网格点的温度变化
Tn(1,:) = Tbottom; % 下边界
Tn(ny,:) = Ttop; % 上边界
Tn(:,1) = Tleft; % 左边界
Tn(:,nx) = Tright; % 右边界

% 更新温度场
T = Tn;

% 绘制温度场
if mod(n,10) == 0 % 每10个时间步绘制一次温度场
figure(1)
imagesc(T);
colorbar;
title(['Temperature Distribution at Time t = ',num2str(n*dt),' s']);
end
end

% 绘制最终温度场
figure(1)
imagesc(T);
colorbar;
title(['Temperature Distribution at Time t = ',num2str(tmax),' s']);