基于边界元法的瞬态热传导仿真分析

张驰，张硕

（沈阳航空航天大学民用航空学院，辽宁沈阳110136）

基于边界元法的瞬态热传导仿真分析

张驰，张硕

（沈阳航空航天大学民用航空学院，辽宁沈阳110136）

首先推导了瞬态热传导的边界积分方程，然后通过一系列变换得到了易求解的矩阵形式，提出用迭代法求解瞬态热传导问题.最后引入数值算例，计算了温度分布及热流密度分布，并与解析解进行比较.结果表明采用边界元法所得的数值仿真解与解析解吻合，证明此方法的有效性.

瞬态热传导；边界元法；热流密度；温度场

边界元法相对有限元法能够大大减少计算时间和存储容量，而且边界元不需要对区域内进行离散，因此温度的测量点位置可以任意选择，尤其是这种方法可以直接给出未知表面的温度和测量困难的热流[1-2]，因此边界元法在传热学领域的工程实际应用具有重大的意义.

1 瞬态热传导的基本方程

瞬态热传导问题与稳态相比要涉及与时间相关的问题，因此相对复杂.对于无内热源点的瞬态热传导问题，有

方程的基本解[3]为

其中d是空间维数.基本解（2）代入到基本方程（1）的积分方程中可得

这就是瞬态热传导问题的边界积分方程.最后一项对应着t=0的初始条件.

2 边界元法

对时间域进行划分，假设函数T及q随时间变化，由于二者比T*和q*变化慢得多，故近似在小的时间间隔内为常数，所以可以对（3）式分段积分，并对时间内层求积分得

其中式中D为点源i到边界单元线的垂直距离，Ei项是指数积分函数，可以通过级数来计算[4].

通过上述计算后（3）式可以写成

对空间域划分，边界Γ划分成N个单元，域Ω划分成M个单元，写成矩阵形式即

式中G、H是系数矩阵，只依赖于几何形状、材料特性及时间步长；Pit1表示整个空间域Ω上的初始温度分布对i点的贡献.

在求解边界积分方程时采用迭代解法，把时间域分成N等分，时间步长Δt=tmax/N.这样再把空间域进行划分，根据初始条件就可以计算出边界上的温度场分布和热流密度分布.

3 数值算例

根据以上推导结果，开发C语言程序包，并采用VC++进行严格调试，通过数值算例对上述方法的正确性加以验证.

方板边长L=1，板沿x方向和y方向的热传导率均为k=1，比热c=1，如图1所示.

图1 二维方板的边界条件

其他边界条件为

仿真结果如下所示.

图2 t=0.75时刻y=0边的温度变化

图3 y=1边不同时刻沿x向的温度变化

以上图中边界元解与有限元解及解析解均能较好的吻合，证明该方法的正确性.

4 结论

从上述的仿真分析中，本文所采用的边界元法结果准确.该方法理论推导简单，所得结果误差小且程序编写相对容易.

这种方法为材料的设计及应用其他相关问题如参数优化、灵敏度分析等提供算法依据，在多维坐标及弹性力学、断裂力学、梯度材料等方面的应用也值得探讨.

〔1〕吴洪潭.边界元法在传热学中的应用[M].北京：国防工业出版社,2008.

〔2〕AlokSutradhar,G.H.Paulino.Thesimpleboundaryelementmethodfortransientheatconductioninfunctionally gradedmaterials[J].ComputerMethodsin AppliedMechanicsandEngineering，2004,193(2004)：4511-4539.

〔3〕A.Sutradhar,G.H.Paulion,L.J.Gray.Transientheatconductioninhomogeneousand non-homogeneousmaterialsbytheLaplace transformGalerkinboundaryelementmethod [J].EngineeringAnalysiswithBoundaryElement,2002,26(5):119-132.

〔4〕李慧.三维非稳态热传导边界元方法研究及数值系统开发[D].山东科技大学,2011.

TB131；O414

A

1673-260X（2013）06-0021-02