周枫林， 袁小涵， 钦宇， 潘先云， 余江鸿

(湖南工业大学机械工程学院， 株洲 412007)

瞬态热传导问题存在于各种工程领域，例如机械、水利、化工、能源、航天航空等，并随着瞬态热传导问题的广泛应用，传热学也在迅速发展，从刚开始的理论与实验结合的分析方法到各式各样的数值分析方法。在实际问题中，由于结构形状以及边界条件的复杂性，很难获得温度分布的解析解[1]，数值分析方法已成为有效解决问题的方法之一。

根据傅里叶热传导定律，瞬态热传导问题通常由抛物线偏微分方程(PDF)控制。瞬态热传导问题最常用的数值方法是有限差分法(finite difference method，FDM)[2]，因为他在求解瞬态热传导问题中显得非常自然。除了FDM外，还有发展了许多求解瞬态热传导问题的数值计算方法，如有限体积法(finite volume method，FVM)[3]、有限元法(finite element method，FEM)[4]、无网格法(meshless method，MM)[5-9]和边界元法(boundary element method，BEM)[10]，其中，边界元法是一种只涉及边界网格化的数值计算方法，优点是可节省计算时间和存储量，计算精度比较高，边界区域的形状可以任意，求解无限大区域的热传导问题比较理想[11]，但是因为问题的复杂程度，一般找不到问题的基本解，这时边界积分方程中会出现与内部热源等相关域积分。为避免域积分的直接计算，Yang等[12-13]、Feng等[14]在边界元法中引入了径向积分方法，将域积分转化为边界积分。周枫林等[15]将双互易边界元法运用于标量波传播问题。Guo等[16]采用了三重互易法来避免瞬态热传导问题中的域积分。对于双互易法的精度问题，黄诚斌等[17]探讨了域内点的布置位置对计算精度的影响，并提出了可让精度得到进一步提高的一种自适应布点方法。

在边界元法的许多计算运用中，通常采用有限差分格式对时间进行离散，此时计算结果的精度与时间步长的选取有着密切的联系[18]，虽然有限差分边界元法是条件稳定的，但是在大时间步长上，数值结果往往不稳定。为了避免解决上述问题，Yu等[19]将精细时域展开法与双互易边界元法相耦合，避免了直接计算域积分；Yao等[20]将精细积分方法[21]运用于瞬态热传导的边界元模型中；Yao等[22]采用精细积分边界元法解决了双重相变问题；陈豪龙等[23]将双互易边界元法和精细积分方法用于求解二维瞬态热传导问题。

对于双互易精细积分方法在三维热传导问题上还有待深入研究。为此，提出一种双互易法(dual reciprocity method，DRM)与精细积分法(precise integration method，PIM)相耦合的方法解决含内部热源的三维瞬态热传导问题。首先将瞬态热传导问题视为准稳态问题，将温度对时间的导数项视为等效热源。通过DRM将热源的域积分转化为边界积分。对空间变量进行离散后，转化为关于域内温度分布的初值问题，最后运用PIM解决该初值问题，在计算了域内节点的温度后，再通过边界积分方程计算边界的温度和热通量。通过边界元方法对结构进行热传导分析，给工程结构的改进提供一种数据支持，并促进边界元分析方法在热传导问题方面的应用。

1 瞬态热传导的双互易边界元法

1.1 控制方程及其边界条件

具有内部热源的各向同性介质瞬态常系数热传导问题可表示为

(1)

1.2 边界积分方程

对于瞬态热传导问题边界元法，控制方程的边界化是关键步骤，借助加权余量法和三维稳态热传导问题的基本解，将控制方程转换为边界积分方程，可表示为

(2)

对于三维问题有

(3)

1.3 径向基函数逼近

双互易边界元法处理瞬态热传导问题的关键在于对相关项的逼近，用径向基函数(radial basis function,RBF)对已知函数(热源密度函数)和未知函数(与时间相关的温度函数)进行近似表示。取径向基函数fi(x)为

(4)

式(4)中：rdis为RBF插值点到源点的距离；s为形状参数；i为内部节点和边界节点数之和。

(5)

故之前所提到的已知函数和未知函数可表示为

(6)

(7)

式中：Q(x，t)为已知的与位置相关的热源密度函数；αi(t)、βi(t)为未知系数。

系数向量β与α可表示为

(8)

F-1Q(t)=β(t)

(9)

并将式(5)、式(6)和式(7)代入式(2)可得

(10)

令：

(11)

然后与推导式(2)相类似，运用分部积分法将式(10)中的域积分转化为边界积分，可表示为

(12)

1.4 边界离散

式(12)中只涉及边界积分，所以在边界离散后，可得到方程的离散形式为

