第九章热传导方程的差分解法(精品).ppt

《第九章热传导方程的差分解法(精品).ppt》由会员分享，可在线阅读，更多相关《第九章热传导方程的差分解法(精品).ppt（29页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第九章 热传导方程的差分解法n n9.1 9.1 热传导方程概述热传导方程概述 考虑三维空间的温度变化情况考虑三维空间的温度变化情况,设设 t t 时刻点时刻点(x,y,z)(x,y,z)处的温度为处的温度为u(x,y,z,t),u(x,y,z,t),则则 t t 时间内通过横截面积时间内通过横截面积 S S 传导的热量为传导的热量为(沿沿 n n 方方向向):):其中其中:K(x,y,z,t):K(x,y,z,t)是介质的热传导系数是介质的热传导系数,为温度梯度的法向量为温度梯度的法向量分量分量.取空间中的一个小区域取空间中的一个小区域 V,V,其边界面其边界面 S S 为一封闭曲面为一封闭

2、曲面.则则 t t1 1 到到 t t2 2 时刻通过包面时刻通过包面 S S 传入传入 V V 的热量为的热量为:由由高斯公式高斯公式:为哈密顿算子为哈密顿算子:设介质的比热容为设介质的比热容为 c,c,密度为密度为 ,则则 V V 内温度变化消耗的内温度变化消耗的热量热量:设设 V V 内部热源密度为内部热源密度为 F(x,y,z,t),F(x,y,z,t),则内部热源产生的热量则内部热源产生的热量为为:根据能量守恒原则根据能量守恒原则:Q:Q2 2=Q=Q1 1+Q+Q3 3即即:亦即亦即:若若 F(x,y,z,t)F(x,y,z,t)0,c,0,c,K,K,为常数为常数,则则:其中其中

3、:为拉普拉斯算子为拉普拉斯算子:所以热传导方程为所以热传导方程为:其中其中:K K c c.n n 9.2 一维热传导方程的差分解法一维热传导方程一维热传导方程:n n初值问题初值问题初值条件初值条件:n n初边值混合问题初边值混合问题初值条件初值条件:边值条件边值条件:(:(关于边界点关于边界点x=0 x=0和和x=l)x=l)第一类第一类.第二类第二类:第三类第三类:其中其中g g1 1(t),g(t),g2 2(t),(t),1 1(t),(t),2 2(t)(t)为给定函数为给定函数,要要求求 1 1(t)(t),2 2(t)(t),且不同时为零且不同时为零.设设空间的步长为空间的步长

4、为 h,h,时间的步长为时间的步长为 .把空间和时间离散化把空间和时间离散化:近似微分近似微分:故可定义故可定义:对空间一阶向前插商对空间一阶向前插商:对空间一阶向后插商对空间一阶向后插商:对空间二阶中心差商对空间二阶中心差商:对时间一阶向前插商对时间一阶向前插商:代入热传导方程代入热传导方程:迭代公式迭代公式:t t i-1 i i+1 x k+1 k第一类初边值条件第一类初边值条件:已知已知:第二类初边值条件第二类初边值条件:已知已知即即:计算过程计算过程:第三类初边值条件第三类初边值条件:已知已知:即即:例例1:1:差分方程差分方程:初边值条件初边值条件:function u=funct

