热传导方程的差分解法省公共课一等奖全国赛课获奖课件.pptx

第九章 热传导方程差分解法第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 t1 1 到到 t t2 2 时刻经过包面时刻经过包面 S S 传入传入 V V 热量为热量为:第2页由高斯公式由高斯公式:n n为哈密顿算子为哈密顿算子:设介质比热容为设介质比热容为 c,c,密度为密度为 ,则则 V V 内温度改变消耗热量内温度改变消耗热量:设设 V V 内部热源密度为内部热源密度为 F(x,y,z,t),F(x,y,z,t),则内部热源产生热量为则内部热源产生热量为:第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页其中其中:为拉普拉斯算子为拉普拉斯算子:所以热传导方程为所以热传导方程为:其中其中:K K c c.第5页n n 9.2 一维热传导方程差分解法一维热传导方程一维热传导方程:n n初值问题初值问题初值条件初值条件:n n初边值混合问题初边值混合问题初值条件初值条件:边值条件边值条件:(:(关于边界点关于边界点x=0 x=0和和x=l)x=l)第一类第一类.第6页第二类第二类:第三类第三类:其中其中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),且不一样时为零且不一样时为零.第7页设空间步长为设空间步长为 h,h,时间步长为时间步长为 .把空间和时间离散化把空间和时间离散化:近似微分近似微分:故可定义故可定义:对空间一阶向前插商对空间一阶向前插商:第8页对空间一阶向后插商对空间一阶向后插商:对空间二阶中心差商对空间二阶中心差商:对时间一阶向前插商对时间一阶向前插商:第9页代入热传导方程代入热传导方程:迭代公式迭代公式:第10页 t t i-1 i i+1 x k+1 k第11页第一类初边值条件第一类初边值条件:已知已知:第12页第二类初边值条件第二类初边值条件:已知已知即即:第13页计算过程计算过程:第14页第三类初边值条件第三类初边值条件:已知已知:即即:第15页例例1:1:差分方程差分方程:初边值条件初边值条件:第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: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页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,Xdata,0:0.1:1,Ydata,u(:,i);set(h1,Xdata,0:0.1:1,Ydata,u(:,i);pause(tao);pause(tao);endend第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(h2,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页n n9.3 二维热传导方程差分解法内部无热源均匀介质中二维热传导方程内部无热源均匀介质中二维热传导方程:初值条件初值条件:边值条件视详细情况而定边值条件视详细情况而定.设空间步长为设空间步长为 h,h,时间步长为时间步长为 .设设Nh=l,Mh=s,Nh=l,Mh=s,把把时间和空间离散化时间和空间离散化:第20页即即:微分近似微分近似:第21页代入热传导方程得代入热传导方程得:k+1ki-1 i i+1j+1jj-1第22页例例:初值条件初值条件:即即:恒温边界恒温边界:绝热边界绝热边界:即即:第23页差分公式差分公式:第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: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页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(:,:,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页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,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页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页rcd2(1,0.001,0.1,0.2,2,1);rcd2(1,0.001,0.1,0.2,2,1);第29页