(13)

式中：j为边界单元个数；m为单元中的节点数，采用八节点二次面单元，m=8；xm为单元中节点的三维空间坐标；Nm为八节点面单元的形函数。

然后在每个边界的插值点上配置场点，将上述方程写为矩阵形式

(14)

δijC(yi)，yi∈Γ

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

将式(8)和式(9)代入式(14)可得

(24)

特别的是，由于Q(x,t)一般为已知的与位置相关的热源密度函数，所以向量Q(t)为已知向量。

1.5 形成微分方程组

式(15)为半离散系统，将其改写为

(25)

假设边界类型不随时间改变，可将边界条件的表达式改写为

c1u(x,t)+c2q(x,t)=c3，x∈Γ

(26)

式(26)中：c3为边界条件；c1和c2为常数值，其值可由边界类型来确定，则有

(27)

所有边界节点上的边界条件可统一表示为

(28)

然后合并式(25)和式(28)可得

(29)

则边界温度和边界热通量可表示为

(30)

令

(31)

(32)

则式(30)简写为

(33)

然后考虑内部点，则有

(34)

(35)

(36)

(37)

(38)

将式(33)代入式(34)得

(39)

令

(40)

(41)

式(39)可简写为

(42)

式(42)为一阶常系数微分方程组。

2 精细积分方法

对于式(42)的一阶常系数微分方程组，可用已知的初始条件求解。

u(x,t0)=u0(x), ∀x∈Ω

(43)

式(42)的通解形式为

(44)

式(44)中：τ为积分时间变量。

(45)

式(45)中：I为单位矩阵。

式(33)称为矩阵指数函数，为使得计算成本最小，采用逐步计算的方法求解式(44)。

(46)

因而，式(46)可化简为

(47)

这是一个代数推进方程，依据式(43)的初始条件，通过该方程计算，可以逐步得到域内温度，在计算完域内温度后，通过式(33)和式(39)可计算出边界上的量。

矩阵指数函数的计算精确度决定了整个系统的一个计算精度，由于矩阵指数函数涉及大量的矩阵乘法，所以计算非常耗时，为了准确计算该矩阵指数函数，采取一种精细计算方法。

式(45)实际上可以看作是泰勒展开式，其收敛速度在很大程度上取决于时间步长。理论上，级数收敛于任意时间步长，然而，与较大的时间步长相比，较小的时间步长下，式(45)可以用很少的截断项来计算。在精细积分方法中，矩阵指数函数中的时间步长被细分为

(48)

式(48)中：M为细分参数。

运用泰勒级数展开，可表示为

=I+Er

(49)

=(I+Er)(I+Er)

=I+2Er+ErEr

(50)

然后重新赋值Er，使得

Er=2Er+ErEr

(51)

可得

=(I+Er)(I+Er)

=I+2Er+Er×Er

(52)

然经过M次重复计算和重新赋值Er后，最终可得

eAΔt=I+Er

(53)

在整个计算过程中，Er的计算涉及了p个矩阵乘法，eAΔt的计算涉及了M个矩阵乘法，特别的是，假如时间步长为常量，那么eAΔt只会被计算一次。并且在时间步推进中，只涉及矩阵的乘法运算。

3 数值算例

为验证所提出方法的有效性、准确性和稳定性，给出3个数值算例，运用双互易精细积分方法获得的计算结果与解析解进行对比。用于评估的相对误差的表达式为

(54)

在这3个算例中，对于双互易方法，使用的近似函数为前面所提到的式(6)，形状参数s为1，对于精细积分方法中的细分参数M和截断参数p，3个案例都取10和6；材料的热导率、热容和密度分别给定为1 W/(m·C)、1 /(kg·C)和1 kg/m-3。

3.1 狄利克雷边界条件下的环面结构

在这个数值算例中，分析图1所示的环形结构上的瞬态热传导问题，圆环的外径和内径分别为26 m和14 m。

狄利克雷边界条件为

(55)

热源项为

Q(x,t)=2 W/m-3

(56)

初始条件为

(57)

解析解为

(58)

在整个计算过程中，如图1所示，总共划分了750个矩形单元，2 381个边界节点，1 473个径向基函数插值点。

考虑0～50 s的时间段，并且选取了6个时间步长用于计算，即0.1、0.2、0.5、1.0、2.5、5 s，图2展示了在不同的时间步长下，内部点计算温度与解析解的相对误差，可以看出，随着时间的变化，相对误差趋近于0，逐渐收敛，并且在时间步比较小的情况下，即使是接近初始时刻，计算结果也有很高的精确度。在本案例中，最佳时间步长为最小的时间步0.1 s，并且在同一时刻，相对误差大小根据时间步长大小排序，越小的时间步，相对误差越小。

