瞬态常系数热传导问题的等几何边界元法研究

Isogeometric Boundary Element Method for Transient Heat Conduction with Constant Coefficients

摘要: 等几何边界元法(IGABEM)在温度场方面的应用研究是边界元法的重要应用之一，相比于传统边界元法，等几何边界元通过样条函数进行几何与物理场的插值近似，可有效降低传统边界元中的离散误差，提高计算精度，因而被广泛应用于数值分析中。采用径向基方法将含有时间梯度域积分项转化成边界积分，并采用NURBS进行几何插值，保持结构外形的精确。通过算例，对比传统BEM与IGABEM的计算结果，验证该算法的正确性与有效性。

Abstract: The application of isogeometric boundary element method (IGABEM) in the temperature field is one of the important applications of the boundary element method. Compared with the traditional boundary element method, the isogeometric boundary element uses the spline function to perform geometric and physical field interpolation approximation. It can effectively reduce the discrete error in the traditional boundary element and improve the calculation accuracy, so it is widely used in numerical analysis. Radial basis method was used to transform the integral term with time gradient domain into boundary integral, and NURBS was used for geometric interpolation to maintain the precise shape of the structure. Through numerical example, the calculation results of traditional BEM and IGABEM are compared to verify the correctness and effectiveness of the algorithm.

1. 引言

边界元法应用于许多方面，主要有热传导、弹性问题、位势问题、声学等。在飞行器的飞行过程中，热防护系统与空气接触的外表面及内表面产生很大的温差，严重影响飞行器的飞行安全，传统边界元法在进行网格离散的过程中会不可避免地产生误差，而且传统边界元计算中会涉及到矩阵的组装与求逆，会大量浪费计算时间，如何降低和解决误差，节省计算时间提升计算效率将是研究人员研究的前沿课题，等几何边界元法(Isogeometric boundary element method, IGABEM)是边界元法与等几何分析相结合，为温度场热传导的发展指明了新的方向，与传统边界元法相比，采用高阶样条插值的等几何边界元法有很多优势，保持几何模型精确同时，高阶近似物理场，提高计算精度，这一优点是传统基于拉格朗日插值算法所不具备的 [1]。等几何边界元的研究起源可以追溯到上世纪60年代，1962年法国雷诺汽车公司的工程师Pierre Bezier提出的Bezier曲线开启了样条函数在曲线曲面建模中的应用，随着工程建模需求的不断提高，非均匀有理B样条(NURBS)在过去数十年间成为了几何造型领域应用最广泛的造型工具 [2]。将等几何边界元与温度场热传导问题相结合，能有效地克服传统边界元的缺点，也能发挥NURBS生成点的优势。近年来高效伟提出的径向积分法将边界元中的域积分转化成边界积分，简化了边界元法的求解过程，提升了边界元法的应用范围 [3] [4]。边界元法的主要步骤可以分为以下几步，第一步是建立边界积分方程，将问题边界化，第二步将边界离散化，采用有限元的离散技巧，由于离散只在边界中进行，误差只产生于边界，而区域内的未知量可以解析的求出，由于基本解存在奇异性使得边界元在解决奇异性问题中精度较高。第三步求解系数矩阵，边界元法计算的优势在于能计算大区域问题，奇异性突出问题。热传导问题广泛存在于航空航天、工业制造、结构材料优化等领域，因而受到许多学者的关注，边界元法已经广泛应用于热传导问题 [5]。因此将等几何方法与边界元方法相结合能更有效服务于温度场热传导问题，最大程度发挥边界元法计算优势 [6]。本文将等几何边界元法用于瞬态温度场分析 [7]，并采用径向基方法将含有时间导数项域积分转化成边界积分，最后通过算例验证该方法的正确性。

2. 控制方程

瞬态问题与时间相关，含有内部热源的各向同性介质热传导控制微分方程表示如下：

$k{\nabla }^{2}T\left(x,t\right)+b\left(x,t\right)=\rho c_x\frac{\partial T\left(x,t\right)}{\partial t}\left(t>t_0,x\in \Omega \right)$ (1)

$T\left(x,t\right)$ 表示t时刻在点x处的温度，k为热导率， $b\left(x,t\right)$ 为热源项， $\rho$ 为材料密度， $c_x$ 为比热容，结构表面接触热源时，一般满足如下三类边界条件：

