层合结构瞬态热传导方程的一种精细解法*

蓝林华，富明慧，刘祚秋

(中山大学应用力学与工程系，广东 广州 510275)

热传导问题一直是热门问题，随着科学技术的发展，更多的领域和方向与之紧密相关，此问题的研究也越来越受到关注。在过去的几十年里，人们进行了大量的研究工作，提出了各种近似求解的方法(有限差分法、边界元法、多重网格法、变分法，无网格法等)。刘继军等[1]针对二维热传导方程问题构造了一种三层显式格式并进行了数值计算。Ochiai[2]采用边界元方法讨论了含热源的稳态热传导问题，同时用区域积分法改进了边界元方法并把它运用到功能梯度材料的热传导问题中。李寿佛等[3]用无网格方法研究了二维、三温热传导方程。陈光南等[4]用变分原理在不规则结构网格上建立热流通量形式的差分格式,研究了二维热传导方程的差分数值模拟。Roy等[5]用多重尺度法计算了双曲热传导问题。Chen等[6]结合直接法和伴随矩阵法得到功能梯度材料稳态和瞬态热传导问题的灵敏度方程并进行了分析，随后利用精细积分法求解了功能梯度材料的瞬态传热问题。精细积分法最初是应用于结构动力分析的方法[7-8]，但由于具有精度高、稳定性好等优点，因此迅速被推广到热传导问题[9-10]。本文针对单层或者多层介质上的二维瞬态热传导方程开展研究，对空间进行离散，并对空间离散后的拉普拉斯差分算子由块三角矩阵转化为块对角形式，结合精细积分算法，建立了求解层合结构瞬态热传导问题的一种有效方法。算例表明，它具有良好的数值稳定性，而且精度高、收敛速度快，便于运用。

1 单层介质的二维瞬态热传导问题

考虑矩形区域上二维瞬态热传导方程的混合问题

(1)

其中φ(x,y,t)是所要求的温度分布函数，a为热扩散系数，Ω为矩形区域内部(lx×ly)，∂Ω为矩形区域边界，g(x,y,t)，φ(x,y)，φ0(x,y)都是已知的连续函数。

用图1所示均匀网格将定解区域离散，阴影部分为一矩形小区域单元(h1×h2)，其中h1=lx/(m+1),h2=ly/(n+1)分别是x和y方向的步长。

图1 二维热传导矩形域的网格划分

(2)

记

式中δij为Kronecker_δ符号。利用上述记号，式(2)可表达为

(4)

其中

(5)

式中In×n是n阶单位矩阵，D是n阶三对角矩阵。

容易求得，矩阵D的特征值λ及相应的特征向量分别为[11]

其中k=1,…,n，η=π/(n+1)。

Q=(q1,q2,…,qn)=

由特征矩阵的性质，有

Q-1CQ=Λ

(6)

由于Q是正交而且对称的，有

Q-1=QT=Q

(7)

由式(6)，式(4)可化为

(8)

其中

(i=1,2,…，m)

(9)

式(8)可记为

(10)

其中

矩阵H具有块三对角形式，这对于建立高效算法是十分有利的。

初始条件也相应地化为

(11)

式(10)的解可表达为

(12)

2 层合结构的瞬态热传导问题

图2所示为一(k+1)层矩形层合结构差分网格示意图，沿x方向均匀离散，y方向分段均匀离散，黑粗线是每两层材料之间的分界线。

图2 层合结构的差分网格示意图

为方便推导，本文假设分界面的热接触良好，没有接触热阻。其导热微分方程(不包括分界线)可分段写为

用与方程(2)类似的符号记法，方程(13)的差分形式可写为

(14)

其中h0=l0/(m+1),hr=lr/(nr-nr-1) 分别是x和y方向的步长。

对于第r条分界线上的某一差分点(xp,ynr)，取如图2b阴影部分的一小单元h0×(hr+hr+1)/2， 根据Fourier定律和能量守恒定理，可以直接写出其差分形式

(15)

记

(16)

式中δij为Kronecker_δ符号。

利用上述记号，式(15)、(16)可合并表达为

(17)

其中

(18)

式中Im是m阶单位方阵，Cr、Cnr、Dr、Dnr均是m阶三对角矩阵。类似于上一节，式(17)的系数矩阵的每一个子块同样可以化为对角矩阵，特别地，当ar=ar+1，hr=hr+1时，式(18)可化为式(5)形式，即复合叠层材料变成单层介质材料。同样，上述方法同样可推广到三维的情况。

3 算 例

图3 是用两种方法得到的温度场等值线分布图，其中图3(a)的在t=1 s后的结果，图3(b)是t=10 s后的结果。此结果中有限元法的时间步长取为Δt=0.005 s，可以看到，当时间步长取得合理时，两种方法的结果吻合得很好。

