基于直线积边界元法的含热源热传导问题研究

高 宇， 刘 彪， 李通盛， 杨新贵，何钱生， 王 桥， 周 伟， 曹 悦

(1. 大唐宣威水电开发有限公司，云南 宣威 655400；2. 武汉大学 水利水电学院，湖北 武汉 430072)

热传导，是介质内无宏观运动时的传热现象，无论在生活中或工程中，都十分常见，简而言之，只要介质内或者介质之间存在温度差，就一定会发生传热。此外，在工程结构中，由于各部分不能自由伸缩的结构出现温度变化或者内部各部分温度不同而产生的温度应力是结构安全性的关键性影响因素。因此研究结构内部的温度分布是十分必要的。进一步的，混凝土结构在浇筑初期会由于水化热而产生大量热量，或者内部发生火灾释放大量热的建筑结构，可以概括为含热源的工程结构，由于热源的存在，结构内部温度场的分布会发生较大变化，导致温度应力急剧增大或者分布发生较大变化，进而对结构安全性带来巨大的安全隐患。

目前研究热传导问题的数值方法很多，例如无网格法[4]，有限元法[5]，边界点法[6]等。边界元法[7,8]作为一种重要的半解析数值方法，可以对研究问题进行降维，只需要对研究域的边界进行离散，无需在内部划分网格，这是边界元法十分显著的优势。尤其是在研究大尺寸，无限域问题时，其优势更加明显。目前，边界元法在位势问题、声学、波的传播、静力学、动力学等方面都已经得到了广泛应用。对于一般热传导问题，边界元法也能非常精确高效地解决。

而当面临考虑热源的热传导问题时，在边界积分方程中，会出现域积分，如果不加以处理，则需要进行单元离散，使得边界元丧失了降维的优势。因此，寻找一种能够将域积分转化为边界积分的方法是十分必要的。

现存的域积分转化方法中，应用最为广泛的方法是双互易方法(Dual Reciprocity Method，DRM)[11,12]。DRM不需要进行研究域内的单元离散化，而是采用内部节点和一系列指定的径向基函数(Radial Basis Function，RBF)进行逼近。DRM的精度和效率在很大程度上取决于RBF的分布和位置，内部点的位置可以任意选择，但域内用于插值的形函数应具有特殊解，不能任意选择，对于复杂的DRM问题，要获得特殊解不是一件容易的事情。此外，还有用于求解泊松方程的蒙特卡罗积分法(Monte Carlo Integration Method，MCM)[13]。此方法简单易行，但不如DRM精确。近年来，还出现了一种十分强大的方法，称为径向积分法(Radial Integration Method，RIM)[14]，其可以将域积分转化为边界积分和径向积分。一般情况下，径向积分可采用解析或数值计算。二维和三维域积分的计算方法是统一的，也不需要对区域进行离散。

本文采用了直线积分法(Line Integration Method，LIM)[15]处理域积分。直线积分法基于散度定理将域积分转化为包含一维线积分的边界积分，无需对研究域内进行离散。目前已经应用于弹性力学和位势问题中的域积分处理，其有效性和精确度也得到了验证。此外，该方法还十分容易和快速多级算法进行结合以适用于计算大规模工程问题。

1 考虑热源的边界积分方程

结合前人研究，可知控制方程为：

(1)

式中：k为热传导系数；B(x)为热源项；θ为温度；x为研究域Ω内的点，其组成元素为xi。

边界情况可设为：

(2)

通过引入权函数、分部积分法和散度定理，控制方程可写成下面的积分方程：

(3)

(4)

(5)

式中：n为边界Γ的外法向量，其组成元素为ni；r为点x和y的距离，即

r=‖x-y‖

(6)

可以看到，式中含有一项域积分：

(7)

需要利用直线积分法将域积分转化为边界积分。

2 直线积分法

本节将介绍利用直线积分法将二维问题中的域积分转化为边界积分，为了将问题一般化，现假设一个域积分为：

