多域边界面法在稳态热传导问题中的应用*

张见明，李湘贺，陆陈俊，李光耀

(湖南大学 汽车车身先进设计制造国家重点实验室,湖南 长沙 410082)

计算机辅助工程(CAE)对于推动产品研发具有重大的意义.目前CAE的发展已比较成熟，许多成熟的商业CAE软件正在被广泛地采用.但CAE仍面临许多难题[1]，比如如何离散复杂几何的单元才能进行有效计算，如何处理大规模工程问题的数值计算.陆续涌现出的有限差分法(FDM)、有限体积法(FVM)、有限元法(FEM)和边界元法(BEM)[2]都未能很好地解决这些问题.而目前大多数CAE软件采用的是有限元法.

有限元法需要对整个求解域进行离散，将会产生一个很大的代数方程组.对于求解三维(3-D)复杂的实体，尤其是含有细小特征时，离散为可以进行有效计算的实体单元往往比较困难.并且，有限元法的实现基于所求物理问题控制方程和边界条件的等效积分“弱”形式，其试函数要求至少具有一阶连续性.导致应力精度总是比位移精度低一阶，但在实际问题中更关注于应力值，比如产生应力集中的部位及其最大值.

相较于FEM，边界元法(BEM)[3]弥补了有限元法的不足，是一种更加有效的数值方法.它具有等几何分析的特点，使其便于模拟复杂的几何形状.边界元法基于边界积分方程，只需对求解模型的边界进行离散，使求解的问题域降低了一级，大大简化了分析和计算.并且边界积分方程采用问题的解析基本解，具有更高的精度.然而，在传统的边界元法中，三维CAD几何模型被离散成边界元分析模型后，CAD模型的原始几何信息基本被丢掉，这会从根本上导致计算精度的问题，有些甚至对计算起着决定性的作用，并且导致设计和分析成为了两个相互独立的过程.由于CAD模型和分析模型的分离，在边界元自适应网格细分过程中，需要反复地与CAD系统进行交互，而每个交互的过程是很繁琐的.

为了克服上述缺点，张见明教授在边界元法[4]和边界点法的基础上，创造性地提出了边界面法(BFM).在该方法中，对边界的数值积分和场变量的插值都在边界曲面的二维参数空间中进行，CAE分析是直接在CAD模型上进行的，实现了复杂结构的CAE分析自动化.由于CAE模型与CAD几何模型融为了一体，不管网格离散有多么粗糙，分析模型在几何上都是精确的.并且，自适应网格细分过程中,不需要再与CAD系统反复地进行交互,使自适应分析变得简单.

目前边界面法已经取得了许多研究成果，覃先云等人开发了一套参数曲面内网格自动生成的算法[5]，使BFM向实际工程应用迈出一大步.在覃先云的方法中，单元定义和网格的划分都在边界表面的参数空间中进行.谷金良等人提出在BFM中用B样条插值方案来做物理变量的近似[6-8]，并将其成功地应用于稳态热传导问题和弹性静力学问题，并且通过在几何模型中采用相同的插值方案，谷金良等人实现了边界面法的等几何分析.庄超等人实现了基于几何模型的BFM[9]，其几何模型通过细分曲面成型技术构建.谢贵重等人通过引入改进的距离变换、指数变换，提出了一套高效的近奇异积分计算方案[10]，解决了薄型结构分析精度下降的问题.张见明等人通过快速多级子算法、自适应交叉拟合(ACA)和分级矩阵(H-matrix)技术，成功地将计算量级从O(N2)降至O(NlogN)，并且基于GPU高性能并行计算，在CUDA编程环境中实现边界面法正则积分的并行加速[11]，使边界面法的大规模工程应用成为可能.

以上所提及的边界面法的应用，都是基于单域模型.然而，在工程实际分析中，往往存在许多复杂的结构，这些结构需要被分割成多个子域.本文将边界面法的应用延伸至多域模型的稳态热传导问题.在多域问题中，矩阵的组装是至关重要的，文中以3个两两相交的域为例，给出矩阵组装的过程.然后，以某大坝为例，按照大坝的真实施工过程、施工参数和现场实验所得的材料参数，并查阅资料，进行仿真分析，研究坝体的温度分布状况，给出大坝采用多域边界面法进行稳态热传导计算的数值结果，并和有限元法(采用Abaqus软件)的结果进行对比.

1 三维稳态热传导问题的边界面法

1.1 边界积分方程

三维稳态热传导问题可用如下表示：

u,ii=0,∀x∈Ω

(1)

把上述问题转化为等效边界积分形式：

(2)

式中y,s分别为边界上的源点和场点；q=∂u/∂n代表边界法向流量；us(s,y) 和qs(s,y) 代表相应的基本解，它们分别满足：

us,ii=δ(y,s),∀s∈Ω

(3)

qs=∂us/∂n(s)

(4)

对于三维问题，

