偏微分一维热传导问题

1、偏微分大作业一维热传导方程问题运用隐式格式求解数值解目录问题描述31解析解分离变量法 32数值解隐式格式 53证明隐式格式的相容性与稳定性 54 数值解分析与 Matlab实现 65数值解与解析解的比较 96随时间变化的细杆上的温度分布情况 117稳定后细杆上的温度分布情况 13参考文献13附录13有限长杆的一维热传导问题问题描述一根单位长度的细杆放入100C的沸水中，当细杆的温度达到100C时取出。 假设细杆四周绝热；在时间t=0时，细杆两端浸入0°c的冰水中。一维热传导方2 2程：Ut a Uxx 0，现在令a 1，从而可知本题：Ut Uxx0。现在要求细杆温度分布：u(x,t)

2、。1解析解一一分离变量法热传导偏微分方程：Ut Uxx 0(1)u(O,t) U(1,t)0.U(x，)(x)其中，0， x 0或x 1(x)"L 100，x (0,1)首先令：U(x,t) X(x)T(t)将(2)式带入(1)式得：X(x)T(t) T(t)X(x) 0于是可得:l(t)T(t)X(x)X(x)可以得到两个微分方程：T(t)L |T(t) 0X(x)X(x) 0先求解空间项：当0时，X(x)Ae xBe x由于 u(0, t) U(1,t)0, t.可知：由于解的收敛性，B 0X(0)=X(1) A Ae 0 A 0 则此时是平庸解。当 0 时，X(x) A BxX

3、(0)=X(1) A A B 0 A 0,B0则此时是平庸解。当 0 时，X (x) A cos kx Bsinkx,其中 k 。X(0) A 0 A 0X(1) Bsink 0 k n ,n 1,2,3所以，X(x)Bnsin(nx)，n 1,2,3因为2 2 n所以，T(t)n2 2tCne，n1,2,3则，u(x,t)nn2 2tDne1sin(nx)初始条件： u(x,)(x)u(x，)Dn sin(n x) (x)n 11Dn 2 0 (x)sin(n x)dxi2100sin(n x)dx 1200 () cosn (1) cosnn当0时，m Dn = 200 (1 cosn )

4、 n / 最终，u(x,t)200(1 1)n)e n t sin(n x)， n 1,2,3n 1 n2数值解隐式格式目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式 格式1。利用u（X,t），关于t进行向前差商：Uk1 Ukk 1j厂；关于X进行二阶中心差:Uk； 2Uf U； （X ；代入偏微分方程可以得到隐式差分格式：U k1 U； U k 1(1)3证明隐式格式的相容性与稳定性(1)相容性根据Taylor展开: 1=u：+ +t 1 uk+*C U：U-1I I k 1 k UI Uj 1 Uj+12t2°12U21Ju22t22 t2（护）x22 x&qu

5、ot;1 2U 2 EX（k3）（k3）3X )2U T + X(水2)(jx)代入隐式格式得：u 1 Jut 2 t2将（2）与原微分方程相减，得到截断误差=1=2所以此隐式格式与原微分方程相容。（2）稳定性首先令:令网格比为r则可以将（1）式改写得到:tx2，rU将（4）代入（3）(1 2r)Ujk1 rU：1 U：Uk1eI(j-DUk；U：U：vU：+1式，根据欧拉公式化简得:Uk 1+2r(3)Uk1eI(4)（j+i）2rcos ) U k(5)G叮 u：故得放大因子是:11+2r（1 cos ）所以根据Fourier方法，隐式格式恒稳定。4数值解一一分析与 Matlab实现（1）

