偏微分一维热传导问题

《偏微分一维热传导问题》由会员分享，可在线阅读，更多相关《偏微分一维热传导问题（9页珍藏版）》请在装配图网上搜索。

1、文档收集于互联网，已重新整理排版.word版本可编辑欢迎下载支持.偏微分大作业一维热传导方程问题运用隐式格式求解数值解目录问题描述31解析解分离变量法32数值解一一隐式格式53证明隐式格式的相容性与稳定性54数值解分析与Matlab实现65数值解与解析解的比较96随时间变化的细杆上的温度分布情况117稳定后细杆上的温度分布情况12参考文献13附录14有限长杆的一维热传导问题问题描述一根单位长度的细杆放入100C的沸水中，当细杆的温度达到100C时取出。 假设细杆四周绝热；在时间t=0时，细杆两端浸入oc的冰水中。一维热传导方 程：ut_a2uxx = 0，现在令a2 =1，从而可知本题：ut_

2、uxx二0。现在 要求细杆温度分布： u(x, t) 。1 解析解分离变量法热传导偏微分方程：(1)(2)u u 0t xxu (x, t) = X (x )T (t)将(2)式带入(1)式得：于是可得：可以得到两个微分方程：先求解空间项：当九 编写矩阵A核心代码：对角线： A(i,i) = i+2r对角线的右方和下方： A(i,i+i) = -r;A(i+i,i) = -r;下面就要运用A * U(k +i,j) = U(k，j)进行迭代。当 k=i 时， A*U(2,j)=U(i,j)当 k=2 时， A*U(3,j)=U(2,j)当 k=3 时， A*U(4,j)=U(3,j)以此迭代下

3、去直到k=M2。就可以得到整个温度随时间和空间的分布矩阵 U。 数值解画图，如图1(a)和图1(b)所示。图 1(a) 数值解的温度分布图现在将着色平稳过渡。图 1(b) 着色平稳过渡的数值解的温度分布图5 数值解与解析解的比较 首先，我们需要将解析解离散化，解析解中有一项e -血2t，当n越来越大 时，会快速趋于0,故我们可以取n=8000。现在来证明可行性，在matlab 里的工作空间运算。将解析解的温度分布画出来，数值解画图切，如图2所示。图 2 解析解的温度分布图将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a) 上可以看出空间上的误差，在边界处误差比较大。图

4、 3(a) 数值解与解析解空间误差我们从图3(a)上可以看出时间的误差，在时间的最开始，处误差最大，然后 又有一个小的波动，最后就误差渐渐变小，最后趋于0。图 3(b) 数值解与解析解时间误差6 随时间变化的细杆上的温度分布情况从数值解的温度分布三维图，如图4(a)和图4(b)可以看出随着时间的增加， 细杆温度下降最后趋于0C。从物理角度来说：细杆的温度会不断地向两端扩散，热量会慢慢散失，最终 随着时间的增加，细杆的温度会趋于0C。图 4(a) 细杆温度随时间的变化图现取细杆中心处一点，观看它随时间的温度变化情况。图4(b)细杆中央(x=0.5)温度随时间的变化图7 稳定后细杆上的温度分布情况

5、从图像上可以看出，最后稳定的情况下，细杆的温度是 0C。参考文献1 冯立伟.热传导方程几种差分格式的MATLAB数值解法的比较J.沈阳化 工大学，辽宁沈阳2011(6)2 一维热传导方程数值解法及Matlab实现EB/OL. 2014-11-20 -945.html附录代码 ：%此程序用于解决一维热传导方程：ut-a2uxx = 0%边界条件：u(0,t) = u(L,t) = 0%初始条件： u(x,0) = 100, x!=0 和 L%u(0,0) = 0%u(L,0) = 0%其中,aA2 = 1, L = 1% % 单位长度的细杆 % 时间% 空间的划分 % % 时间的划分% 网格比c

6、lc; clear all; %区域及划分网格L = 1; T = 2;h = 0.1; t = 0.01;r = t/(h*h);构%造的矩阵： U( 时间，空间) %编程包含边值，如 U(k,1)=u(0,t) % 时间划分了 M2 份，有 M2+1 个节点 %两个边界处温度恒为零%设计步长 M1 = L/h; M2 = T/t;%构造边界条件U = zeros(M2+1,M1+1); for k = 1：M2+1U(k,1) = 0; U(k,M1+1) = 0;end;%构造初始条件for j = 2：M1% 位置划分了 M1 份，有 M1+1 个节点U(1,j) = 100;end;

7、U(1,1) = 0;U(1,M1+1) = 0;%差分格式的矩阵形式 A*U(k+1,j)=U(k,j)%构造矩阵AA = zeros(M1-1);for i = 1:M1-1A(i,i) = 1+2*r;end;for i = 1:M1-2A(i,i+1) = -r;A(i+1,i) = -r;end;%构造AU=B中的B本题边值的特殊，矩阵B大大简化了B = zeros(M1-1,1);for k = 1:M2j = 2:M1;B(j-1,1) = U(k,j);x = AB;for j = 2:M1U(k+1,j) = x(j-1); %k+1 时刻的不同位置的温度分布end;end;

8、%作图x = 0:h:1;y = 0:t:2;xx,yy=meshgrid(x,y);figure(1); surf(xx,yy,U);shading flattitle (维热传导方程-数值解一温度分布图);xlabel( 位置 x);ylabel(时I、可 t);zlabel( 温度 T);figure(2)s = 0;for i= 1:8000s = s+(200*(1-(-1)Ai)/(i*pi)*sin(i*pi*xx).*exp(-iA2*piA2*yy); end;surf(xx,yy,s);title (维热传导方程-解析解-温度分布图); xlabel(位置 x);ylabe

9、l( 时、 t); zlabel( 温度 T);figure(3)x = 0:h:1;y = 0:t:2;xx,yy = meshgrid(x,y);dd = U-s; surf(xx,yy,dd);title( 一维热传导方程-误差-温度分布图); xlabel( 位置 x);ylabel(时I、可 t);zlabel( 误差(数值解减解析解) ) ;figure(4)z = zeros(M2+1,1);if mod(M1+1),2)=0i = 1:M2+1; z(i,1) = U(i,M1/2);elsei = 1:M2+1; z(i,1) = U(i,(M1+1)/2);end; plot(z); title(温度随时间增加的趋势图); xlabel( 时、 t); ylabel( 温度 T);