%% 二維常系數瞬態導熱計算，時間顯示格式
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
關鍵詞：