6、边值与初值离散化将边值与初值离散化，与式（3）联立得差分线性方程组:r rukld2r)Ujk1 rU：11 U：，j (0,1,2k (0,1,2,|,M|,N-1)1)U0 =(Xj),j(0,1,2, |,M)U :U00,k(0,1,2, (j|,N)< u：0，k(0,1,2, |，N)再将方程组改写成AUB的形式：1+2r rr 1 2r rr 1 2rUik1Uik+rU0Tu；1u：*114r 1 2rruk 1M 2ukM2r1 2r(M 1)(m 1)uk 1M 1ukM 1k 1rUM本题的边界条件均为零。所以可以将上式改写。/1+2rruk1/r1 2rrurr1

7、 2r41 1u：1*I1r 1 2rk 1rU M 2r 12ru k 1(M 1)(M 1)M 1u2ku：uM 2uM 1(2) Matlab的实现? 杆长1米，时间2秒。设计空间步长h=和时间步长t=，网格比是rh2从而得到划分的空间网格点数是 M1+1，时间网格点数是 M2+1。先设初 始的温度矩阵U(M2+1,M1+1)。再将边界条件和初始条件编写到表示温度分 布的矩阵中。具体代码可见最后附录。? 编写矩阵A核心代码：对角线：A(i,i) = 1+2r对角线的右方和下方：A(i,i+1) = -r;/A(i+1,i) = -r;/? 下面就要运用A*u(k 1,j) u (k, j

8、)进行迭代。当 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。u：1u：? 数值解画图，如图1(a)和图1(b)所示。维熱传导方程-数煩解-温度分布團.50.图1(a)数值解的温度分布图现在将着色平稳过渡。维如专导方园魏值解温度分布團_|snsn4020n-60.2CT-图1(b)着色平稳过渡的数值解的温度分布图5数值解与解析解的比较2 2首先，我们需要将解析解离散化，解析解中有一项e n t，当n越来越大 时，会快速趋于0，故我们可以取

9、n=8000。现在来证明可行性，在 matlab 里的工作空间运算。将解析解的温度分布画出来，数值解画图2，如图2所示。维热传导方程解折解温度分布图图2解析解的温度分布图_|I130504020O 2将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a)上可以看出空间上的误差，在边界处误差比较大维热传导方程.误差“温度分布圉位亂图3( a)数值解与解析解空间误差我们从图3(a)上可以看出时间的误差，在时间的最开始，处误差最大，然后 又有一个小的波动，最后就误差渐渐变小，最后趋于 0-维熱伟导方程“误差温度分布图,JLibJLjJJ00.20.40.60.S11.21.41

10、.6132时间t图3 (b) 数值解与解析解时间误差6随时间变化的细杆上的温度分布情况从数值解的温度分布三维图，如图4(a)和图4(b)可以看出随着时间的增加， 细杆温度下降最后趋于0°C。从物理角度来说：细杆的温度会不断地向两端扩散， 热量会慢慢散失，最终 随着时间的增加，细杆的温度会趋于 0C。錐熬传导方程数值解-温度分布團00.20.4 D6 0.811.21.41.6132时间t图4(a)细杆温度随时间的变化图现取细杆中心处一点，观看它随时间的温度变化情况。图4(b)细杆中央(x=)温度随时间的变化图7稳定后细杆上的温度分布情况从图像上可以看出，最后稳定的情况下，细杆的温度是

11、0C参考文献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) = 0%u(L,0) = 0%其中,aA2 = 1, L = 1%clc;clear all ;%区域及划分网格L = 1;%单位长度的细杆？T = 2;%时间h = ;%空间的划分 %t = ;%时间的

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

13、 = 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); end ;%k+1时刻的不同位置的温度分布end ;%乍图x = 0:h:1;y = 0:t:2;xx,yy=meshgrid(x,y);figure(1);

14、surf(xx,yy,U);shad ingflattitle(' 一维热传导方程-数值解-温度分布图);xlabel('位置x');ylabel('时间 t');zlabel('温度T');figure(2)s = 0;for i= 1:8000s = s+(200*(1-(-1)Ai)/(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(' 一维热传导方程-误差-温度分布图)；xlabel( '位置 x');ylabel('时间