热传导方程的差分解法c语言,九热传导方程的差分解法.PPT

九热传导方程的差分解法

第九章 热传导方程的差分解法 9.1 热传导方程概述 考虑三维空间的温度变化情况, 设 t 时刻点(x,y,z)处的温度为u(x,y,z,t), 则 ?t 时间内通过横截面积 ?S 传导的热量为(沿 n 方向): 其中: K(x,y,z,t) 是介质的热传导系数, 为温度梯度的法向量分量. 取空间中的一个小区域 V, 其边界面 S 为一封闭曲面. 则 t1 到 t2 时刻通过包面 S 传入 V 的热量为: 由高斯公式: 为哈密顿算子: 设介质的比热容为 c, 密度为 ?, 则 V 内温度变化消耗的热量: 设 V 内部热源密度为 F(x,y,z,t), 则内部热源产生的热量为: 根据能量守恒原则: Q2 = Q1 + Q3 即: 亦即: 若 F(x,y,z,t)?0, c, ?, K, 为常数,则: 其中: ? 为拉普拉斯算子: 所以热传导方程为: 其中: ? ? K?c?. 9.2 一维热传导方程的差分解法 一维热传导方程: 初值问题 初值条件: 初边值混合问题 初值条件: 边值条件:(关于边界点x=0和x=l) 第一类. 第二类: 第三类: 其中g1(t), g2(t), ?1(t), ?2(t) 为给定函数, 要求?1(t)??, ?2(t) ??, 且不同时为零. 设空间的步长为 h, 时间的步长为 ?. 把空间和时间离散化: 近似微分: 故可定义: 对空间一阶向前插商: 对空间一阶向后插商: 对空间二阶中心差商: 对时间一阶向前插商: 代入热传导方程: 迭代公式: t 第一类初边值条件: 已知: 第二类初边值条件: 已知 即: 计算过程: 第三类初边值条件: 已知: 即: 例1: 差分方程: 初边值条件: function u = rcd(lamda,tao,h,H,T) x = 0:h:H; t = 0:tao:T; a = tao*lamda/h^2; N = length(x); M = length(t); u(:,1) = (4*x.*(1-x))'; u(1,2:M) = 0; u(N,2:M) = 0; for k=1:M-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); end end h1=line('Color',[1 0 0],'Marker','.','MarkerSize',20,'EraseMode','xor'); for i=1:length(t) set(h1,'Xdata',[0:0.1:1],'Ydata',u(:,i)); pause(tao); end [X,Y]=meshgrid(x,[0:0.01:0.2]); Z = repmat(u(:,1)',size(X,1),1); h2 = surface(X,Y,Z); shading interp,axis equal; set(h2,'EraseMode','xor'); for i=1:length(t) CD = repmat(u(:,i)',size(X,1),1); set(h2,'Cdata',CD); pause(tao); end 9.3 二维热传导方程的差分解法 内部无热源均匀介质中二维热传导方程: 初值条件: 边值条件视具体情况而定. 设空间的步长为 h, 时间的步长为 ?. 设Nh=l, Mh=s, 把时间和空间离散化: 即: 微分近似: 代入热传导方程得: 例: 初值条件: 即: 恒温边界: 绝热边界: 即: 差分公式: function u = rcd2(lamda,tao,h,T,L,S) % 二维热传导方程 t = 0:tao:T; x = 0:h:L; y = 0:h:S; a = tao*lamda/h^2; if a>0.25 error('lamda*tao/h^2>0.25'); end D = length(t); N = length(x); M = length(y); M1 = ceil(M/2)-3; M2 = ceil(M/2)+3; u = zeros(N,M,D); u(:,:,1) = 0; u(:,1,:) = 0; u(:,M,:) = 0; u(1,M