同时，对该模型进行有限元分析，运用软件为workbench，计算模型如图3所示，设置单元大小为0.5 m，节点数为169 388，单元数为119 334。

此时考虑一般情况，即解析解不易求得的情况，给定狄利克雷边界条件为

u(x,t)=60 ℃，x∈Γ

(59)

图1 圆环边界元模型Fig.1 Ring boundary element model

图2 温度的相对误差随时间变化情况Fig.2 Variation of relative error of temperature along time

热源项为

Q(x,t)=1 W/m3

(60)

初始条件为

u(x,t0)=u0(x)=0 ℃，x∈Ω

(61)

在两个模型中，取同样的内部点A(7.927，-7.444，0.634 8)，运用相同的时间步长0.1 s，考虑0～10 s的时间段，将边界元温度计算结果与有限元分析温度结果进行对比，如图4所示。

图3 圆环的有限元模型Fig.3 Finite element model of ring

图4 PI-DRBEM与FEM结果比较Fig.4 Comparison of results between PI-DRBEM and FEM

图4表明，双重互易法和精细积分法耦合拥有很高的精度，论证了该方法的有效性和准确性。

3.2 混合条件下的山型结构

此算例考虑在混合边界条件下，一个山型结构的瞬态热传导问题，三维模型及其尺寸如图5所示。边界划分如图6所示。

图5 山型结构几何示意图Fig.5 Geometric diagram of mountain structure

Γ1和Γ2为诺依曼型边界图6 边界划分几何示意图Fig.6 Geometric diagram of boundary division

边界条件为

(62)

q(x,t)=62x2+16t，x1,x2,x3∈Γ1

(63)

其余边界Γ3为狄利克雷型边界，边界条件为

8t(x1+2x2+1.5x3)

(64)

热源项为

Q(x,t)=2x1+4x2+3x3

(65)

初始条件为

u(x,t0)=u0(x)

(66)

解析解为

1.5x3)，x1,x2,x3∈Ω

(67)

在此算例的计算中，总共划分了500个连续的随机四边形单元，1 920个边界节点，164个径向基函数插值点。使用不同的时间步长来验证方法对时间变量的收敛性,相对误差结果如图7所示。

在本算例中，考虑同样的6个时间步长，在混合边界条件下，时间步长趋近于0时，计算结果从一开始就具有很高的精度。稍大的时间步长中，尽管刚开始误差偏大，但随着时间的推移，相对误差快速减小。可以看出，双互易精细积分法的精确性和有效性。

3.3 狄利克雷边界条件下的弯管结构

此算例考虑空心弯管结构，模型结构及尺寸如图8所示。所有表面为狄利克雷边界，其条件为

图7 温度的相对误差随时间变化情况Fig.7 Variation of relative error of temperature along time

R为半径图8 弯管几何结构示意图Fig.8 Geometric diagram of elbow

(68)

热源项、初始条件、解析解与算例2保持一致。在整个计算过程中，总共划分了448个矩形单元，1 718个边界节点，326个径向基函数插值点，插值点分布如图9所示。

图9 弯管RBF插值点分布Fig.9 Distribution of RBF interpolation points of elbow

在本算例中，同样考虑0～50 s时间段，结果如图10所示，可以看出，相对误差在任何一个时间步中都是随着时间推移而减小，并且时间步长取0.1 s时，内部点热通量计算结果的相对误差最小，并且在初始时刻附近也具有较小的相对误差。有趣的是，内部点热通量的计算精度，并不是随着时间步长减小而变高，在选取的6个时间步中，在取中间大小的时间步1.0 s时，相对误差在任何时刻都是最大的。

图10 热通量随时间的相对误差变化Fig.10 Relative error variation of heat flux along time

4 结论

提出了求解含有内部热源的三维瞬态热传导问题的双互易精细积分方法。在处理该热学问题中，首先运用双互易法处理域积分，再利用精细积分方法处理所得到的初始条件问题。得出如下结论。

(1)双互易边界元法和精细积分方法的耦合具有较高的精度和稳定性。

(2)对于不同的结构体和边界条件，并不是时间步长越小越好。存在时间步长减小，但是计算结果的相对误差却增大的情况，但通常情况下，时间步长越接近零，越有可能获得高精度结果。

(3)在无解析解的情况下，该方法依旧有很高的精确性。

不足是该方法的计算效率有待进一步提高，特别的是，在第一个时间步内，需要进行大量的矩阵与矩阵的乘法运算及矩阵求逆运算。而在后续时间步中，只需进行已有矩阵与向量的乘法运算，容易实现并行计算。