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

clc

clear

%% 初始參數(shù)

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

c = 1000; ? %比熱容

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

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

Length = 3; %控制體長度3m

T_Left = 3; %左邊界溫度℃

T_Right = 5;%右邊界溫度℃

T_initial = 3; %初始溫度

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

nx = 5; ? ?%x方向網(wǎng)格數(shù)

dx = Length/nx; %x方向距離步長 ?

dt = 100; ?%時間步長100s

maxstep= 200; ?%時間計算步

t_sum = dt*maxstep; ? ? ?%總計算時間s

%% 初始賦值

T0 = zeros(maxstep+1,nx+2);%初始溫度值

T = zeros(maxstep+1,nx+2);%計算溫度值

ae0 = zeros(1,nx+2);

aw0 = zeros(1,nx+2);

ap0 = zeros(1,nx+2);

ap1 = zeros(1,nx+2);

b = zeros(1,nx+2);

%% 初始條件和邊界條件賦值

for t = 1:maxstep+1

? ?for i = 2:nx+1

? ? ? ?T0(t,i) = T_initial;

? ?end

end

T0(:,1) = T_Left;

T0(:,nx+2) = T_Right;

T= T0;

%% 系數(shù)計算

%內(nèi)部網(wǎng)格系數(shù)

for i = 3:nx

? ?ae0(i) = k/dx;%e右，w左

? ?aw0(i) = k/dx;

? ?ap0(i) = rou*c*dx/dt-ae0(i)-aw0(i);

? ?ap1(i)= ae0(i)+aw0(i)+ap0(i);

? ?b(i) = s*dx;

end

%邊界網(wǎng)格系數(shù)

i = 2;

? ?ae0(i) = k/(dx);%e右，w左

? ?aw0(i) = k/(dx/2);

? ?ap0(i) = rou*c*dx/dt-ae0(i)-aw0(i);

? ?ap1(i)= ae0(i)+aw0(i)+ap0(i);

? ?b(i) = s*dx;

i = nx+1;

? ?ae0(i) = k/(dx/2);%e右，w左

? ?aw0(i) = k/(dx);

? ?ap0(i) = rou*c*dx/dt-ae0(i)-aw0(i);

? ?ap1(i)= ae0(i)+aw0(i)+ap0(i);

? ?b(i) = s*dx;

%% 時間步推進(jìn)

for t = 2:maxstep+1

? ?% 顯式計算

? ?for i = 2:nx+1

? ? ? ?T(t,i) = (ae0(1,i)*T0(t-1,i+1) + aw0(1,i)*T0(t-1,i-1) + ap0(1,i)*T0(t-1,i) + b(1,i))/ap1(1,i);

? ?end

? ?T0(t,:) = T(t,:);

end

%% 畫圖

x = zeros(1,nx+2);%距離儲存矩陣矩陣

x(1,1) = 0;

x(1,end) = Length;

x(1,2:end-1) = dx/2:dx:Length-dx/2;

for i=1:maxstep

? ?grid on

? ?ylim([3 6]);

? ?pause (0.005);

? ?a = {['運行時間',num2str(i),'s']};

? ?plot(x,T(i,:));

? ?title(a)

end