还剩14页未读,继续阅读
本资源只提供10页预览,全部文档请下载后查看!喜欢就下载吧,查找使用更方便
文本内容:
偏微分大作业一维热传导方程问题——运用隐式格式求解数值解目录TOC\o1-3\h\z\u问题描述31解析解——分离变量法32数值解——隐式格式53证明隐式格式的相容性与稳定性54数值解——分析与__tlab实现65数值解与解析解的比较96随时间变化的细杆上的温度分布情况117稳定后细杆上的温度分布情况12____13附录14有限长杆的一维热传导问题问题描述一根单位长度的细杆放入100℃的沸水中,当细杆的温度达到100℃时取出假设细杆四周绝热;在时间t=0时,细杆两端浸入0℃的冰水中一维热传导方程,现在令,从而可知本题现在要求细杆温度分布1解析解——分离变量法热传导偏微分方程1其中,首先令2将2式带入1式得于是可得可以得到两个微分方程先求解空间项当时,由于可知:由于解的收敛性,则此时是平庸解当时,则此时是平庸解当时,,其中所以,,因为所以,,则,初始条件最终,,2数值解——隐式格式目前,研究热传导问题特别是非稳态热传导问题十分重要这里使用隐式格式利用,关于t进行向前差商;关于x进行二阶中心差;代入偏微分方程可以得到隐式差分格式
(1)3证明隐式格式的相容性与稳定性
(1)相容性代入隐式格式得2将
(2)与原微分方程相减,得到截断误差所以此隐式格式与原微分方程相容
(2)稳定性令网格比为,则可以将
(1)式改写得到
(3)首先令
(4)将
(4)代入
(3)式,根据欧拉公式化简得
(5)故得放大因子是所以根据Fourier方法,隐式格式恒稳定4数值解——分析与__tlab实现
(1)边值与初值离散化将边值与初值离散化,与式
(3)联立得差分线性方程组,,再将方程组改写成的形式本题的边界条件均为零所以可以将上式改写
(2)__tlab的实现杆长1米,时间2秒设计空间步长h=
0.1和时间步长t=
0.01,网格比是从而得到划分的空间网格点数是M1+1,时间网格点数是M2+1先设初始的温度矩阵UM2+1M1+1再将边界条件和初始条件编写到表示温度分布的矩阵中具体代码可见最后附录编写矩阵A核心代码对角线Aii=1+2r对角线的右方和__Aii+1=-r;Ai+1i=-r;下面就要运用进行迭代当k=1时,A*U2j=U1j当k=2时,A*U3j=U2j当k=3时,A*U4j=U3j以此迭代下去直到k=M2就可以得到整个温度随时间和空间的分布矩阵U数值解画图,如图1a和图1b所示图1a数值解的温度分布图现在将着色平稳过渡图1b着色平稳过渡的数值解的温度分布图5数值解与解析解的比较首先,我们需要将解析解离散化,解析解中有一项,当n越来越大时,会快速趋于0,故我们可以取n=8000现在来证明可行性,在__tlab里的工作空间运算将解析解的温度分布画出来,数值解画图,如图2所示图2解析解的温度分布图将数值解与解析解相减,得到误差图如图3a和图3b,我们从图3a上可以看出空间上的误差,在边界处误差比较大图3(a)数值解与解析解空间误差我们从图3a上可以看出时间的误差,在时间的最开始,处误差最大,然后又有一个小的波动,最后就误差渐渐变小,最后趋于0图3(b)数值解与解析解时间误差6随时间变化的细杆上的温度分布情况从数值解的温度分布三维图,如图4a和图4b可以看出随着时间的增加,细杆温度下降最后趋于0℃从物理角度来说细杆的温度会不断地向两端扩散,热量会慢慢散失,最终随着时间的增加,细杆的温度会趋于0℃图4a细杆温度随时间的变化图现取细杆中心处一点,观看它随时间的温度变化情况图4b细杆__(x=
0.5)温度随时间的变化图7稳定后细杆上的温度分布情况从图像上可以看出,最后稳定的情况下,细杆的温度是0℃____
[1]冯立伟.热传导方程几种差分格式的__TLAB数值解法的比较[J].沈阳化工大学,辽宁沈阳.20116.
[2]一维热传导方程数值解法及__tlab实现[EB/OL].2014-11-20http://___.doc
88.com/p-65221__
1945.html附录代码%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%此程序用于解决一维热传导方程:ut-a^2uxx=0%%边界条件u0t=uLt=0%%初始条件ux0=100x!=0和L%%u00=0%%uL0=0%%其中a^2=1L=1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;clearall;%区域及划分网格L=1;%单位长度的细杆T=2;%时间h=
0.1;%%%%空间的划分%%%%t=
0.01;%%%%时间的划分%%%%r=t/h*h;%网格比%设计步长M1=L/h;M2=T/t;%构造边界条件%构造的矩阵U时间,空间)U=zerosM2+1M1+1;%编程包含边值,如Uk1=u0tfork=1:M2+1%时间划分了M2份,有M2+1个节点Uk1=0;%两个边界处温度恒为零UkM1+1=0;end;%构造初始条件forj=2:M1%位置划分了M1份,有M1+1个节点U1j=100;end;U11=0;U1M1+1=0;%差分格式的矩阵形式A*Uk+1j=Ukj%构造矩阵AA=zerosM1-1;fori=1:M1-1Aii=1+2*r;end;fori=1:M1-2Aii+1=-r;Ai+1i=-r;end;%构造AU=B中的B%本题边值的特殊,矩阵B大大简化了B=zerosM1-11;fork=1:M2j=2:M1;Bj-11=Ukj;x=A\B;forj=2:M1Uk+1j=xj-1;%k+1时刻的不同位置的温度分布end;end;%作图x=0:h:1;y=0:t:2;[xxyy]=meshgridxy;figure1;surfxxyyU;shadingflattitle一维热传导方程--数值解--温度分布图;xlabel位置x;ylabel时间t;zlabel温度T;figure2s=0;fori=1:8000s=s+200*1--1^i/i*pi*sini*pi*xx.*exp-i^2*pi^2*yy;end;surfxxyys;title一维热传导方程--解析解--温度分布图;xlabel位置x;ylabel时间t;zlabel温度T;figure3x=0:h:1;y=0:t:2;[xxyy]=meshgridxy;dd=U-s;surfxxyydd;title一维热传导方程--误差--温度分布图;xlabel位置x;ylabel时间t;zlabel误差(数值解减解析解);figure4z=zerosM2+11;ifmodM1+12~=0i=1:M2+1;zi1=UiM1/2;elsei=1:M2+1;zi1=UiM1+1/2;end;plotz;title温度随时间增加的趋势图;xlabel时间t;ylabel温度T;。