一维常系数瞬态导热计算_时间显式格式

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

clc

clear

%% 初始参数

rou = 100;  %内能密度

c = 1000;   %比热容

k = 10;     %热流密度系数

s = 10;     %热源强度

Length = 3; %控制体长度3m

T_Left = 3; %左边界温度℃

T_Right = 5;%右边界温度℃

T_initial = 3; %初始温度

%% 划分网格  

nx = 5;    %x方向网格数

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;

%% 系数计算

%内部网格系数

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

%边界网格系数

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;

%% 时间步推进

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