一维非稳态导热方程-热流绝热边界(附Matlab程序)

functiononed%根据网上一维算例改了下边界条件，边界条件为绝热和恒定热流%byhxgclearall;clc;%%%%%%%%%%%%%%需要输入的物性参数Lambda=10;%导热系数cp=440;%热容rou=7800;%密度qw=500000;%壁面热流%%%%%%%%%%%%%a=Lambda/rou/cp;%定义中间系数c=qw/Lambda;%定义中间系数xspan=[00.012];%轴向坐标起止位置tspan=[010];%仿真时间起止ngrid=[100020];%空间网格数和时间网格数%%%%%%%%%%%%%%调用子函数[T,x,t]=rechuandao(a,c,xspan,tspan,ngrid);%画图[x,t]=meshgrid(x,t);figure(1)mesh(x,t,T);xlabel('x')ylabel('t')zlabel('T')function[U,x,t]=rechuandao(a,c,xspan,tspan,ngrid)%热传导方程：%Ut(x,t)=c^2*Uxx(x,t)axbtsttf%一维非稳态导热方程%参数说明%c：方程中的系数%f：初值条件%g1,g2：边值条件%xspan=[a,b]：x的取值范围%tspan=[ts,tf]：t的取值范围%ngrid=[n,m]：网格数量，m为x网格点数量，n为t的网格点数量%U：方程的数值解%x,t：x和t的网格点n=ngrid(1);m=ngrid(2);h=range(xspan)/(m-1);x=linspace(xspan(1),xspan(2),m);k=range(tspan)/(n-1);t=linspace(tspan(1),tspan(2),n);r=a*k/h^2;%%ifr0.5error('为了保证算法的收敛，请增大步长h或减小步长k!')end%两个中间系数s=1-2*r;s2=c*h;U=zeros(ngrid);%初步赋初值%初值条件，给了一个均温场U(1,:)=20;%差分计算forj=2:nfori=2:m-1U(j,i)=s*U(j-1,i)+r*(U(j-1,i-1)+U(j-1,i+1));end%%%%%%%下面是边界条件U(j,1)=s2+U(j-1,2);%定热流边界条件U(j,m)=U(j,m-1);%绝热边界条件end