(5)

其中r是源点和场点之间的距离.

1.2 边界积分方程的离散

本节介绍一种基于表面单元的Lagrange近似.面单元定义在表面的二维参数空间而不是在物理空间或者其他的参数空间，这样所考虑的边界积分的几何数据可以通过参数转换直接进行计算，这样的参数变换与参数空间的映射方案相同.换言之，可以精确得到积分中的几何信息.下面以四节点四边形单元为例(如图1所示)，进行物理变量的近似.

图1 四节点面单元及其坐标变换

在这里构建的形函数如下：

(6)

这些形函数只用于边界表面上的物理变量的近似，几何数据保持精确值，这是和传统边界元法的主要区别.

(7)

其中x=x(u,v),y=y(u,v),z=z(u,v),uk和qk分别是边界节点上的温度值和法向流量，N为所有插值点的个数.采用这样的近似方案后，边界积分方程(2)离散为：

Nk(y))ukdΓ(s)=0

(8)

把场点分布到每一个插值点后，边界积分方程可组装成：

Gq-Hu=0

(9)

其中

(10)

(11)

上式中，Γj为形函数Nk(s)值不为0的边界单元，数量记为Num.

将方程(9)进行变换，使未知量移到左边，已知量移到右边，形成线性组，

AX=b

(12)

式中X是uk或qk的未知数向量.求解方程(12)，就可以得到所有节点k(k=1,2，…，N)上未知量uk或qk的值.根据节点上的值，相应地可以求得边界上和域内任意一点的势和流量值.

2 多域边界面法在三维稳态热传导问题中的应用

本章首先介绍单域问题的矩阵组装，然后仿照单域问题和文献[12]，推导稳态热传导问题的多域边界面法的边界积分方程及其矩阵组装.由于给定第一类(温度、位移)或者第二类(热流密度、面力)边界条件时矩阵的组装较容易，这里我们只介绍含有对流边界的矩阵组装方法.

2.1 单域稳态热问题矩阵的组装

单域问题的边界积分方程为：

(13)

其中Hij,Gij分别为式(9)中的矩阵块，上式可写为：

(14)

令：

(15)

其中，

(16)

βi=-hu，αi=-h

(17)

组装之后为：

(18)

2.2 多域边界面法中矩阵的组装

以3个两两相邻的立方体为例，推导多域问题的边界积分方程，如图2所示．

图2 3个两两相邻立方体示意图

在这个模型中，立方体1和立方体2的边界相交于Γ12,立方体1和立方体3边界相交于Γ13，立方体2和立方体3边界相交于Γ23.两个域相交处的温度边界和热流密度边界分别为：

(19)

(20)

立方体1的边界积分方程如式(21)所示：

(21)

其中，下标d,n,r分别表示第一个域温度边界、热流密度边界、对流边界所对应的系数矩阵块.下标2，3表示第一个域分别与第2、第3个域相交边界所对应的系数矩阵块.加入立方体1的边界条件：

(22)

并且将式(15)代入式(21)，之后把未知量移到左边，已知量移到右边，可得：

(23)

类似立方体1，立方体2和立方体3也分别有如下重新组装的边界积分方程：

(24)

(25)

(26)

求解该方程组得到边界上节点的温度和流量，然后利用边界积分方程(8)，取y为任意域内点，则得到域内任意点的温度.再对边界积分方程(8)两边同时求y点任意方向的导数，即得到y点在任意方向的流量值.由于这种方向求导不是对形函数而是对基本解求的，不会降低拟合精度，因此使得温度和流量的结果具有同阶的计算精度.

3 数值算例

本章分析图3所示水坝结构在固定水温及固定空气温度的边界条件下，坝体内温度的分布情况.

3.1 混凝土水坝模型

如图3所示，考虑对称性，模型厚度取为20 m.基岩高270 m，宽471 m.基岩上的坝体高度为171 m，不均匀地分为62层(与建设过程保持一致)，第一浇筑层底部与基岩接触部分宽157 m.大坝模型左侧为上游面，右侧为下游面，上游面假设蓄水高度为375 m，距坝体最顶端9 m.

图3 某重力坝段结构示意图

3.2 大坝的材料参数

坝段建造采用两种不同的混凝土材料：坝体最底端两层采用常态混凝土，基岩材料视为与第一、二浇筑层相同，坝体其他层采用碾压混凝土.材料参数见表1．

表1 大坝材料参数

3．3 边界条件

基岩左右两侧及底部绝热，模型厚度方向前后面绝热.上游面375 m高度以上、顶部表面、下游面添加空气边界，当地年平均气温为20.7 ℃.上游面蓄水区添加水温边界条件，其中水温边界如下：

1)水温在深度0～123.4 m随深度拟合为线性变化，其拟合函数为：

Tw=20.7-0.059 157 2×h

(27)

2)在距水面123.4 m以下，温度取深水温度13.4 ℃.