$\begin{array}{l}T\left(x,t\right)=\bar{T}\left(x,t\right)\left(t>t_0,x\in {\Gamma }_T\right)\\ q\left(x,t\right)=-k\frac{\partial T\left(x,t\right)}{\partial n\left(x\right)}=\bar{q}\left(x,t\right)\left(t>t_0,x\in {\Gamma }_q\right)\\ T_f\left(x,t\right)=\overline{T_f}\left(x,t\right)\left(t>t_0,x\in {\Gamma }_f\right)\end{array}$ (2)

(2)式中， ${\Gamma }_T$ 表示温度恒定边界； ${\Gamma }_q$ 表示温度通量恒定边界； ${\Gamma }_f$ 表示边界接触外部环境温度(对流边界条件)； $n\left(x\right)$ 为边界点处单位外法线方向向量。

瞬态问题与时间相关，需给定初始条件，如下：

$T\left(x,t_0\right)=T_0\left(x\right)\left(x\in \Omega \right)$

$q\left(x,t_0\right)=q_0\left(x\right)\left(x\in \Gamma \right)$

在上式中， $T_0$ 表示域的初始温度； $q_0$ 表示结构边界； $\Gamma$ 是初始温度通量。

3. 瞬态温度场边界元法

通过加权余量操作和分部积分运算，可得到如下的边界–域积分方程，表达如下：

$c\left(x\right)\hat{T}\left(x,t\right)+\left(H^{*}\hat{T}\right)\left(x\right)=\left(G^{*}q\right)\left(x\right)+\left(\tilde{G}b\right)\left(x\right)-\left(\tilde{G}\stackrel{\dot{}}{T}\right)\left(x\right)$ (3)

在式(3)中，当源点x位于结构域内时系数 $c\left(x\right)=1$，位于光滑结构边界上时 $c\left(x\right)=0.5$。另外，

