%% 二维常系数瞬态导热计算，时间显示格式
clc
clear
(资料图片)
%% 初始参数--------------
rou = 100;%内能密度
c = 1000;%比热容
k = 10;%热流密度系数
s = 10;%热源强度
%% 划分网格--------------
%空间网格
LENGTH_x = 3; %总长度3m
nx = 5; %距离网格数
dx = LENGTH_x/nx;%距离步长
LENGTH_y = 2;
ny = 5;
dy = LENGTH_y/ny;
%时间网格
dt = 100;%时间步长100s
nt = 2;%时间网格数，步数
t = dt*nt;%总计算时间
%% 初始赋值--------------
%创建温度储存矩阵矩阵
T = zeros(ny+2,nx+2,nt);
%边界温度℃
T_above = 4;
T_down = 2;
T_Left = 3;
T_Right = 5;
%边界温度赋值
for i=1:nt
T(1,:,i) = T_above;
T(end,:,i) = T_down;
T(:,1,i) = T_Left;
T(:,end,i) = T_Right;
T(1,1,i) = 0;
T(1,end,i) = 0;
T(end,end,i) = 0;
T(end,1,i) = 0;
end
%初始时刻温度赋值，假设为1℃
T(2:(end-1),2:(end-1),1) = 3;
%% 计算参数--------------
for time = 2:nt
% w左、e右、n上、s下
%内部节点
ae0_in = k*dy/dx;
aw0_in = k*dy/dx;
an0_in = k*dx/dy;
as0_in = k*dx/dy;
ap0_in = rou*c*dx*dy/dt-ae0_in-aw0_in-an0_in-as0_in;
ap1_in = ae0_in+aw0_in+an0_in+as0_in+ap0_in;
b = s*dx*dy;
for i = 3:ny
for j = 3:nx
T(i,j,time) = (ae0_in*T(i,j+1,time-1) + aw0_in*T(i,j-1,time-1) + ...
an0_in*T(i-1,j,time-1) + as0_in*T(i+1,j,time-1) + ...
ap0_in*T(i,j,time-1) + b)/ap1_in;
end
end
%左侧单边界节点
ae0_outl_1 = k*dy/dx;
aw0_outl_1 = k*dy/(dx/2);
an0_outl_1 = k*dx/dy;
as0_outl_1 = k*dx/dy;
ap0_outl_1 = rou*c*dx*dy/dt-ae0_outl_1-aw0_outl_1-an0_outl_1-as0_outl_1;
ap1_outl_1 = ae0_outl_1+aw0_outl_1+an0_outl_1+as0_outl_1+ap0_outl_1;
b = s*dx*dy;
for i = 3:ny
T(i,2,time) = (ae0_outl_1*T(i,3,time-1) + aw0_outl_1*T(i,1,time-1) + ...
an0_outl_1*T(i-1,2,time-1) + as0_outl_1*T(i+1,2,time-1) + ...
ap0_outl_1*T(i,2,time-1) + b)/ap1_outl_1;
end
%右侧单边界节点
ae0_outR_1 = k*dy/(dx/2);
aw0_outR_1 = k*dy/dx;
an0_outR_1 = k*dx/dy;
as0_outR_1 = k*dx/dy;
ap0_outR_1 = rou*c*dx*dy/dt-ae0_outR_1-aw0_outR_1-an0_outR_1-as0_outR_1;
ap1_outR_1 = ae0_outR_1+aw0_outR_1+an0_outR_1+as0_outR_1+ap0_outR_1;
b = s*dx*dy;
for i = 3:ny
T(i,end-1,time) = (ae0_outR_1*T(i,end,time-1) + aw0_outR_1*T(i,end-2,time-1) + ...
an0_outR_1*T(i-1,end-1,time-1) + as0_outR_1*T(i+1,end-1,time-1) + ...
ap0_outR_1*T(i,end-1,time-1) + b)/ap1_outR_1;
end
%上侧单边界节点
ae0_outa_1 = k*dy/dx;
aw0_outa_1 = k*dy/dx;
an0_outa_1 = k*dx/(dy/2);
as0_outa_1 = k*dx/dy;
ap0_outa_1 = rou*c*dx*dy/dt-ae0_outa_1-aw0_outa_1-an0_outa_1-as0_outa_1;
ap1_outa_1 = ae0_outa_1+aw0_outa_1+an0_outa_1+as0_outa_1+ap0_outa_1;
b = s*dx*dy;
for j = 3:nx
T(2,j,time) = (ae0_outa_1*T(2,j+1,time-1) + aw0_outa_1*T(2,j-1,time-1) + ...
an0_outa_1*T(1,j,time-1) + as0_outa_1*T(3,j,time-1) + ...
ap0_outa_1*T(2,j,time-1) + b)/ap1_outa_1;
end
%下侧单边界节点
ae0_outd_1 = k*dy/dx;
aw0_outd_1 = k*dy/dx;
an0_outd_1 = k*dx/dy;
as0_outd_1 = k*dx/(dy/2);
ap0_outd_1 = rou*c*dx*dy/dt-ae0_outd_1-aw0_outd_1-an0_outd_1-as0_outd_1;
ap1_outd_1 = ae0_outd_1+aw0_outd_1+an0_outd_1+as0_outd_1+ap0_outd_1;
b = s*dx*dy;
for j = 3:nx
T(end-1,j,time) = (ae0_outd_1*T(end-1,j+1,time-1) + aw0_outd_1*T(end-1,j-1,time-1) + ...
an0_outd_1*T(end-2,j,time-1) + as0_outd_1*T(end,j,time-1) + ...
ap0_outd_1*T(end-1,j,time-1) + b)/ap1_outd_1;
end
%左上角两边界节点
ae0_outwn_2 = k*dy/(dx);
aw0_outwn_2 = k*dy/(dx/2);
an0_outwn_2 = k*dx/(dy/2);
as0_outwn_2 = k*dx/(dy);
ap0_outwn_2 = rou*c*dx*dy/dt-ae0_outwn_2-aw0_outwn_2-an0_outwn_2-as0_outwn_2;
ap1_outwn_2 = ae0_outwn_2+aw0_outwn_2+an0_outwn_2+as0_outwn_2+ap0_outwn_2;
b = s*dx*dy;
T(2,2,time) = (ae0_outwn_2*T(2,3,time-1) + aw0_outwn_2*T(2,1,time-1) + ...
an0_outwn_2*T(1,2,time-1) + as0_outwn_2*T(3,2,time-1) + ...
ap0_outwn_2*T(2,2,time-1) + b)/ap1_outwn_2;
%左下角两边界节点
ae0_outws_2 = k*dy/(dx);
aw0_outws_2 = k*dy/(dx/2);
an0_outws_2 = k*dx/(dy);
as0_outws_2 = k*dx/(dy/2);
ap0_outws_2 = rou*c*dx*dy/dt-ae0_outws_2-aw0_outws_2-an0_outws_2-as0_outws_2;
ap1_outws_2 = ae0_outws_2+aw0_outws_2+an0_outws_2+as0_outws_2+ap0_outws_2;
b = s*dx*dy;
T(end-1,2,time) = (ae0_outws_2*T(end-1,3,time-1) + aw0_outws_2*T(end-1,1,time-1) + ...
an0_outws_2*T(end-2,2,time-1) + as0_outws_2*T(end,2,time-1) + ...
ap0_outws_2*T(end-1,2,time-1) + b)/ap1_outws_2;
%右上角两边界节点
ae0_outne_2 = k*dy/(dx/2);
aw0_outne_2 = k*dy/(dx);
an0_outne_2 = k*dx/(dy/2);
as0_outne_2 = k*dx/(dy);
ap0_outne_2 = rou*c*dx*dy/dt-ae0_outne_2-aw0_outne_2-an0_outne_2-as0_outne_2;
ap1_outne_2 = ae0_outne_2+aw0_outne_2+an0_outne_2+as0_outne_2+ap0_outne_2;
b = s*dx*dy;
T(2,end-1,time) = (ae0_outne_2*T(2,end,time-1) + aw0_outne_2*T(2,end-2,time-1) + ...
an0_outne_2*T(1,end-1,time-1) + as0_outne_2*T(3,end-1,time-1) + ...
ap0_outne_2*T(2,end-1,time-1) + b)/ap1_outne_2;
%右下角两边界节点
ae0_outes_2 = k*dy/(dx/2);
aw0_outes_2 = k*dy/(dx);
an0_outes_2 = k*dx/(dy);
as0_outes_2 = k*dx/(dy/2);
ap0_outes_2 = rou*c*dx*dy/dt-ae0_outes_2-aw0_outes_2-an0_outes_2-as0_outes_2;
ap1_outes_2 = ae0_outes_2+aw0_outes_2+an0_outes_2+as0_outes_2+ap0_outes_2;
b = s*dx*dy;
T(end-1,end-1,time) = (ae0_outes_2*T(end-1,end,time-1) + aw0_outes_2*T(end-1,end-2,time-1) + ...
an0_outes_2*T(end-2,end-1,time-1) + as0_outes_2*T(end,end-1,time-1) + ...
ap0_outes_2*T(end-1,end-1,time-1) + b)/ap1_outes_2;
end
关键词：