5、ion u=rcd(lamda,tao,h,H,Trcd(lamda,tao,h,H,T)x=0:h:H;x=0:h:H;t=0:tao:T;t=0:tao:T;a=a=taotao*lamda/h2;*lamda/h2;N=length(x);N=length(x);M=length(t);M=length(t);u(:,1)=(4*x.*(1-x);u(:,1)=(4*x.*(1-x);u(1,2:M)=0;u(1,2:M)=0;u(N,2:M)=0;u(N,2:M)=0;for k=1:M-1for k=1:M-1 for i=2:N-1 for i=2:N-1 u(i,k+1)=a*u

6、(i+1,k)+(1-2*a)*u(i,k)+a*u(i-1,k);u(i,k+1)=a*u(i+1,k)+(1-2*a)*u(i,k)+a*u(i-1,k);end endendendh1=line(Color,1 0 h1=line(Color,1 0 0,Marker,.,MarkerSize,20,EraseMode,xor);0,Marker,.,MarkerSize,20,EraseMode,xor);for i=1:length(t)for i=1:length(t)set(h1,Xdata,0:0.1:1,Ydata,u(:,i);set(h1,Xdata,0:0.1:1,Yd

7、ata,u(:,i);pause(taopause(tao););endendX,Y=meshgrid(x,0:0.01:0.2);X,Y=meshgrid(x,0:0.01:0.2);Z=repmat(u(:,1),size(X,1),1);Z=repmat(u(:,1),size(X,1),1);h2=surface(X,Y,Z);h2=surface(X,Y,Z);shading shading interp,axisinterp,axis equal;equal;set(h2,EraseMode,xor);set(h2,EraseMode,xor);for i=1:length(t)f

8、or i=1:length(t)CD=repmat(u(:,i),size(X,1),1);CD=repmat(u(:,i),size(X,1),1);set(h2,Cdata,CD);set(h2,Cdata,CD);pause(taopause(tao););endendn n9.3 二维热传导方程的差分解法内部无热源均匀介质中二维热传导方程内部无热源均匀介质中二维热传导方程:初值条件初值条件:边值条件视具体情况而定边值条件视具体情况而定.设空间的步长为设空间的步长为 h,h,时间的步长为时间的步长为 .设设NhNh=l,=l,MhMh=s,=s,把时间和空间离散化把时间和空间离散化:即即

9、:微分近似微分近似:代入热传导方程得代入热传导方程得:k+1ki-1 i i+1j+1jj-1例例:初值条件初值条件:即即:恒温边界恒温边界:绝热边界绝热边界:即即:差分公式差分公式:function u=rcd2(lamda,tao,h,T,L,S)function u=rcd2(lamda,tao,h,T,L,S)%二维热传导方程二维热传导方程t=0:tao:T;t=0:tao:T;x=0:h:L;x=0:h:L;y=0:h:S;y=0:h:S;a=a=taotao*lamda/h2;*lamda/h2;if a0.25if a0.25 error(lamdaerror(lamda*tao

10、/h20.25);*tao/h20.25);endendD=length(t);D=length(t);N=length(x);N=length(x);M=length(y);M=length(y);M1=ceil(M/2)-3;M1=ceil(M/2)-3;M2=ceil(M/2)+3;M2=ceil(M/2)+3;u=zeros(N,M,D);u=zeros(N,M,D);u(:,:,1)=0;u(:,:,1)=0;u(:,1,:)=0;u(:,1,:)=0;u(:,M,:)=0;u(:,M,:)=0;u(1,M1:M2,2:D)=1;u(1,M1:M2,2:D)=1;for k=1:D-

11、1for k=1:D-1 for i=2:N-1 for i=2:N-1 for j=2:M-1 for j=2:M-1 u(i,j,k+1)=(1-4*a)*u(i,j,k)+a*(u(i+1,j,k)+u(i-u(i,j,k+1)=(1-4*a)*u(i,j,k)+a*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k);1,j,k)+u(i,j+1,k)+u(i,j-1,k);u(N,j,k+1)=u(N-1,j,k+1);u(N,j,k+1)=u(N-1,j,k+1);if(jM2)if(jM2)u(1,j,k+1)=u(2,j,k+1);u(1,j

12、,k+1)=u(2,j,k+1);end end end end end endendendX,Y=X,Y=meshgrid(x,ymeshgrid(x,y););Z=u(:,:,1);Z=u(:,:,1);h2=surface(X,Y,Z);h2=surface(X,Y,Z);shading shading interp,axisinterp,axis equal;equal;set(h2,EraseMode,xor);set(h2,EraseMode,xor);for i=1:length(t)for i=1:length(t)CD=u(:,:,i);CD=u(:,:,i);set(h2,Cdata,CD);set(h2,Cdata,CD);pause(taopause(tao););endendrcd2(1,0.001,0.1,0.2,2,1);rcd2(1,0.001,0.1,0.2,2,1);