(8)

假设函数f(x,y)可以从矢量场的散度得到，具体表达为：

f(x,y)=divF=∇·F

(9)

式中：div为散度运算符；F为向量场，可表达为：

F=f1(x,y)e1+f2(x,y)e2

(10)

式中：ei为直角坐标系中的一组基向量。

通过运用散度定理，可以将域积分转化为以下形式：

(11)

为了寻找一个合适的函数F满足公式，可以采用特解法，故而式可以改写为：

(12)

或

(13)

式中：nx，ny为向量n在x，y轴方向的分量，并且有

(14)

(15)

通过使用式，可以将域积分式转化为包含直线积分的边界积分，为了简单起见，可以将任意积分线的起始点设为一个常数，x=a，其中a为任意值。并且可以将函数f(x,y)的定义扩展至R2，这样就不需要去对积分点在积分线上的位置做判断。

假设将边界划分为N个边界单元，这样式可改写为：

(16)

(17)

式(17)的计算一般都可以采用数值计算的方式，现在可以将计算域用背景网格覆盖并逐级等分为小网格，进而可以将积分线截断成为一段段的积分线，则式(17)可以改写为：

(18)

式中：Lk为背景网格划分后的直线积分线段。

引入背景网格之后，利用每个网格内的积分线段上的高斯积分点就可以直接进行计算并且可以获得很高的精度。进一步地，域积分可转化为：

(19)

最后，将直线积分法应用于边界积分方程中的域积分的计算，具体转换后的表达式为：

(20)

3 数值验证

为了验证直线积分法在边界元法中域积分转换的有效性和精确度，本文采用了两个例子进行验证，第一个例子为一个承受温度应力的矩形梁，并已知其解析解，进而可以验证本文所采用方法的精度。第二个例子为一个混凝土重力坝，计算结果将与有限元法模拟结果进行对比以验证本方法在大坝温度场分析中的有效性。

3.1 矩形柱的含热源热传导分析

如图1所示，该矩形梁的宽度W为0.2 m，高度H为0.5 m，热传导系数k为0.0025，上下边界温度为0，热源项的具体函数表达式为：

图1 矩形柱的几何模型

HS(y)=2e-y

(21)

相应的温度场的解析解为：

(22)

本算例中将梁的边界离散为200个单元，并将利用直线积分法和边界元法进行计算得到的结果与解析解进行对比。选取了直线x=0.1 m上的9个内部点作为对比，具体结果见表1，结合图2，显而易见，数值解与解析解高度吻合，进而证明了直线积分边界元法的有效性和精确性。

表1 温度值对比

图2 温度的解析解与数值解对比

3.2 坝体结构的含热源热传导分析

大坝模型具体几何尺寸如图3所示，单位为m，热传导系数为1000 W/(m·K)，大坝顶部温度设为38.5 ℃，底部温度为10 ℃，其余边界均为绝热状态。热源假设为常数0.5 W/m3。

图3 重力坝几何模型

为了验证本文方法的正确性，由于缺乏解析解，本文采用有限元法计算的结果与本文方法进行对比。如图4所示，边界元计算结果与有限元计算结果十分相近。

图4 计算结果对比

4 结 论

本文主要基于边界元法，进行了含热源的热传导研究，由于考虑热源导致边界元积分方程中出现了域积分项，为了保证边界元法的降维优点，即只需要对边界进行离散，本文提出了使用直线积分法进行域积分的处理，直线积分法可以基于散度定理，在直角坐标系中，将域积分转化为包含直线积分的边界积分。为了提高直线积分的精度，还可以使用背景网格将直线积分分割为线积分段，之后便可以分别对直线积分段上的积分点进行计算。结合数值验证的第一个例子可以看出，本文所提出方法的计算结果和解析解高度契合，精度很高。此外，还将本文方法运用到大坝结构的热传导分析当中，并与有限元计算结果进行了对比，同样成功地证明了该方法的广泛适用性。