图3 温度场等值线分布图

表1 中心点在t=0.5 s的温度值

Table 1 The temperature of center point att=0.5 s

时间步长/s温度值有限元法本文方法0.005 99.236 506 912 023 6099.211 297 459 871 740.01 99.261 546 525 305 5599.211 297 459 872 300.02 —99.211 297 459 872 610.1 —99.211 297 459 872 780.5 —99.211 297 459 873 08

表1 给出了中心点在t=0.5 s时的温度值，分别用本文方法和有限元方法进行了计算，两种方法均采用相同的时间、空间网格划分。数字的下划线部分是两种方法在不同时间步长下得到结果的不同有效位数，从结果可以看出，对于不同的时间步长，即使时间步长等于时间t，本文方法的计算结果有长达11位的相同有效位数，而有限元法的结果对时间步长的变化情况比较敏感，特别是当步长增大到一定程度时还会出现数值不稳定的情况(表中的虚横线)。事实上，用有限元法计算二维瞬态热传导方程时其差分格式分显式格式和隐式格式，隐式格式能够保证收敛，但在高维的情况下其计算量很大，相当于显式格式计算量的平方陪，故对于高维的情况一般采用显式格式(本文的有限元法也采用显式格式)，但显式格式为了保证收敛性，时间、空间的网格划分须满足一定条件

表2 有限元法在不同时间步长下的稳定性分析

表2 给出了有限元法在不同时间步长下的稳定性分析情况，其中红色的大于0.5表示该时间步长不满足稳定行的充要条件，显然表2 的分析结果跟表1是一致的。

而对于本文采用的精细积分(PTI)方法，即使时间步长取整个时间t，而精细积分算法中的dt=t/220将是一个很小的量(根据需要还可以取得更小)，它总能够满足上面的稳定性条件。也就是说，本文的求解过程可以无条件地满足算法的稳定性要求。

4 结 论

本文针对单层或者多层介质上的二维瞬态热传导方程，给出了对空间的详细离散过程，并对空间离散后的拉普拉斯差分算子由块三角矩阵转化为块对角形式，随后结合已有的精细积分算法，建立了求解层合结构瞬态热传导方程的一种有效方法。算例表明，本文方法不仅具有良好的精度(相同的有效位数长达11位)，还可以无条件地满足算法的稳定性要求，并且该方法还可以推广到三维的情况，因此可以成功地解决层合结构的瞬态热传导问题，为国内、外同类研究提供一定的参考。

参考文献：

[1]刘继军. 二维热传导方程的三层显式差分格式[J]. 应用数学和力学,2003,24(5):537-544.

[2]OCHIAI Y. Two-dimensional steady heat conduction in functionally gradient materials by triple-reciprocity boundary element method[J]. Eng Analysis with Boundary Elements,2004,28(12):1445-1453.

[3]李寿佛,张瑗,刘玉珍. 热传导方程的一类无网格方法[J]. 计算物理,2007,24(5):573-580.

[4]陈光南,张永慧. 基于变分原理的二维热传导方程差分格式[J]. 计算物理,2002,19(4):299-304.

[5]ROY S,VASUDEVA MURTHY A S,KUDUNATTI R B. A numerical method for the hyperbolic-heat conduction equation based on multiple scale technique[J]. Applied Numerical Mathematics,2009,59:1419-1430.

[6]GU Y X,CHEN B S,ZHANG H W,et al. A sensitivity analysis method for linear and nonlinear transient heat conduction with precise time integration[J]. Struct Multidisc Optim,2002,24:23-37.

[7]钟万勰. 一类非线性热传导方程的单点精细积分[J]. 自然科学进展,1996,6(3):313-316.

[8]富明慧,林敬华. 精细积分法在非线性动力学问题中的应用[J]. 中山大学学报：自然科学版,2008,47(3):1-5.

[9]蓝林华,富明慧,程正阳. 功能梯度材料瞬态热传导问题的降维精细积分法[J]. 固体力学学报,2010,31(4): 406-410.

[10]蓝林华,富明慧,高文乐. 功能梯度材料稳态热传导方程的分层精细指数法[J]. 中山大学学报：自然科学版,2011,50(4): 1-6.

[11]李荣华,冯果沈. 微分方程数值解法[M]. 北京:人民教育出版社,1980.

[12]富明慧,刘祚秋,林敬华. 一种广义精细积分法[J]. 力学学报,2007,39(5):672-677.

[13]富明慧,梁华力. 一种改进的精细-龙格库塔法[J].中山大学学报：自然科学版,2009,48(5):1-5.

[14]戴嘉尊. 微分方程数值解法[M]. 南京:东南大学出版社,2002.