《计算物理(本科)》[第9章].ppt

《计算物理(本科)》[第9章].ppt

《《计算物理(本科)》[第9章].ppt》由会员分享，可在线阅读，更多相关《《计算物理(本科)》[第9章].ppt（26页珍藏版）》请在冰豆网上搜索。

第第九九章章热热传传导导方方程程的的数数值值解解法法第九章第九章热传导方程的数值解法热传导方程的数值解法9.1热传导方程的引入热传导方程的引入物理学中对热传导现象和扩散现象等物理过程的描述，物理学中对热传导现象和扩散现象等物理过程的描述，都归结为的二阶偏微分方程，统称为都归结为的二阶偏微分方程，统称为热传导方程热传导方程，也，也称为称为抛物线型方程抛物线型方程。

本章以热传导现象进行介绍：

本章以热传导现象进行介绍：

当物体内部各处温度不同时，就会有热量从温度较高当物体内部各处温度不同时，就会有热量从温度较高处向温度较低处传递，这就是处向温度较低处传递，这就是热传导现象热传导现象。

合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法一、传热过程一、传热过程：

在：

在t时刻，时刻，dt时间内通过时间内通过dS截面的热量截面的热量其中其中K为热传导系数；为热传导系数；是在是在S面的法向上温度梯度。

面的法向上温度梯度。

对对上式上式进行闭曲面积分，并用矢量积分定理得进行闭曲面积分，并用矢量积分定理得其中其中V为为闭曲面包括的体积；闭曲面包括的体积；t1,t2表示表示t2t1时间段的时间段的积分；积分；合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法三、闭体能量守恒三、闭体能量守恒：

若体积：

若体积v内无内无热源，热传导系数热源，热传导系数K是是常数时，根据能量守恒常数时，根据能量守恒Q=Q得三维齐次热传导方程得三维齐次热传导方程。

二、体内耗热过程二、体内耗热过程：

设物体的比热容为：

设物体的比热容为c，密度为密度为，则，则v内温度变化所消耗的热量为内温度变化所消耗的热量为合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法9.2一维热传导方程的差分法一维热传导方程的差分法一、热传导方程：

一、热传导方程：

物体内无热源的一维热传导方程为物体内无热源的一维热传导方程为二、初始条件（二、初始条件（t=0的温度分布情况）：

的温度分布情况）：

三、边界条件（物体边界的温度分布情况）：

三、边界条件（物体边界的温度分布情况）：

第一类边界条件第一类边界条件;第二类边界条件第二类边界条件合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法第三类边界条件第三类边界条件其中其中h1（t）、h2（t）、g1（t）、g2（t）为给定的函数，为给定的函数，h1（t）0、h2（t）0且不同时为且不同时为0。

合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法并使并使Nh=lM=T则有则有t=kk=0,1,2,M（时间步序号时间步序号）x=ihi=0,1,2,N（空间空间x方向的步序号方向的步序号）四、差分方法求解四、差分方法求解建立差分格式建立差分格式通常分别取空间步长和时间步长为通常分别取空间步长和时间步长为h和和，均为常数；计均为常数；计算时的步序号空间用算时的步序号空间用i表示，时间用表示，时间用k表示。

表示。

合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法定义：

一阶向前差商近似定义：

一阶向前差商近似一阶向后差商近似一阶向后差商近似二阶二阶中心差商近似中心差商近似用二用二阶阶中心差商作为二中心差商作为二阶微商近似是阶微商近似是对时间的一阶差商对时间的一阶差商近似是近似是合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法由上两式，并由上两式，并令令=/h2，得得一维差分格式一维差分格式其中其中N=l/h,M=T/,“”表示取整。

表示取整。

定解条件定解条件几何关系几何关系把某一时刻计算的各点称为一层，计算过程是一层层把某一时刻计算的各点称为一层，计算过程是一层层地进行的，各点关系如图。

由图可见，地进行的，各点关系如图。

由图可见，k+1时刻时刻i点的点的值是由值是由k时刻的时刻的i-1、i和和i+1三点的值推算出来的。

三点的值推算出来的。

txk+1ki+1ii-1i0,N为什么?

合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法可以证明，上可以证明，上式的差分式的差分格式收敛稳定条件为格式收敛稳定条件为=/h21/2五、归纳上述一维热传导的差分格式为五、归纳上述一维热传导的差分格式为合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法六、差分格式计算步骤六、差分格式计算步骤1.给定给定，l，h，T；2.计算计算N=l/h，M=T/，=h2/；计算边界值：

计算边界值：

3.计算初始值：

计算初始值：

4.差分格式的计算差分格式的计算合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法9.3二维热传导方程的差分法二维热传导方程的差分法二维热传导方程求解与一维的类似，在确定了差分格二维热传导方程求解与一维的类似，在确定了差分格式和给出定解条件后，按时间序号分层计算，每一层式和给出定解条件后，按时间序号分层计算，每一层由二维点阵组成，通常称为由二维点阵组成，通常称为网格网格。

一、热传导方程：

一、热传导方程：

均匀物体内部无热源的二维热传导方程均匀物体内部无热源的二维热传导方程为为二、初始条件二、初始条件合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法建立差分格式建立差分格式设空间步长为设空间步长为h，时间步长为时间步长为。

