《偏微分一维热传导问题.doc》由会员分享，可在线阅读，更多相关《偏微分一维热传导问题.doc（20页珍藏版）》请在三一文库上搜索。
1、偏微分大作业一维热传导方程问题运用隐式格式求解数值解目录问题描述 31 解析解分离变量法 32 数值解隐式格式 53 证明隐式格式的相容性与稳定性 54 数值解分析与 Matlab 实现 75 数值解与解析解的比较 106 随时间变化的细杆上的温度分布情况 127 稳定后细杆上的温度分布情况 14参考文献 14附录 16有限长杆的一维热传导问题问题描述100 C时取一维热传uxx 0。一根单位长度的细杆放入100 C的沸水中，当细杆的温度达到 出。假设细杆四周绝热；在时间t=0时，细杆两端浸入0 C的冰水中2 2导方程：ut a Uxx 0，现在令a 1，从而可知本题：ut 现在要求细杆温度分
2、布：u(x,t)。1解析解一一分离变量法热传导偏微分方程：(1)Ut Uxx 0u(0,t)u(1,t)0I u(x，)(x)其中，0， x 0或x 1 (x)<100，x (0,1)首先令：u(x,t) X(x)T(t)将式带入式得：X(x)T&t) T(t)X(x) 0于是可得：T&t)T(t) X(x)可以得到两个微分方程：T(t) T(t) 0瓠)X(x) 0先求解空间项:当 0 时，X(x) Ae X Be X由于 u(0, t) u(1,t)0, t.可知：由于解的收敛性，B 0X(0)=X(1) A Ae 0 A 0则此时是平庸解。当 0 时，X(x) A
3、BxX(0)=X(1) A A B 0 A 0,B0则此时是平庸解。当 0 时，X (x) A cos kx Bsinkx,其中 k 。X(0) A 0 A 0所以，X(x)Bnsin(n x)， n 1,2,3因为n2 2所以，T(t)n2 2tCne，n 1,2,3则，X(1) Bsi nk 0k n ,n 1,2,3n2 2tu(x,t)Dnesin(n x)n 1初始条件：u(x，)(x)u(x?0)Dn sin(n x) (x)n 11Dn 2 0 (x)sin(n x)dx100sin(n x)dx最终,200 (丄)cosn (1n)cosn当0时 Jim Dn 二200 cos
4、n )n2002 2u(x,t) (1 ( 1)n)e n t sin(n x), n 1,2,3n 1 n2数值解隐式格式目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式 格式1。利用u(x,t)，关于t进行向前差商:uk1u：Vt；关于x进行二阶中心k 1k 1k 1差:U j 1 2U j U j 1 (Vxp代入偏微分方程可以得到隐式差分格式:u： 1 u：Vtuk2Ujk 1uk1(Vx)23证明隐式格式的相容性与稳定性(1)相容性根据Taylor展开:1 k U1=U:+ Vt J t1 UkVt J tu：; uk+-VtJ1 J tU：U：12t212U2t21
5、2U2t2Vt2Vt21 2Uo(Vt2)上Vxx卫Vxx1 2U22Vx22 x21 皀 Vx22 x2(Vt3)(Vt3)代入隐式格式得:(Vt2)=2Ux(Vt2)(Vx)(Vx3)(Vx3)将(2)与原微分方程相减，得到截断误差Vt- (Vx)所以此隐式格式与原微分方程相容。(2)稳定性Vt令网格比为r2，则可以将(1)式改写得到：Vx2k 1k 1k 1krU j 1(1 2r)U j rU j 1 U j(3)首先令:ruk； Uk1eI(M U k 1 U k 1eI j U k U keI jLuk+； Uk1eI(j+1)将(4)代入(3 )式，根据欧拉公式化简得:Uk (1
6、+2r2rcos )Uk(5)故得放大因子是:Uk1Uk11+2r(1 cos )所以根据Fourier方法，隐式格式恒稳定4数值解一一分析与 Matlab实现（1）边值与初值离散化k 1k 1k 1krU j 1(1 2r)U j rU j 1 U j将边值与初值离散化，与式（3）联立得差分线性方程组:j (0,1,2,L ,M -1)k (0,1,2,L ,N-1)U0=：(Xj),j(0,1,2,L,M)Uo0,k(0,1,2,L,N)U：0，k(0,1,2 丄,N)1+2r1 2r2rU：1u；1Uk1Uk+rU0k1U：U：1 2r1 2r(MU k 1M 2U k 11)(M 1)
7、 MU kM 2kk 1U M 1 rU M本题的边界条件均为零。所以可以将上式改写。1+2rrUk1U：r 12rru；1U；r 1 2rU：1U：OOOMMr 1 2rruM12UM 2r 12r(M 1)(M 1)uM11uM 1(2) Matlab的实现? 杆长1米，1时间2秒。设计空间步长h=0.1和时间步长t=0.01，网格比是tr 2。h从而得到划分的空间网格点数是 M1+1 ,时间网格点数是M2+1。先设初 始的温度矩阵U(M2+1,M1+1)。再将边界条件和初始条件编写到表示温度 分布的矩阵中。具体代码可见最后附录。编写矩阵A核心代码：对角线：A(i,i) = 1+2r对角线
8、的右方和下方：A(i,i+1) = -r;A(i+1,i) = -r;F面就要运用A*U(k 1,j) U (k, j)进行迭代。 当 k=1 时，A*U(2,j)=U(1,j)当 k=2 时，A*U(3,j)=U(2,j)当 k=3 时，A*U(4,j)=U(3,j)以此迭代下去直到k=M2。就可以得到整个温度随时间和空间的分布矩阵U。? 数值解画图，如图1(a)和图1(b)所示。一维熱传导方程数值解-温度分布團o o o o o O 21000 6 A 2图1(a)数值解的温度分布图现在将着色平稳过渡。图1(b)着色平稳过渡的数值解的温度分布图维热去导万瑋魏殖解”温度兮布用OO8O.S 位
9、.6 O2 D_|504020O 2-1n越来越在 matlab5数值解与解析解的比较2 2?首先，我们需要将解析解离散化，解析解中有一项e n t，当大时，会快速趋于0，故我们可以取n=8000。现在来证明可行性， 里的工作空间运算。将解析解的温度分布画出来，数值解画图，如图2所示图2 解析解的温度分布图维热去导万瑋.餅析觥”温度分布用_|30504020O 21将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a)上可以看出空间上的误差，在边界处误差比较大。一维热传导方程误差-温度分布團O 5 O £ 2图3( a) 数值解与解析解空间误差我们从图3(a)上
10、可以看出时间的误差，在时间的最开始，处误差最大，然后又有一个小的波动，最后就误差渐渐变小，最后趋于0。一维熱伟导方稈“误差温度分布團10； I-Be4200 0.2D.81 1J 1.41.6132时间t图3 (b)数值解与解析解时间误差6随时间变化的细杆上的温度分布情况从数值解的温度分布三维图，如图 4(a)和图4(b)可以看出随着时间的增加， 细杆温度下降最后趋于0 °C。从物理角度来说：细杆的温度会不断地向两端扩散， 热量会慢慢散失，最终 随着时间的增加，细杆的温度会趋于 0 C錐熬传导方程数值解“温度分布團4000.20.40.60.S1 1J 1.41.6132时间t图4(
11、a)细杆温度随时间的变化图现取细杆中心处一点，观看它随时间的温度变化情况图4(b)细杆中央(x=0.5 )温度随时间的变化图7稳定后细杆上的温度分布情况从图像上可以看出，最后稳定的情况下，细杆的温度是0 c参考文献1 冯立伟热传导方程几种差分格式的 MATLAB数值解法的比较J 沈阳化 工大学，辽宁沈阳.2011 (6).2 一维热传导方程数值解法及 Matlab实现EB/OL . 2014-11-20 附录代码：%此程序用于解决一维热传导方程:ut-aA2uxx = 0%边界条件： u(0,t) = u(L,t) = 0%初始条件： u(x,0) = 100, x!=0 和 L%u(0,0)
12、 = 0%u(L,0) = 0% 其中 ,aA2 = 1, L = 1%clc;clear all ;% 区域及划分网格L = 1; %单位长度的细杆 ?T = 2;%时间h = 0.1;%空间的划分%t = 0.01;%时间的划分%r = t/(h*h);%网格比%设计步长M1 = L/h;M2 = T/t;%构造边界条件U = zeros(M2+1,M1+1); for k = 1:M2+1U(k,1) = 0;%构造的矩阵： U( 时间，空间)% 编程包含边值，如 U(k,1)=u(0,t)% 时间划分了 M2 份，有 M2+1 个节点% 两个边界处温度恒为零U(k,M1+1) = 0;
13、 end ;%构造初始条件for j = 2:M1U(1,j) = 100; end ;U(1,1) = 0;U(1,M1+1) = 0;% 位置划分了 M1 份，有 M1+1 个节点% 差分格式的矩阵形式 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;
14、B(j-1,1) = U(k,j);x = AB;for j = 2:M1U(k+1,j) = x(j-1); %k+1 时刻的不同位置的温度分布 end ;end ;%作图x = 0:h:1;y = 0:t:2;xx,yy=meshgrid(x,y);figure(1);surf(xx,yy,U); shading flat title( '一维热传导方程 -数值解 - 温度分布图 '); xlabel( '位置 x');ylabel( '时间 t');zlabel( '温度 T');figure(2)s = 0;for i=
15、1:8000s = s+(200*(1-(-1)")/(i*pi)*s in (i*pi*xx).*exp(-22*piA2*yy); end ;surf(xx,yy,s);title( '一维热传导方程 -解析解 -温度分布图 '); xlabel( '位置 x');ylabel( '时间 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( '一维热传导方程
16、-误差 -温度分布图 '); xlabel( '位置 x');ylabel( '时间 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');