1、第九章 热传导方程差分解法第1页第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 为一封闭曲面为一封闭曲面.则则 t
2、 t1 1 到到 t t2 2 时刻通过包面时刻通过包面 S S 传入传入 V V 热量为热量为:第2页第2页由高斯公式由高斯公式:n n为哈密顿算子为哈密顿算子:设介质比热容为设介质比热容为 c,c,密度为密度为 ,则则 V V 内温度改变消耗热量内温度改变消耗热量:设设 V V 内部热源密度为内部热源密度为 F(x,y,z,t),F(x,y,z,t),则内部热源产生热量为则内部热源产生热量为:第3页第3页依据能量守恒原则依据能量守恒原则: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,为常数为常数,则则:第4页第
3、4页其中其中:为拉普拉斯算子为拉普拉斯算子:因此热传导方程为因此热传导方程为:其中其中:K K c c.第5页第5页n n 9.2 一维热传导方程差分解法一维热传导方程一维热传导方程:n n初值问题初值问题初值条件初值条件:n n初边值混合问题初边值混合问题初值条件初值条件:边值条件边值条件:(:(关于边界点关于边界点x=0 x=0和和x=l)x=l)第一类第一类.第6页第6页n n第二类:n n第三类:n n其中g1(t),g2(t),1(t),2(t)为给定函数,要求1(t),2(t),且不同时为零.第7页第7页设空间步长为设空间步长为 h,h,时间步长为时间步长为 .把空间和时间离散化把
4、空间和时间离散化:近似微分近似微分:故可定义故可定义:对空间一阶向前插商对空间一阶向前插商:第8页第8页对空间一阶向后插商对空间一阶向后插商:对空间二阶中心差商对空间二阶中心差商:对时间一阶向前插商对时间一阶向前插商:第9页第9页代入热传导方程代入热传导方程:迭代公式迭代公式:第10页第10页 t t i-1 i i+1 x k+1 k第11页第11页第一类初边值条件第一类初边值条件:已知已知:第12页第12页第二类初边值条件第二类初边值条件:已知已知即即:第13页第13页计算过程计算过程:第14页第14页第三类初边值条件第三类初边值条件:已知已知:即即:第15页第15页例例1:1:差分方程差
5、分方程:初边值条件初边值条件:第16页第16页function u=rcd(lamda,tao,h,H,T)function u=rcd(lamda,tao,h,H,T)x=0:h:H;x=0:h:H;t=0:tao:T;t=0:tao:T;a=tao*lamda/h2;a=tao*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:
6、M-1 for i=2:N-1 for i=2:N-1 u(i,k+1)=a*u(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 endendend第17页第17页h1=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,Xd
7、ata,0:0.1:1,Ydata,u(:,i);set(h1,Xdata,0:0.1:1,Ydata,u(:,i);pause(tao);pause(tao);endend第18页第18页X,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 interp,axis equal;shading interp,axis equal;set(h
8、2,EraseMode,xor);set(h2,EraseMode,xor);for i=1:length(t)for 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(tao);pause(tao);endend第19页第19页n n9.3 二维热传导方程差分解法内部无热源均匀介质中二维热传导方程内部无热源均匀介质中二维热传导方程:初值条件初值条件:边值条件视详细情况而定边值条件视详细情况而定.设空间步长为设空间步长为
9、 h,h,时间步长为时间步长为 .设设Nh=l,Mh=s,Nh=l,Mh=s,把把时间和空间离散化时间和空间离散化:第20页第20页即即:微分近似微分近似:第21页第21页代入热传导方程得代入热传导方程得:k+1ki-1 i i+1j+1jj-1第22页第22页例例:初值条件初值条件:即即:恒温边界恒温边界:绝热边界绝热边界:即即:第23页第23页差分公式差分公式:第24页第24页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
10、:L;x=0:h:L;y=0:h:S;y=0:h:S;a=tao*lamda/h2;a=tao*lamda/h2;if a0.25if a0.25 error(lamda*tao/h20.25);error(lamda*tao/h20.25);endend第25页第25页D=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(
11、:,:,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;第26页第26页for k=1:D-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,
12、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,k+1)=u(2,j,k+1);end end end end end endendend第27页第27页X,Y=meshgrid(x,y);X,Y=meshgrid(x,y);Z=u(:,:,1);Z=u(:,:,1);h2=surface(X,Y,Z);h2=surface(X,Y,Z);shading interp,axis equal;shading interp,axis 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(tao);pause(tao);endend第28页第28页rcd2(1,0.001,0.1,0.2,2,1);rcd2(1,0.001,0.1,0.2,2,1);第29页第29页
©2010-2025 宁波自信网络信息技术有限公司 版权所有
客服电话:4008-655-100 投诉/维权电话:4009-655-100