物理课件第四章导热问题的数值求解.docx
《物理课件第四章导热问题的数值求解.docx》由会员分享，可在线阅读，更多相关《物理课件第四章导热问题的数值求解.docx（21页珍藏版）》请在冰豆网上搜索。
物理课件第四章导热问题的数值求解
第四章导热问题的数值求解
随着计算机的普及应用和性能的不断改善，以及相关的数值计算方法的发展和应用程序的开发，传热学数值计算方法作为数值求解传热问题的有效工具也得到了相应的发展，利用计算机求解传热学问题愈来愈受到人们的普遍重视，而且在计算复杂传热问题中显示出它的优越性，因而成为传热学的一个重要的分支。
数值传热的相关内容也很自然地成为工程类学生学习传热学课程的不可缺少的部分。
为了使学生能简要地掌握传热学数值计算的基本方法，在这里我们以导热问题为例对传热学数值计算方法做一个简单的介绍。
4-1导热问题数值解概述
在第二章和第三章中我们对较为简单的导热问题，如一维、二维简单几何形状和边界条件的稳态导热和非稳态导热、以及一些特殊导热问题，象通过肋片的导热和忽略内阻的集总导热系统，进行了分析求解。
然而对于一些更为复杂的导热问题，如复杂的几何形状和边界条件以及物性变化较大的情况，分析求解往往很复杂或者根本不可能。
此时求解问题的唯一途径是利用数值分析的办法获得数值解。
数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再求解。
这里要介绍的是后一种方法。
如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。
作为一本入门的教材，这里仅向读者简要地介绍用有限差分析方法从微分方程确立代数方程的处理过程。
有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限个离散点上的数值集合来近似表示。
有限差分的数学基础是用差商代替微商（导数），而几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。
在图4-1中可以看出有限差分表示的温度场与真实温度场的区别。
图中用T0、T1、T2…表示连续的温度场T；Δx为步长，它将区域的x方向划分为有限个数的区域，Δx0、Δx1、Δx2…，它们可以相等，也可以不相等。
当Δx相等时，T1处的真实变化率a可以用平均变化率b、c或d来表示，其中b、c和d分别表示三种不同差分格式下的温度随时间的变化率，即：
b为向后差分格式
；
c为向前差分格式
；
d为中心差分格式
。
这种差分格式也可以推广到高阶微商的情形。
对于二阶微商的差分格式可以在一阶差分格式的基础上得出：
。
采用这样的处理之后，反映温度场随时间、空间连续变化的微分方程就可以用反映离散点间温度线性变化规律的代数方程来表示。
当利用相应的数学办法求解这些代数方程组之后，我们就能获得离散点上的温度值。
这些温度值就可以近似表示温度场的连续的温度分布。
从上面的分析不难看出，当我们要对导热问题进行数值求解时一定要采取三个大的步骤，即研究区域的离散化；离散点（节点）差分方程的建立；节点方程（代数方程）的求解。
下面我们将导热问题的数值求解进行较为详细的讨论。
4-2研究区域温度场的离散化
导热问题的温度场是假设为时间和空间的连续函数，当进行数值求解时首先要做的事情是在所研究的时间和空间区域内把时间和空间分割成为有限大小的小区域，尤如地球被人为地划分为不同的地域且冠以不同的名称，时间被年、月、日和时、分、秒分割。
如果在所分割的每一个时间间隔和空间区域内均用同一个温度值来表示，那么原来连续变化的温度场就被一个离散的阶跃变化的温度分布所代替。
这就是连续变化的温度场离散化处理的基本思路。
这里我们以一个矩形长柱体的非稳态导热过程为例来讨论区域离散化问题。
如果不考虑矩形长柱体长度方向上的温度变化，那么它是一个二维非稳态导热问题，图4-2表示了长柱体矩形截面上区域离散化的情况。
图中可见，对于给定的空间区域，在x方向上的步长为Δx，在y方向上的步长为Δy，用它们作为空间尺度可以将矩形区域划分成纵横交错的网格，交点称为节点。
然后以节点为中心，在两个节点的中心处划分界限，定出节点的控制面积，对于三维情况则为控制体积或控制容积，因而常在一般意义上称之为节点的控制体。
控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。
网格的步长在每一个方向上可以均匀划分，也可以不均匀的划分。
因此，选用不同的步长和不同的划分方法，可以将同一区域划分出不同大小、不同数目的控制区域，以及不同数目的节点数。
显然，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散的温度场就会更加接近连续的温度场，但计算工作量也会随之增加。
在时间方向上离散化的步长常用Δτ来表示，Δτ的选取也是可大可小的，也可以随时间的进程而变化。
显然，无限小的时间步长Δτ亦会使得离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量的增加。
4-3温度场节点方程的建立
为了得出所研究区域的节点温度，必须建立相应的节点方程。
建立节点方程可以采用不同的方法，为了更好地理解节点方程的物理意义和掌握节点方程的建立方法，我们采用控制体热平衡法来建立节点方程。
下面我们以实例的形式介绍不同节点的节点方程的建立过程。
1控制体的内节点方程
控制体热平衡法建立节点方程的过程是将能量守恒方程应用于控制体，建立该节点与周围节点之间的能量平衡关系式，再利用傅立叶的导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。
考察图4-2中的节点P及其控制体，由能量平衡关系应有
，4-1
式中，QW、QE、QS和QN分别为邻近节点W、E、S和N通过传导方式传给节点P的热流量；QV为单位时间控制体内热源的发热量；ΔΕ为控制体单位时间内热能的增加量。
由导热傅立叶定律，在线性温度分布的假设下，时刻K周围节点传给节点P的热流量分别为：
以及控制体的发热流量
（qV为内热源强度，即单位时间单位体积的内热源发热量。
）
控制体单位时间的内能增加量为
，前者为时间上的向前差分，而后者为时间上向后差分。
以上关系式中温度T的上标为所在时刻，下标为所在空间位置。
将以上关系式一并代入方程4-1中，且假设Δx=Δy，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即
显示差分格式
4-2
和隐示差分格式
。
4-3
比较上面两种差分格式可以看出，显示差分格式最突出的优点是节点温度表达式的右边只涉及K时刻的节点温度值，那么只要知道这一时刻周围节点的温度值就可以求出该节点的下一时刻的温度值；而隐示差分格式却相反，温度表达式的两端都是同一K时刻的节点温度值，这就意味着必须同时计算同一时刻所有节点的温度值，即必须联立求解K时刻所有节点的差分方程组，增大计算工作量是显而易见的。
虽然显示差分格式计算比较方便，但它却存在着一个缺点，即计算式中aΔτ/Δx2值必须满足一定的条件才不至于引起数值计算出现不收敛的问题，这在数值计算中称为差分格式的不稳定性。
这里差分方程稳定性的条件是方程4-2中的变量T前面的系数必须大于或等于零，分析一下差分方程中的各项系数，在TPK前的系数应为
，改写成为
。
4-4
此式称为显示差分格式的稳定性判据，从中看出时间步长和空间步长是相互制约的。
为了获得较为精确的节点温度值，空间步长Δx的选择不能太小，按照稳定性判据的要求势必会使时间步长Δτ也要相应地不能太大，因而必须在增加节点数目的同时增多时间间隔，从而使计算工作量加大。
与显示差分格式相反，由于隐示差分格式的节点方程中没有会使方程系数成为负值的系数项，因而不存在方程求解的不稳定性的问题。
也就是说，对于隐示差分格式，无论aΔτ/Δx2中的Δx和Δτ取什么样的数值，均不会出现数值计算结果的不收敛问题，因而是无条件稳定的。
这样就使得我们能在满足一定精度的情况下尽可能地加大时间步长或空间步长，亦可以在计算过程中随意改变步长，而不必担心会造成计算结果的不收敛。
2控制体的边界节点方程
在数值计算中所研究的区域的边界条件是通过边界节点的节点方程来反映的，因而边界节点的差分方程的建立十分重要。
这里同样采用边界节点的控制体热平衡来确立边界节点的差分方程。
下面将具体讨论一些典型的边界节点的节点方程。
对流换热边界条件下的节点方程。
如图4-3所示的边界节点P，其控制体的热平衡关系式为
，4-5
式中各项可以写成：
将上述关系式一并代入4-5式中，经整理可得出在二维非稳态导热过程中处于平直边界上的边界节点在对流换热边界条件下的节点差分方程（当Δx=Δy），即：
对于显示差分格式有
；4-6
对于隐示差分格式有
。
4-7
分析公式4-6，可以得到其稳定性判据为
，不难看出条件要比内节点的情形更加苛刻一些。
绝热边界条件下的边界节点方程。
如图4-4所示的绝热边界节点P，其控制体的热平衡关系式为
，4-8
式中各项可写成：
将上述关系式代入4-8式中，可以整理得到绝热边界条件下的节点方程（当Δx=Δy），即：
对于显示差分格式有
；4-9
对于隐示差分格式有
。
4-10
对于其它边界节点也可以采用上述的控制体热平衡的办法建立对应的节点差分方程式。
从以上的节点方程中可以发现两个特别的系数项，即
它们分别是控制体（元体）的傅立叶数和毕欧数，其物理意义是前者表征控制体的导热性能与热储蓄性能的对比关系，反映控制体温度随时间变化的动态特性；而后者则体现控制体和环境间的换热性能与其导热性能的对比关系。
另外，还要一个常数项
它则反映了内热源产生的热流引起的节点温度升高值。
下面我们归纳性地给出二维非稳态导热问题的两种差分格式的节点方程式，见表4-1。
表4-1二维非稳态导热节点差分方程（式中Fo=aΔτ/Δx2,Bi=αΔx/λ）
节点内型
差分方程（Δx=Δy）
稳定性条件
内节点
显示差分格式
Ti,jK+1=Fo（Ti+1K,j+Ti-1K,j+TiK,j-1+Ti,Kj+1）+（1-4Fo）Ti,jK
Fo≤1/4
隐示差分格式
Ti,jK=[Fo（Ti+1K,j+Ti-1K,j+TiK,j-1+Ti,Kj+1）+Ti,jK-1]/（1+4Fo）
无条件
平直绝热边界节点
显示差分格式
Ti,jK+1=Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+（1-4Fo）Ti,jK
Fo≤1/4
隐示差分格式
Ti,jK=[Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+Ti,jK-1]/（1+4Fo）
无条件
平直对流边界节点
显示差分格式
Ti,jK+1=Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+
2FoBiT∞K+（1-2FoBi-4Fo）Ti,jK
Fo≤
1/[2（2+Bi）]
隐示差分格式
Ti,jK=[Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+2FoBiT∞K
+Ti,jK-1]/（1+2FoBi+4Fo）
无条件
平直辐射边界节点
显示差分格式
Ti,jK+1=Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+（1-4Fo）Ti,jK-
2σ0εS[（Ti,jK）4-（T∞K）4]/（ρcΔx）
隐示差分格式
Ti,jK（1-4Fo）+2σ0εS（Ti,jK）/（ρcΔx）=
Fo（2Ti-1K,j+TiK,j-1+Ti,Kj+1）+Ti,jK-1+2σ0εS（T∞K）4/（ρcΔx）
无条件
4-4节点方程的求解
由上面的讨论可以看出，对应于离散温度场的每一个节点均可以列出相应的差分方程，这样就可以得出与节点数目相同的一组代数方程组。
当联立求解这个代数方程组时，最后就可以得出每一个节点的温度值。
一般情况下，差分方程组是线性代数方程组，而线性代数方程组是可以用直接法和迭代法求解的。
常用的直接法有高斯消元法、列主元素消去法和矩阵求逆法，而迭代法常用的有高斯－赛德尔迭代和超（欠）松弛迭代。
下面我们仅介绍迭代法求解代数方程组的过程，并在例题中加以应用。
今有一线性代数方程组：
迭代求解该方程组的思路为，寻找一个由（T1，T2，…，Tn）组成的列向量，使其收敛于某一个极限向量（T1*,T2*,…，Tn*）,且该极限向量就是该方程的精确解。
当这个线性代数方程组的系数项aii≠0（i=1,2,…,n）时，可将其改写成迭代形式，有：
以上各式可以用一个通用的形式来表示：
i=1,2,…,n
利用上式就可以进行迭代求解了，其步骤是，合理选择（假设）各节点的初始温度，将其作为第零次迭代的近似温度值，记为Ti（0）（i=1,2,…,n）;将Ti（0）代入上式的右端，得到第一次迭代的近似值Ti
（1）；之后将Ti
（1）再代入上式的右端，则得出第二次的近似值Ti
（2）；如此反复进行下去，直至进行到K次，使相邻的两次近似解Ti（K+!
）和Ti（K）（i=1,2,…,n）之间的偏差小于预先设定的小量ε时，即满足∣Ti（K+1）－Ti（K）∣≤ε（i=1,2,…,n）或∣（Ti（K+1）－Ti（K））/Ti（K）∣≤ε（i=1,2,…,n）。
此时各节点的温度值[T1（K）,T2（K）,…,Tn（K）]已经有足够的精度用来表示代数方程组的解，从而可以结束方程求解的迭代过程。
从上述的迭代过程不难发现，当我们用第零次迭代值去进行第一次迭代时，Ti
（1）的值已经不断地产生出来，当计算Tr
（1）时，到r-1的Ti
（1）已经求出。
如果此时在计算Tr
（1）时涉及到的Ti（0）（i=1,2,…,r-1）全部用已求出的Ti
（1）（i=1,2,…,r-1）代替，这势必会加快迭代收敛的速度。
这种改进后的迭代方法被称之为高斯－赛德尔迭代法。
高斯－赛德尔迭代法的方程组迭代形式为：
i=1,2,…,n。
归纳起来，高斯－赛德尔迭代法的求解步骤可表述为，将代数方程组写成迭代形式；设初始值经迭代得出节点新值；有新值则去掉旧值，不断以新换旧，且在迭代过程中应用；在迭代获得满足给定精度的节点温度值后结束方程组的迭代。
4-5导热问题数值求解举例
为了使读者对导热问题的数值求解过程有一个完整的认识，下面将以无限大平板对称加热的非稳态导热过程为例，阐明显示差分格式和隐示差分格式的计算求解过程。
有一厚度为0.12m的无限大平板，初始温度为20℃,两侧表面同时受到温度为150℃的流体加热，流体与平板表面之间的对流换热系数为24W/（m2℃）,平板材料的导热系数为0.24W/（m℃）,而热扩散系数为0.147×10-6m2/S。
试计算平板温度分布随时间的变化。
由于大平板是对称受热，利用温度场的对称性可以只对平板的一半进行数值计算。
于是就可以给出该导热问题的边界条件为一侧是平板与环境进行对流换热的第三类边界条件，另一侧（即对称中心面）是处于绝热状态的第二类边界条件。
下面我们分别用显示差分格式和隐示差分格式对其进行数值分析求解。
1显示差分格式求解过程介绍
针对该问题可从表4-1中简化出中心节点和边界节点的差分方程，即
相应的稳定性条件在目前情况下有两个，一个是内节点和绝热边界节点应满足的F0=aΔτ/Δx2≤1/2；一个是对流边界节点应满足的F0≤1/（2+Bi）。
在这两个条件中后一个要求要苛刻一些，故只需满足后一个稳定性条件即可。
由稳定性条件可见，空间步长Δx和时间步长Δτ是相互关联的，当选定Δx之后Δτ也随之被确定下来。
在此我们选取Δx=0.006m，从稳定性条件计算出Δτ≤76.3S,计算中选取Δτ=50S。
图4-6显示差分格式计算机程序框图
计算程序的框图如图4-6所示，而程序中的变量标识符则说明于表4-2中。
表4-2计算机程序中的标识符（含显示和隐示差分格式用到的变量和常数）
A
平板材料的导温系数
K
平板材料的导热系数
AB（I,J）
节点温度中间变量
I
表示空间的下标变量
AA
计算的累计时间
IJ
表示空间的下标变量
BIO
网格毕欧数
J
表示时间的下标变量
B1
中间变量
L
平板厚度的一半
B2
中间变量
M
计算的时间间隔数目
DT
时间步长
N
计算的空间间隔数目
DX
空间步长
T（I,J）
节点温度变量
FO
网格傅立叶数
TF
环境流体温度
H
对流换热系数
T0
平板初始温度
计算机程序采用FORTRAN语言编写，详细清单如下：
REALK，L
DIMENSIONT（50,50）
READ（*,*）H,K,A,DT,TF,L,T0,N,M
DX=L/FLOAT（N）
BIO=（H*DX）/K
FO=（A*DT）/（DX**2）
IF（FO.GT.（1.0/（2.0+2.0*BIO）））GOTO60
DO10I=1,N+1
10T（I,1）=T0
DO20J=2,M+1
T（1,J）=2.0*FO*T（2,J-1）+（1.0-2.0*FO）*T（1,J-1）
DO30IJ=2,N
T（IJ,J）=FO*（T（IJ-1,J-1）+T（IJ+1,J-1））
30T（IJ,J）=T（IJ,J）+（1.0-2.0*FO）*T（IJ,J-1）
T（N+1,J）=2.0*FO*（T（N,J-1）+TF*BIO）
T（N+1,J）=T（N+1,J）+（1.0-2.0*FO-2.9*FO*BIO*T（N+1,J-1）
AB=FLOAT（J-1）*DT
WRITE（*,40）AB
40FORMAT（F9.2）
WRITE（*,50（T（II,J）,II=1,N+1）
50FORMAT（11F7.2）
20CONTINUE
60STOP
END
利用以上程序经过49次运算，即计算到非稳态导热过程进行2400S后的结果，将其显示于表4-3中表中x为平板的坐标位置，x=0为对称中心平面，x=0.06m为平板的外表面。
从表中可见经过2400S后，温度为150℃的流体仍然在向平板加热。
进一步的计算表明，平板被均匀加热到150℃的时间大约到96000S（26.7h），共需计算1920次。
表4-3算例进行了49次迭代后所得的结果
x[m]
0.0
0.006
0.012
0.018
0.024
0.030
0.036
0.042
0.048
0.056
0.06
T[℃]
23.02
23.58
25.25
28.34
33.17
40.16
49.69
62.05
77.29
95.21
115.31
2隐示差分格式计算机程序的编制
从表4-1中同样可以简化出该算例的隐示差分格式的内节点及边界节点的差分方程式如下：
TiK=[Fo（Ti+1K+Ti-1K）+Ti-1K-1]/（1+2Fo）（内节点）；
T1K=2FoT2K+T1K-1）/（1+2Fo）（对称中心边界节点）；
TN+1K=（2FoTNK+2FoBiT∞K+TN+1K-1）/（1+2FoBi+2Fo）（对流换热边界节点）。
由于隐示差分格式是无条件稳定的，故在编制程序时不必设立稳定性判别语句，但确要联立求解所有的节点方程组。
这里我们采用高斯－赛德尔迭代法，程序用FORTRAN语言编写。
程序中的EPS是精度控制的允许误差，计算中取为0.01。
若时间步长也取Δτ＝50S，计算到2000S时可以得到与表4-3所示的相近似的结果；若取Δτ＝2000S，只需计算48次即可接近平板被均匀加热到150℃的稳定工况。
当然随着Δτ的增大，计算结果的误差也会相应增大。
计算程序框图如图4-7所示。
按上面的框图编写的计算机程序如下：
REALK,L
DIMENSIONT（50,50）,AB（50,50）
READ（*,*）H,K,A,DT,DF,L,T0,EPS,N,M
DX=L/FLOAT（N）
BIO=（H*DX）/K
FO=（A*DT）/DX**2
B1=1.0+2.0*FO
B2=1.0+2.0*FO+2.0*FO*BIO
DO10I=1,N+1
10T（I,1）=T0
DO20J=2,M+1
DO30IJ=1,N+1
T（IJ,J）=T0
30AB（IJ,J）=T0
35T（1,J）=（T（1,J-1+2.0*FO*T（2,J））/B1
DO40IK=2,N
40T（IK,J）=（T（N+1,J-1）+FO*（T（IK+1,J）+T（IK-1,J）））/B1
T（N+1,J）=（T（N+1,J-1）+2.0*FO*BIO*TF+2.0*FO*T（N,J））/B2
DO50II=1,N+1
IF（ABS（T（II,J）-AB（II,J））•••GT•EPS）GOTO65
50CONTINUE
GOTO75
65DO70IM=1.N=1
70AB（IM,J）=T（IM,J）
GOTO35
75AA=FLOAT（J-1）*DT
WRITE（*,80）
80FORMAT（F9.2）
WRITE（*,90）（T（IL,J）,IL=1,N+1）
90FORMAT（11F7.2）
20CONTINUE
60STOP
END
思考题
1.试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。
2.试说明用节点控制体热平街法建立节点温度差分方程的基本思想。
3.推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散的差分方程的过程十分相似，什么前者得到的是温度场的精确描写，而由后者解出的却是温度场的近似分布。
4.第三类边界条件的边界节点也可以采用将第三类边界条件表达式中的一阶导数用差分公式表示来建立其节点差分方程。
试比较这样建立起来的差分方程与用热平衡法建立起来的差分方程的异同与优劣。
5.什么是显式格式，而什么是隐示差分格式？
为什么显式格式在计算中存在稳定性问题而隐示差分格式却不存在?
6.用高斯-赛德尔迭代法求解代数方程组时是否一定可以得到收敛的解?
在不能得出收敛的解时是否是因为初场的假设不合适造成的?
习题
一般性数值计算
4-1试用数值计算证实，对方程组
用高斯-赛德尔迭代法求解，其结果是发散的，并分析其原因。
4-2试对附图所示的常物性、无内热源的三维稳态导热问题用高斯-赛德尔迭代法计算t1，t2，t3，t4之值。
4-3试对附图所示的等截面直肋的稳态导热问题，用数值方法求解节点2、3的温度。
图中t0=85℃、tf=25℃、h=30W/（m2·K）。
肋高H=4cm，纵剖面面积AL=4cm2，导热系数λ=20W/（m·K）。
离散方程建立
4-4试将直角坐标中的常物性、无内热源的三维非稳态导热微分方程化为显式差分格式，并指出其稳定性条件（Δx≠Δy）。
4-5极坐标中常物性、无内热源的非稳态导热方程为
试利用本题附图中的符号，列出节点（i，j）的差分方程式。
4-6一金属短圆柱在炉内受热后被竖直地移置到空气中冷却，底面可认为是绝热的。
为用数值法确定冷却过程中柱体温度的变化，取中心角为1rad的区域来研究（如本题附图所示）。
已知柱体表面发射率ε、自然对流表面传热系数从环境温度t∞、金属的热扩散率a,试列出图中节点（1，1）、（M,l）、（M，n）及（M,N）的离散方程式。
在r及z方向上网格是各自均分的。
4-7一个三维物体的竖直表面受液体自然对流冷却，为考虑局部表面传热系数的影响，表面传热系数采用h=c（t-tf）1.25来表示。
试列出附图所示的稳态、无内热源物体边界节点（M,n）的温度方程，并对如何求解这一方程组提出你
升级会员 