%% 二维常系数瞬态导热计算，时间显示格式

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