边界节点法模拟瞬态热传导及Matlab工具箱开发

边界节点法模拟瞬态热传导及Matlab工具箱开发

万长风;林继;洪永兴

【摘要】针对高温条件下材料热传导问题,采用边界节点法数值模拟二维和三维瞬态热传导问题并编写了相应的计算软件。边界节点法是无网格方法的一种,采用径向基函数的非奇异通解来满足控制方程。对于热传导方程,采用Implicit-Euler差分离散其时间项,将不含时间项的非齐次亥姆赫兹方程以特解与通解的和的形式给出。使用双重互易法在计算区域配点求出方程的特解,使用边界节点法在边界配点来求出方程的齐次解。使用结果表明：该方法计算精度高,数学简单,收敛快,能很好地数值模拟二维及三维瞬态热传导问题。基于上述的算法,使用Matlab语言开发了一个热传导问题的数值模拟工具箱,该工具箱具有简单易操作、代码开源、界面友好等特点。

【期刊名称】《能源与环保》

【年(卷),期】2017(039)002

【总页数】6页(P12-17)

【关键词】瞬态热传导边界节点法双重互易法差分格式 Matlab工具箱

【作者】万长风;林继;洪永兴

【作者单位】河海大学力学与材料学院,江苏南京211100

【正文语种】中文

【中图分类】O241.82

某些工程材料在一些热处理过程中需要在高温条件下长期工作，因此研究高温环境中材料的传热有着重要的意义[1]。由于实验过程耗时费力而且有时对设备要求较高，因此研究一种能够简单效率地计算热传导问题，能为热处理工艺和新兴材料的研发与应用提供一定的参考价值。

国内外很多学者已经使用一些常见算法例如有限元方法来模拟热传导问题[2-7]。

有些网格类算法需要进行复杂的网格划分，在解决一些复杂问题时常会遇到网格划分困难，耗时费力的网格重构等问题，这些问题直接影响着计算的效率和计算精度。为了克服网格类算法的这些不易解决的缺点，无网格方法受到学术界的广泛重视，文献[8]中给出了近年来无网格粒子类方法的研究进展。相比于网格类算法，此类

算法使用边界或区域内配点的信息来求得插值矩阵，不需要在区域内构建网格，因此能够快速收敛且后处理简单。无网格算法已经被广泛地应用于各种科学和工程问题的研究[5-6，9-10]。例如基本解方法是由Kurprade和Aleksidze在1964年

提出[11]，基本解法只需得到配置点的信息并求得插值矩阵，因此比网格类算法的收敛速度快。但是在源点上基本解会出现奇异性问题，为了避免该问题，该方法在物理边界外设置了虚假边界，而虚假边界的设置一般凭借经验，导致数值结果变得不稳定[12]。

为了避免虚假边界的困扰，陈文教授提出了仅在物理边界上离散而不设置虚假边界的边界节点法(Boundary knot method，BKM)[13-14]，该方法使用径向基函数

非奇异通解来满足控制方程，所以该方法不需要在计算区域边界外设置额外的虚假边界且保持了计算的稳定性与收敛性。本文采用边界节点法并结合双重互易法(Dual reciprocity method，DRM)[15-17]模拟瞬态热传导问题。首先，基于Implicit-Euler差分格式离散热传导方程中的时间项，那么原热传导问题的解可写成特解和通解2部分。齐次解可由满足控制方程的非奇异一般解的线性组合来近似，特解可以用双重互易法选取适当的径向基函数并在计算区域配点求解方程来得

到，源项由径向基函数近似从而可以有效的给出热传导问题的数值解。最后基于上述算法，本文介绍了以上述算法为原理编写的Matlab工具箱，此工具箱可解决瞬态热传导问题，为数值方法的比较研究、新的算法讨论提供参考数值结果；也可用于工程问题的计算，为各种工程问题提供较为准确的参考。

假定瞬态热传导方程的计算区域为Ω，边界为Γ=∂Ω，热传导方程可表示为式(1)：边界条件约束如下：

初始条件如下：

T(x，0)=T0，x∈Ω

其中，T(x，t)为在t时刻x点的温度；k为材料导热系数；ρ为材料的密度；c为

比热容；Q(x，t)为材料内部热源；为边界上的温度；为边界上的法向热流量；T0

为任意区域内一点的初始时刻的温度；Γ=ΓD+ΓN。

使用Implicit-Euler差分格式来离散方程中的时间项，在任意一段时间间隔[tn，

tn+1]⊂[0，t]内，物体任一点的温度T(x，t)由T(x，tn+1)来近似，任一点的温度对时间的一阶导数由整个时间段内温度的平均变化率来近似，任一点的内热源Q(x，t)由Q(x，tn+1)来近似[18]：

T(x，t)=θT(x，tn+1)+(1-θ)T(x，tn)

Q(x，t)=θQ(x，tn+1)+(1-θ)Q(x，tn)

其中，n为步数；dt=tn+1-tn为步长；θ为可指定的参数，这里取θ=1，即Implicit-Euler方式。

将方程代入原始方程中：

方程左边的温度t=tn+1的求解需要前一个时间步t=tn的近似解Tn(x)，所以使

用方程的初始条件进行迭代。方程也可简化为：

(▽2-λ2)T(x)=f(x)，x∈Ω

其中

方程的解可写成以下特解与齐次解和的形式：

T(x)=Tp(x)+Th(x)

其中，Tp(x)为方程的特解；Th(x)为齐次解。特解Tp(x)满足式(13)：

(▽2-λ2)Tp(x)=f(x)，x∈Ω

需要注意的是，Tp(x)只需要满足方程而不需要满足边界条件，而齐次解Th(x)满足边界条件而不需要满足方程：

(▽2-λ2)Th(x)=0，x∈Ω

2.1 双重互易法

为了求解得到特解Tp(x)，需引入双重互易法用如下的线性组合来近似方程的右端项f(x)：

其中，αi为待求的插值矩阵；NI为配置点个数；φi(x)为选取的径向基函数。通过在区域内配点可以求解出未知系数αi，从而方程的特解Tp(x)可以表示为式(18)：其中，ψi(x)满足：

(▽2-λ2)ψi(x)=φi(x)

本文选取如下的径向基函数φ(x)：

φ=r2lnr

根据方程，相应的ψ(x)的取值如下：

其中，γ≈0.577 215 664 901 532 8，为欧拉常数；k0为零阶第二类修正的贝塞尔函数。

2.2 边界节点法

本节求解齐次解Th(x)，采用边界节点法引入控制方程的非奇异一般解的线性组合来近似[14]

式中，{βi}为待求系数；为修正的Helmholtz算子(Δ-λ2)的非奇异通解；{xi}为源点集合；Ns为源点总数。