$\{\begin{array}{l}\left(H^{*}\hat{T}\right)\left(x\right):={\displaystyle {\int }_{\Gamma }H\left(x,y\right)\hat{T}\left(y,t\right)\text{d}\Gamma \left(y\right)}\hfill \\ \left(G^{*}q\right)\left(x\right):={\displaystyle {\int }_{\Gamma }G\left(x,y\right)q\left(y,t\right)\text{d}\Gamma \left(y\right)}\hfill \\ \left(\tilde{G}b\right)\left(x\right):={\displaystyle {\int }_{\Omega }G\left(x,\tilde{y}\right)b\left(\tilde{y},t\right)\text{d}\Omega \left(\tilde{y}\right)}\hfill \\ \left(\tilde{G}\stackrel{\dot{}}{T}\right)\left(x\right):={\displaystyle {\int }_{\Omega }G\left(x,\tilde{y}\right)\tilde{\rho }\left(\tilde{y}\right)\stackrel{\dot{}}{\hat{T}}\left(\tilde{y},t\right)\text{d}\Omega \left(\tilde{y}\right)}\hfill \end{array}$ (4)

在式(4)中， $\hat{T}=kt$ 表示标准化温度， $\stackrel{\dot{}}{\hat{T}}=\frac{\partial \hat{T}}{\partial t}$ 表示标准化温度， $G\left(x,y\right)$ 和 $H\left(x,y\right)$ 分别为格林函数及其导数，满足如下表达：

$\{\begin{array}{l}G\left(x,y\right)=\frac{1}{2\pi }\ln \left(\frac{1}{r}\right)\hfill \\ H\left(x,y\right)=\frac{\partial G\left(x,y\right)}{\partial n\left(y\right)}=\frac{-1}{2\pi r}\frac{\partial r}{\partial n\left(y\right)}\hfill \end{array}$ (5)

观察上式(3)可知，最后两项是域积分非边界积分，直接使用边界元法难以进行域积分计算，虽然采用域离散法，比如有限元法、有限差分法等可以有效进行该类积分计算，但失去了边界元法只需在结构表面进行积分计算的优势。本文采用径向基方法将域积分转化成边界积分，保持了边界元法的降维计算优点。

$\left(\tilde{G}b\right)\left(x\right):={\displaystyle {\int }_{\Gamma }\frac{B\left(x,y\right)}{r\left(x,y\right)}}\frac{\partial r}{\partial n\left(y\right)}\text{d}\Gamma \left(y\right)$ (6)

式(6)中， $B\left(x,y\right)$ 为径向基方法处理后得到的基函数，表达如下：

$\begin{array}{l}B\left(x,y\right)={\displaystyle {\int }_0^{r\left(x,y\right)}G}\left(x,\tilde{y}\right)b\left(\tilde{y},t\right)r\left(x,\tilde{y}\right)\text{d}r\left(\tilde{y}\right)\hfill \end{array}$ (7)

上述结果详细推导过程高效伟已经给出，对于简单热源分布函数，基函数B可以被解析求出。对于复杂热源分布函数，可以通过高斯数值积分近似计算得到。然后对由温度对时间导数引起的域积分，即等式右边最后一项进行径向基处理，以转换成边界积分。但由于该项含有未知温度对时间导数项 $\stackrel{\dot{}}{\hat{T}}$，无法直接使用径向积分法将其转换成边界积分，在此我们可以通过插值计算近似得到。

4. 边界积分方程的离散

为了进行数值积分，我们将边界离散成包含Ne个NURBS单元的不重叠集合，表达如下：

$\Gamma ={\displaystyle \underset{e=1}{\overset{N_e}{\cup }},}{\Gamma }_i\cap {\Gamma }_j=0,i\cancel{=}j$ (8)

使用NURBS基函数进行插值近似，在任一NURBS单元e上的变量可近似表达如下：

$\{\begin{array}{l}\hat{T}\left(\xi \right)={\displaystyle \underset{l=1}{\overset{m}{\sum }}{R}_{l,p}}\left(\xi \right){\hat{T}}^l\hfill \\ q^{\left(\xi \right)}={\displaystyle \underset{l=1}{\overset{m}{\sum }}{R}_{l,p}}\left(\xi \right)q^l\hfill \end{array}$ (9)

式(9)中l代表控制节点编号，m为控制点数， $\xi \in \left[-1,1\right]$ 表示局部参数坐标， ${\hat{T}}^l$ 和 $q^l$ 分别为控制节点处规则化温度及其梯度值。将上式带入到式，可重新得到前两项表达，如下：

$\{\begin{array}{l}\left(H^{*}\hat{T}\right)\left(x\right):={\displaystyle \underset{e=1}{\overset{N_e}{\sum }}{\displaystyle \underset{l=1}{\overset{m}{\sum }}\left[{\displaystyle {\int }_{-1}^{1}H}\left(x,y\left(\xi \right)\right)R_{l,p}\left(\xi \right)J^e\left(\xi \right)\text{d}\xi \right]{\hat{T}}^l}}\hfill \\ \left(G^{*}q\right)\left(x\right):={\displaystyle \underset{e=1}{\overset{N_e}{\sum }}{\displaystyle \underset{l=1}{\overset{m}{\sum }}\left[{\displaystyle {\int }_{-1}^{1}G}\left(x,y\left(\xi \right)\right)R_{l,p}\left(\xi \right)J^e\left(\xi \right)\text{d}\xi \right]q^l}}\hfill \end{array}$ (10)

由于源项为0，故可直接移项得：

$\left(H^{*}\hat{T}\right)\left(x\right)-\left(G^{*}q\right)\left(x\right)\text{+}c\left(x\right)\hat{T}\left(x,t\right)=-\left(\tilde{G}\stackrel{\dot{}}{T}\right)\left(x\right)$ (11)

最后将式(11)合并同类项可写为：

$H^{*}\hat{T}-G^{*}q=-C\stackrel{\dot{}}{\hat{T}}$ (12)

可得到两个矩阵方程相减与时间项的关系，式子左侧为边界积分相减，右侧为域积分，这样就将域积分转化成边界积分，通过一步转化将不可积的积分转化成可积的积分，我们一般采用时间推进法进行求解，是一种插值的时间方法。

5. 时间推进法

下式(13)是一组关于时间的积分方程，我们可以假设时间为 $t_0\to t_m$ 并将其分成m份则第i时的温度表达如下：

$t_i=t_0+\frac{t_m-t_0}{m}i$ (13)

则 $t_i$ 时刻温度值为 ${\hat{T}}^i$， $t_{i+1}$ 时刻温度值为 ${\hat{T}}^{i+1}$。

我们采用时间推进法进行方程组的求解，对 $\stackrel{\dot{}}{\hat{T}}$ 的时间导数向前差分格式为：

$\stackrel{\dot{}}{\hat{T}}=\frac{{\hat{T}}^{i+1}-{\hat{T}}^i}{\Delta t}$ (14)

同时，令：

$\{\begin{array}{l}\hat{T}=\beta {\hat{T}}^{i+1}+\left(1-\beta \right){\hat{T}}^i\hfill \\ q=\beta q^{i+1}+\left(1-\beta \right)q^i\hfill \end{array}$ (15)

其中 $\beta$ 取值为0~1，当 $\beta$ 为0时为向前差分格式，当 $\beta$ 为1时为向后差分格式， $\beta$ 为0.5时为中心差分格式。 $i+1$ 与i表示时间步。

$\{\begin{array}{l}\left(H^{*}\hat{T}\right)\left(x\right):=H^{*}\left[\beta {\hat{T}}^{n+1}+\left(1-\beta \right){\hat{T}}^n\right]\hfill \\ \left(G^{*}q\right)\left(x\right):=G^{*}\left[\beta q^{n+1}+\left(1-\beta \right)q^n\right]\hfill \end{array}$ (16)

将上面三个式子(13)、(14)、(15)整理可得：

$\left(H^{*}\hat{T}\right)\left(x\right)-\left(G^{*}q\right)\left(x\right)=-c\left(x\right)\frac{{\hat{T}}^{n+1}-{\hat{T}}^n}{\Delta t}$ (17)

将式子(17)展开后合并同类项，令：

$\{\begin{array}{l}P=\gamma H+c\left(x\right)/\Delta t\hfill \\ Q=\left(1-\gamma \right)H-c\left(x\right)/\Delta t\hfill \end{array}$ (18)

计算第 $n+1$ 个时间步时，第n时刻的温度与通量是已知的，通过前一步迭代出后一步的结果即为时间推进法，将 $n+1$ 时间步的边界条件应用到上式，最终可得矩阵方程：

$B{x}^{n+1}=y^{n+1}$ (19)

式中 $x^{n+1}$ 为第 $n+1$ 个时间步边界热通量及温度组成的 $N_A$ 阶列向量，每一次的时间推进法都会形成线性方程组，对其进行求解可得到时间步的未知量，当边界条件为热通量q已知时，求解后得到的未知量为kt，需要除以热导率k才能得到实际温度值。

6. 椭圆算例

考虑一个长半径为2，短半径为1的椭圆板，初始温度为100℃，板边界接触热源为400℃。材料密度 $\rho =271\text{kg}/{\text{m}}^{\text{3}}$，比热容 $c_x=871\text{J}/\left(\text{kg}\cdot \text{K}\right)$，热导率 $k=202.4\text{W}/\left(\text{m}\cdot \text{K}\right)$，边界对流换热系数 $h=80\text{W}/\left({\text{m}}^2\cdot \text{K}\right)$。我们选取X轴上9个点和Y轴正半轴上9个点进行边界元温度计算。X轴9个点的坐标分别为(−1.6, 0)、(−1.2, 0)、(−0.8, 0)、(−0.4, 0)、(0, 0)、(0.4, 0)、(0.8, 0)、(1.2, 0)、(1.6, 0)，Y轴9个点的坐标分别为(0, 0.9)、(0, 0.8)、(0, 0.7)、(0, 0.6)、(0, 0.5)、(0, 0.4)、(0, 0.3)、(0, 0.2)、(0, 0.1)，由于椭圆模型为轴对称模型，所以下部受热情况类似上部，故下部未选点，两坐标轴均有选点，点的选取具有一般性。我们可以通过细分程序将椭圆周围的点细化成不同的份数，只需要确定初始点数然后进行细分带入到边界元程序中进行计算，通过给出的控制点，不断地插入节点进而拟合给出的椭圆曲线，这种生成点的方法为细分曲线，对于不规则的模型来说，不像椭圆一样有完整的轨迹方程，我们只能通过确定控制点，进行点的细分从而生成所需要的边界节点，即为等几何边界元(IGABEM)相应的NURBS是一种生成节点的方法，这是通过Bezier操作不改变原有控制点的位置，在部分单元上中间点处插入新的控制点，从而逼近模型情况。图1为生成的椭圆模型图。

Figure1. Ellipse model

图1. 椭圆模型图

表1为200S时间下不同步长对不同点处温度值的影响图，其中左边的NSTEP代表的是不同的计算步长，NODE1-NODE9为X轴上的不同坐标点，NTYP (径向基函数的类型) = 3，横向来看越靠近原点(0, 0)温度值越低，其主要原因是受热不均匀，越靠近内部点接受热源越少传热速度越慢，所以温度值变化速率越慢，而关于原点对称的两个点其温度值的相对误差接近10⁻⁷，焦点在X轴上的椭圆，越靠近长轴两端升温速率越快。虽然选取不同时间差分步长会得到不同计算结果，但总体影响不大。一般来说，较小的步长会带来更可靠的结果，对二维问题步长取5是合适的。

Table 1. The effect of asynchronous length on temperature at different points in 200 s

表1. 200 s下不同步长对不同点处温度值的影响

表2为不同的径向基函数类型表，径向基函数的不同类型影响着差值的类型，最常用的是NTYP = 3的逆复合二次插值径向基函数和NTYP = 8的高斯插值函数，未作特殊说明均采用NTYP = 3的插值函数，常用的差值还有薄板样条插值，复合二次插值，四阶样条函数插值等，选取不同的径向基函数类型对温度值也存在一定的影响，但径向基函数只是一种插值函数，它影响的主要是被替代函数温度倒数，其逼近程度决定了计算结果的精确与否，此外多项式的阶数(NPOLY)也是一个影响因素，多项式阶数的大小与计算时间成反比，与计算精度成正比。

Table 2. Different types of radial basis functions

表2. 不同的径向基函数类型表

表3是步长为5，计算时长为200S，不同的径向基函数NTYP对应的不同的温度值，由表可知越靠近内部受热越慢温度值越小，原因类似于步长图内部温度小于外部温度的原因，其变化趋势如图2所示由边界到内部下降速率增加，因为受热程度相似导致对称点上温度值的误差非常小，其温度值高于100℃且小于400℃，由此可见采用不同的径向基函数对温度的影响并不大，径向基函数的变化影响的是插值函数的变化，即为与原函数的误差变化，虽然无法精确的表示出原函数但采用径向基函数可以逼近原函数，所以此解为数值解非解析解，温度值并没有出现大的波动或者跳跃，说明算法的可靠性与稳定性得到验证。

Table 3. The effect of asynchronous length on temperature at different points in 200 s

表3. 200 s下不同的径向基函数对温度的影响图

Figure 2. Comparison of temperature values obtained by traditional BEM and IGABEM

图2. 传统BEM与IGABEM获得的温度值的对比

图2中X轴和Y轴正半轴上点的温度值传统BEM计算结果与IGABEM计算结果对比图，其中左图为不同时刻沿X轴温度分布图，由图可见IGABEM与传统BEM的方法相比几乎没有误差，或者在二维问题中误差极小，由此可验证本文算法的正确性，右图为沿Y轴正半轴温度分布图，越靠近边界温度越高。通过二图不难看出IGABEM的方法在与传统边界元的计算方法相比较情况下还是有一定优势，减少了矩阵的组装求逆过程，节省了计算时间，提高了计算精度，此外还有FLUENT方法也可验证IGABEM算法的正确性。

7. 结论

随着时间趋于无穷大时，整个椭圆板会一直受热，椭圆板内部任意一点都将充分受热，并且温度值都将会趋近于接触热源的温度值。径向基函数的类型是影响温度值的因素之一，但由于径向基函数的选取不同，温度值的相对误差较小，反映了算法的正确性和稳定性，步长对温度值的影响较为明显，对于二维问题来说，步长取5~10是合理的，但步长问题对温度值的影响不可忽略，合理的步长选择是计算温度值正确稳定的前提。相邻两点越靠近边界温度值变化速率越大，因为靠近边界温度传入的速度快，升温速率明显快于靠近内部的点，这对结构的优化设计，飞行器外边界与空气的摩擦降温设计都有重要的意义。

NOTES

^*通讯作者。