3.4 网 格

为做对比分析，除了使用边界面方法分析之外，还采用基于有限元方法的商用软件ABAQUS 11.0对此结构进行稳态热分析，由于本模型中大坝每层的厚度较薄，因此网格较密.在分析过程中一共使用了47 432个二次六面体单元，共计272 300个计算节点，网格模型如图4所示.

图4 有限元网格

在使用边界面法进行分析的过程中，共划分6 494个二次单元(包括三角形单元和四边形单元)，共29 139个计算节点，网格模型如图5所示.

3.5 求解结果及对比

图6和图7分别为ABAQUS与BFM温度分布的计算结果对比.

图5 BFM网格

图6 有限元计算结果

图7 BFM计算结果

对比上面两图，可以看出边界面法与有限元法分析结果基本相同，温度最大值均为20.7 ℃，且温度等值线高度一致.为做更加清晰详细的对比，我们取如图8所示的两个截面上的节点温度，绘成温度曲线进行对比．

图8 截面示意图

以x轴为横坐标，绘出上面两个截面上各个节点的温度值，组成温度曲线，如图9所示．

x

可见，对于有限元法和边界面法的计算结果，截面1和截面2上的温度曲线都保持高度一致.上述两方面的对比可说明边界面法的计算结果与有限元基本一致.

因此，边界面法应用于大规模工程结构的稳态热分析，具有良好的精度.而且，在边界面法的计算中，计算时长共用了83 s，有限元法用了257 s.并且，由于边界面法是直接基于几何模型进行操作，无需定义多个域的装配以及面间的接触等，因此边界面法消耗较少的人力和计算机资源，具有更高的计算效率.

4 结 论

本文用边界面法解决了含有62个浇筑层的大坝的稳态热传导问题.通过和有限元法的对比，证明了计算结果的正确性.说明本文所提出的多域边界面法可以应用在大规模工程问题中的稳态热分析中，并且更加便捷，具有更高的效率.

在后续工作中，我们将实现多域边界面法在瞬态热传导问题中的应用，并且考虑通过优化矩阵的组装方式，来进一步提升计算效率.

[1]ZHANG J M,TANAKA M.Fast HdBNM for large-scale thermal analysis of CNT-reinforced composites[J].Computational Mechanics,2008,41:777-787.

[2]龙述尧.边界单元法概论[M].北京：中国科学文化出版社,2002.

LONG Shu-yao.The introduction of boundary element method[M].Beijing:China Science and Culture Press,2002.(In Chinese)

[3]HANG J M,YAO Z H,LI H.A hybrid boundary node method[J].International Journal for Numerical Methods in Engineering,2002,53:751-763.

[4]ZHANG J M,YAO Z H,TANAKA M.The meshless regular hybrid boundary node method for 2-D linear elasticity[J].Engineering Analysis with Boundary Elements,2003,27:259-268.

[5]QIN X Y,ZHANG J M,LI G Y,etal.A finite element implementation of the boundary face method for potential problems in three dimensions[J].Engineering Analysis with Boundary Elements,2010,34:934-943.

[6]GU J L,ZHANG J M,LI G Y.Isogeometric analysis in BIE for 3-D potential problem[J].Engineering Analysis with Boundary Elements,2012,36:858-865.

[7]GU J L,ZHANG J M,SHENG X M,etal.B-spline approximation in boundary face method for three dimensional linear elasticity[J].Engineering Analysis with Boundary Elements,2011,35:1159-1167.

[8]GU J L,ZHANG J M,SHENG X M.The boundary face method with variable approximation by b-spline basis function[J].International Journal of Computational Methods,2012,9:9-18.

[9]ZHUANG C,ZHANG J M,QIN X Y,etal.Integration of subdivision method into boundary element analysis[J].International Journal of Computational Methods,2012,9:19-29.

[10]XIE G Z,ZHANG J M,QIN X Y,etal.New variable transformations for evaluating nearly singular integrals in 2D boundary element method[J].Engineering Analysis with Boundary Elements,2011,35:811-817.

[11]张见明，余列祥，刘路平.基于GPU加速的边界面法正则积分的研究[J].湖南大学学报:自然科学版,2013,40(3):41-45.

ZHANG Jian-ming,YU Lie-xiang,LIU Lu-ping.Research on regular integration in boundary face method based on GPU-acceleration[J].Journal of Hunan University:Natural Sciences,2013,40(3):41-45.(In Chinese)

[12]肖金生,钱耀南,冉一元.柴油机活塞热负荷的轴对称分区边界元分析[J].武汉水运工程学院学报,1986,10(2):69-78.

XIAO Jin-sheng,QIAN Yao-nan,RAN Yi-yuan.Axisymmetric subregion boundary element analysis of thermal loading of the disel engine’s piston[J].Journal of Wuhan University of Water Transportation Engineering,1986,10(2):69-78.(In Chinese)