matlab pde热传导,有限差分和Matlabpde求解一维稳态传热问题

《有限差分和Matlabpde求解一维稳态传热问题》由会员分享，可在线阅读，更多相关《有限差分和Matlabpde求解一维稳态传热问题(8页珍藏版)》请在人人文库网上搜索。

1、有限差分和pde函数求解一维定态热传导方程分别用有限差分方法和pde函数求解一维定态热传导方程，初始条件和边界条件，热扩散系数=0.00001，(1)求解过程： 1. 用Tylaor展开法推导出FTCS格式的差分方程首先对T进行泰勒展开得到如下两式子：上述两个方程变换得：(2)(3)将上述式子(2)(3)代入(1)得：(4)2. 方程的相容性和稳定性讨论：上述方程截项为：由于所以方程有相容性其经过傅里叶变换后，只有当满足下列条件时，方程具有较好的稳定性：其中N为节点个数，L为边界长度由于：所以当时方程具有稳定性3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。由式(4)得：要想知。

2、道下一个时刻下的温度，则必须知道在时刻下该点以及相邻两点的温度，如下图所示：所以需要知道的是在各个时刻x轴两侧的温度，以及初始时刻x轴上各点温度。下述各题求解过程中，均采用 ，作为边界条件。边界条件是由外界给定，本身对过程求解，精度无影响4. 编写M文件求解上述方程M文件的编写见附录所示。在matlab中调用heat_conduct(0.00001,15000,500,0.1)，得到如下结果：rms error isrms =0.7959说明：1.上述结果采用的分别为0.1和500，时间最终给定为15000当为(0.01,5)时的结果附录1，传热方程有限差分的 M文件编写function he。

3、at_conduct(alpha,Timemax,dt,dx)%方程需要输入的参数有，最大时刻数，alpha，dt，dx,实际Tmax=Timemax/dt+1;%时刻数JMAX=1/dx+1;%x轴上的差分点个数J=JMAX-1;MAXEX=J;T=zeros(Tmax,JMAX);%构造温度矩阵TE=zeros(Tmax,JMAX);%实际温度矩阵%设定边界条件T(:,1)=100;T(:,JMAX)=100;s=alpha*dt/dx/dx;%计算各个点的温度for i=1:(Tmax-1)for j=2:JT(i+1,j)=s*T(i,j+1)+(1-2*s)*T(i,j)+s*T(i。

4、,j-1);endend%计算实际值%画图z=0:dx:1;t=0:(Tmax-1);a=3000/dt+1;b=9000/dt+1;subplot(1,2,1)plot(z,T(a,:),g*,z,T(a,:),g,z,T(b,:),o,z,T(b,:),r,z,T(Tmax,:),b)xlabel(x方向)ylabel(温度T)title(3000,9000,15000时刻的温度分布)%画整体分布图for i=1:JMAXX(:,i)=t;end;for j=1:TmaxY(j,:)=z;endsubplot(1,2,2)mesh(X,Y,T);view(-1 1 1);xlabel(时刻。

5、t);ylabel(方向x);zlabel(温度T);title(整体分布)二、调用MATLAB函数完成上述计算1. 编写M文件求解上述方程，并用适当的文字对程序做出说明。Matlab中自带的pdepe函数可直接求解一般偏微分方程。其调用格式为：T=pdepe(m,pdefun,pdeic,pdebc,x,t);其中T为输出结果，在本题中，其输出格式为一个二维数组，分别表示一维空间和时间上的温度；x，t，分别为所需输入的时间和空间数组；为达到该pdepe函数求解要求，将上述方程转化为如下函数：即pdefun所要求的函数的描写格式：本题中，在描述pdefun时，；pdeic为初值条件，其描述的标。

6、准格式为：pdebc为边界条件，本题要求边界温度为100，其描述的标注格式为：所以左右边界的与分别表示为：；上述程序编写见附录2heat_conduct1(0.1,500,15000)得到如下结果所示：附录2:Matlab PDE求解上述一维传热方程%带求解函数function c,f,s=pdefun(x,t,T,dT)c=1;f=0.00001*dT;s=0;%边值条件，a为左边界，b为右边界，都为100function pa,qa,pb,qb=pdebc(xa,Ta,xb,Tb,t)pa=100-Ta;qa=0;pb=100-Tb;qb=0;%初值设定为0function T0=pdei。

7、c(x)T0=0;调用求解函数function heat_conduct1(dx,dt,Timemax)clcx=0:dx:1;t=0:dt:Timemax;m=0;T=pdepe(m,pdefun,pdeic,pdebc,x,t);figure(numbertitle,off,name,PDE Demoby Matlabsky)subplot(1,2,2)surf(x,t,T(:,:)title(solution)xlabel(X)ylabel(t)zlabel(T)subplot(1,2,1)plot(x,T(7,:),b,x,T(19,:),g*,x,T(19,:),g,x,T(31,:),ro,x,T(31,:),r)title(3000,9000,15000时刻温度)xlabel(X)ylabel(T)end。