如图所示，将如图所示，将xoy平面均分成平面均分成NM的网格，并使的网格，并使Nh=lMh=s三、差分方法求解三、差分方法求解则有则有t=kk=0,1,2,（时间时间t序号序号）x=ihi=0,1,N（空间空间x序号序号）y=jhj=0,1,M（空间空间y序号序号）yxMNM2M1合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法对节点对节点（i，j），在，在k时刻时刻（即即t=k）的差分式的差分式将上式代入二维热传导方程中，并将上式代入二维热传导方程中，并令令=/h2即可得即可得二维二维热传导的差分格式热传导的差分格式i=1,2,N-1；j=1,2,M-1。

（注意：

注意：

i0,N；j0,M）合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法几何关系几何关系用上式的二维热传导的差分格式，计算用上式的二维热传导的差分格式，计算k+1时刻点时刻点（i,j）处的温度值时，使用了处的温度值时，使用了k时刻点时刻点（i,j）及其相邻的及其相邻的四个点的温度值，如图所示相关各点之间的关系。

四个点的温度值，如图所示相关各点之间的关系。

txk+1ki+1ii-1j+1j-1j合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法四、边界条件四、边界条件（如图所示的具体问题如图所示的具体问题）绝热边界：

粉红色部分绝热壁，绝热边界：

粉红色部分绝热壁，x=0,y（0,M1h）.and.（M2h,Mh）;x=Nh,y（0,Mh）yxMNM2M1上式上式改变为差分形式，整理得改变为差分形式，整理得绝热位置应满足绝热位置应满足j0,M此处由下面此处由下面边界条件定边界条件定合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法恒温热源边界：

红色部分，即恒温热源边界：

红色部分，即x=0,yM1h，M2h归一化后高温源温度取归一化后高温源温度取“1”，即，即温度为温度为0的边界的边界：

实际也是与恒温源相连，白色部分，：

实际也是与恒温源相连，白色部分，即即y=0和和y=Mh两个边。

两个边。

归一化后低温源温度取归一化后低温源温度取“0”，即，即合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法可以证明，上可以证明，上式的差分式的差分格式收敛稳定条件为格式收敛稳定条件为=/h21/4五、归纳上述二维热传导的差分格式五、归纳上述二维热传导的差分格式合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法六、差分格式计算步骤六、差分格式计算步骤1.给定给定,h,T及及XN和和YM；2.计算计算N=XN/h，M=YM/h，=h2/，k的上界的上界T/；3.计算计算xi=ih,yi=jh,tk=k（i=0,1,，N；j=0,1,M；k=0,1,T/）；4.计算初始值和边界值；计算初始值和边界值；5.差分格式的计算差分格式的计算ui,j,k+1（i=1,2,N-1；j=1,2,M-1；k=0,1,T/）合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法9.3MATLAB实现热传导方程数值解法实现热传导方程数值解法Matlab提供了一个求解此热传导方程函数，形式为提供了一个求解此热传导方程函数，形式为u=parabolic（u0,tlist,b,p,e,t,c,a,f,d）其中其中u0为初始值；为初始值；tlist为时间序列；为时间序列；b是边界条件，是边界条件，可以用矩阵形式，也可以用可以用矩阵形式，也可以用m文件格式，主要依赖时文件格式，主要依赖时间间t；p,e,t为网格参数；为网格参数；c,a,f,d是方程系数，可以是时是方程系数，可以是时间的函数；另外可以在函数末加入误差限来控制。

对间的函数；另外可以在函数末加入误差限来控制。

对于标量形式的热传导方程，函数返回值为一个矩阵。

于标量形式的热传导方程，函数返回值为一个矩阵。

设热传导方程形式为设热传导方程形式为合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法解：

在解：

在MATLAB中建立中建立m文件文件jswlx_9_3_1.m【例例1】求热传导方程求热传导方程求解的求解的范围为正方形区域：

范围为正方形区域：

-1x,y1，t0；初始条件：

当初始条件：

当x2+y20.42时，时，u（x,y,0）=1；其它处为其它处为0。

边界条件：

当边界条件：

当t0时，时，u（1,y,t）=0，u（x,1,t）=0。

t=0t0tT合肥工业大学电子科学与应用物理学院第第九九章章热热传传导导方方程程的的数数值值解解法法%求解求解标准热传导方程标准热传导方程%du/dt-div（grad（u）=0g=squareg;%定义单位方形区域函数定义单位方形区域函数b=squareb1;%定义方形定义方形0边界条件函数边界条件函数c=1;a=0;f=0;d=1;p,e,t=initmesh（g）;%初始的粗糙网格初始的粗糙网格化化P_（X，Y）%定义初条件定义初条件,变量在半径为变量在半径为0.4的的%圆周内值为圆周内值为1,其它处则为其它处则为0u0=zeros（size（p,2）,1）;%size（p,2）取取p列长，列长，zeros是取是取1列列size行的行的0阵阵ix=find（sqrt（p（1,:

）.2+p（2,:

）.2）0.4）;%find（X）是取是取X的序号的序号u0（ix）=ones（size（ix）;%在时间段为在时间段为0,0.1的的20个点上求解个点上求解nframes=20;tlist=linspace（0,0.1,nframes）;%以下求此抛物线型方程问题以下求此抛物线型方程问题u1=parabolic（u0,tlist,b,p,e,t,c,a,f,d）;x=linspace（-1,1,31）;y=x;unused,tn,a2,a3=tri2grid（p,t,u0,x,y）;%矩形网格插值，为了加速绘图矩形网格插值，为了加速绘图%以下动画图形显示以下动画图形显示newplot;%建立新的坐标系建立新的坐标系Mv=moviein（nframes）;%影片帧面的内存预置影片帧面的内存预置umax=max（