%% 二維常系數(shù)瞬態(tài)導(dǎo)熱計(jì)算，時(shí)間顯示格式

clc

clear

%% 初始參數(shù)--------------

rou = 100;%內(nèi)能密度

c = 1000;%比熱容

k = 10;%熱流密度系數(shù)

s = 10;%熱源強(qiáng)度

%% 劃分網(wǎng)格--------------

%空間網(wǎng)格

LENGTH_x = 3; ? %總長度3m

nx = 5; ? ? ? ? %距離網(wǎng)格數(shù)

dx = LENGTH_x/nx;%距離步長

LENGTH_y = 2;

ny = 5;

dy = LENGTH_y/ny;

%時(shí)間網(wǎng)格

dt = 100;%時(shí)間步長100s

nt = 2;%時(shí)間網(wǎng)格數(shù)，步數(shù)

t = dt*nt;%總計(jì)算時(shí)間

%% 初始賦值--------------

%創(chuàng)建溫度儲存矩陣矩陣

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

%初始時(shí)刻溫度賦值，假設(shè)為1℃

T(2:(end-1),2:(end-1),1) = 3;

%% 計(jì)算參數(shù)--------------

for time = 2:nt

? ?% w左、e右、n上、s下

? ?%內(nèi)部節(jié)點(diǎn)

? ?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

? ?%左側(cè)單邊界節(jié)點(diǎn)

? ?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

? ?%右側(cè)單邊界節(jié)點(diǎn)

? ?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

? ?%上側(cè)單邊界節(jié)點(diǎn)

? ?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

? ?%下側(cè)單邊界節(jié)點(diǎn)

? ?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

? ?%左上角兩邊界節(jié)點(diǎn)

? ?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;

? ?%左下角兩邊界節(jié)點(diǎn)

? ?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;

? ?%右上角兩邊界節(jié)點(diǎn)

? ?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;

? ?%右下角兩邊界節(jié)點(diǎn)

? ?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