1维非稳态导热问题求解的Matlab实现

%按照@计算传热学大叔的视频指导写了一个小程序，欢迎各位指正。

clc;

clear;

tic

format longG

rho=100;

cp=1000;

k=10;

T0=3;

TW=3;

TE=5;

s=10;

t=60000;

x=3;

dt=1;

dx=0.015;

sizet=t/dt;

sizex=x/dx;

Tmatrix=zeros(sizet,sizex);

for i=1

   Tmatrix(i,:)=T0;

end

%先用显格式算一算

for i=2:sizet

   for j=1:1

       ap1=rho*cp*dx/dt;

       aw0=2*k/dx;

       ae0=k/dx;

       ap0=ap1-aw0-ae0;

       Tmatrix(i,j)=ap0/ap1*Tmatrix(i-1,j)+aw0/ap1*TW+ae0/ap1*Tmatrix(i-1,j+1)+s*dx/ap1;

   end

   for j=2:sizex-1

       ap1=rho*cp*dx/dt;

       aw0=k/dx;

       ae0=k/dx;

       ap0=ap1-aw0-ae0;

       Tmatrix(i,j)=ap0/ap1*Tmatrix(i-1,j)+aw0/ap1*Tmatrix(i-1,j-1)+ae0/ap1*Tmatrix(i-1,j+1)+s*dx/ap1;

   end

   for j=sizex:sizex

       ap1=rho*cp*dx/dt;

       aw0=k/dx;

       ae0=2*k/dx;

       ap0=ap1-aw0-ae0;

       Tmatrix(i,j)=ap0/ap1*Tmatrix(i-1,j)+aw0/ap1*Tmatrix(i-1,j-1)+ae0/ap1*TE+s*dx/ap1;

   end

end

   A=Tmatrix(sizet,:);

   Ans=A'

   toc