二维功能梯度材料热传导格点型有限体积法研究

龚京风， 徐宗著， 宣领宽， 江致远， 郑文利

(武汉科技大学 汽车与交通工程学院， 湖北 武汉 430065)

功能梯度材料(functionally graded materials，FGM)被广泛应用于航空航天、核反应堆、内燃机等领域。由于FGM常被应用于高温高压的恶劣工作环境，FGM结构内的温度场、应力场变化剧烈，其热力性能得到了广泛关注。刘文光等[1]研究了热流密度、陶瓷体积分数指数和热环境等参数对功能梯度圆柱壳热传导的影响；Fu等[2]分析了在弹性基础上，考虑非线性热环境的多孔功能梯度材料圆柱壳的热声响应；许新等[3]基于Eluer-Bernoulli梁理论和单向耦合的热传导理论，研究了FGM微梁的热弹性阻尼；梅靖等[4]研究了考虑厚度方向稳态温度分布的石墨烯增强功能梯度复合材料板条热弹性问题。

Steinberg等[5]提出为承受复杂的温度场，要求材料参数性能在2个甚至3个方向变化，因此，研究二维FGM热传导问题是有必要的。Rahul等[6]在一阶剪切变形理论基础上，对二维温度变化作用下的双向功能梯度圆板进行了自由振动分析。Tang等[7]研究了具有横向和纵向位移耦合的双向功能梯度梁的非线性湿热动力学。这些研究都考虑多方向功能梯度材料的温度场影响，研究热环境下的二维FGM结构性能，准确预测其温度场是前提。

许杨健等[8]推导了复杂边界条件下的二维FGM板的稳态温度场解析解；蒋科[9]在离散域采用数值法，非离散域采用解析法的混合方法，研究了FGM结构的热传导特性。解析方法虽然能获得精确解，但难以应用于复杂工程问题，而数值模拟方法具有更广泛的适用性。Sladek等[10]采用无网格局部Petrov-Galerkin(meshless local Petrov-Galerkin,MLPG)方法分析了非均匀各向异性FGM三维轴对称结构的瞬态热传导问题；Ching等[11]采用MLPG方法求解了功能梯度梁的瞬态热弹性变形。MLPG方法计算量较大，处理结构复杂网格量大的热机械问题存在一定困难。仝国军等[12]采用有限元方法(finite element methed,FEM)研究了非均匀温度场下变物性二维FGM板的瞬态热传导分布问题。针对FGM热传导问题，Gong等[13-14]提出非结构时域有限体积法(finite volume methed,FVM)，该方法将四边形单元处理为双线性单元，有效避免由于将四边形单元作为线性单元处理的数值振荡问题；随后进一步发展该方法，提出交错格点型有限体积法(cell-vetex finite volume methed,CV-FVM)，用于研究活塞、功能梯度多孔材料等热传导问题[15]。当前，基于FEM的数值计算方法在固体热传导分析领域占据主导地位。考虑到FVM在流体力学里的广泛应用及其在固体力学中的成功应用，为了发展统一的多物理场数值求解方法，拟基于FVM进一步探究其在FGM热传导问题中的应用。

为了研究二维FGM材料热传导问题，本文进一步发展交错CV-FVM，将材料物性作为空间坐标的函数，考虑全域的物性参数变化。基于线性三角形单元和双线性四边形单元，建立二维FGM热传导数值求解方法，并验证其正确性，讨论其适用性；同时，在处理复杂FGM热传导问题，考虑大温度梯度区域难以预知情况时，探究CV-FVM在不加密该区域网格情况下的数值稳定性。本文通过研究求解二维FGM问题的CV-FVM解法，为开发拥有自主知识产权的CAE分析软件奠定基础。

1 热传导数值方法

1.1 控制方程

在控制体M上建立导热方程，控制体体积为V，边界面积为S。根据能量守恒定律建立热传导方程：

(1)

式中：ρ、c、T分别是密度、比热容和温度；q为热流密度；n为控制体边界面的外法线矢量；Q为热源在单位体积内产生的热量。

由傅里叶定律可知-q=k▽T，则可将式(1)写成：

(2)

式中k为导热系数，考虑各项同性问题，k为标量。

1.2 数值离散

1.2.1 控制体的建立

对于二维问题，CV-FVM依次将围绕节点的相邻单元中心、边中点连接，构造控制体。采用三角形单元或四边形单元离散计算域。图1为局部网格示意图，空心点表示单元中心或边中点，实心点为单元节点。

图1 二维CV-FVM离散单元示意Fig.1 Schematic diagram of two-dimensional CV-FVM discrete unit

以节点1为例，构造的控制体为图1中虚线围成的区域O1-L2-O2-L3-O3-L4-O4-L1-O1，该区域即围绕节点1的控制体M，其中，L1-O1-L2表示节点1在相邻单元1内的控制体边界面，节点2～7皆为节点1的相邻节点。

采用交错网格的思想，将待解温度定义在网格节点，并假设温度在控制体内均匀分布；将材料参数、体积热源定义在单元中心，假设其在网格单元内均匀分布。由于材料物性的空间分布，在控制体内，物性参数非均匀，从而将物性参数的空间分布引入离散方程。

1.2.2 导热方程的离散

为了方便方程的离散，将式(2)改写为：

(3)

式中nx、ny分别表示控制体边界面的法向矢量在x、y方向的分量。

方程(3)右端第1项借用FEM的思想，采用形函数对该项进行离散：

(4)

由于体积源项定义在单元中心并假设在单元内均匀分布，因此方程(3)右端第2项可离散为：

(5)

方程(3)左端为一阶时间项，采用向后差分的方式离散：

(6)

式中：t为当前时刻;Δt为时间步长。

对于边界上的节点，相应的导热方程离散需要考虑恒温边界SD，热流密度边界SN和换热边界SR的影响：

SD:T=TB

(7)

(8)

(9)

式中：TB为边界温度;qB为热流密度;hB为对流换热系数;Tf为环境温度。

对于边界SD上的节点，采用式直接给定温度。对于边界SN和边界SR上的节点，将式(8)和式(9)代入式(4)得：

(10)

将方程的最终离散形式整理为:

(11)

式中：ABi表示与节点相邻的第i个边界网格面的面积矢量；nN表示与当前节点相邻的位于SN上的边界网格面的个数;nR表示与当前节点相邻的位于SR上的边界网格面的个数。

本文只考虑各项同性且无内热源的稳态问题，因此，式(11)可以简化为：

(12)

1.3 程序设计

为实现上述数值算法，采用C++语言编写程序，图2为程序流程图。

图2 程序流程Fig.2 Program flow chart

为了考虑FGM热传导问题，在物性参数初始化时，根据其变化规律，基于单元中心坐标计算网格单元内的物性参数。为了考虑非均匀温度边界条件，在初始化时，根据边界面单元中心坐标计算非均匀边界参数。同时，计算程序不采用迭代法求解，而采用直接法求解。

将离散后的控制方程整理为式的形式：

AX=B

(13)

(14)

(15)

(16)

式中：e表示节点总数；A中的元素φlm为：

(17)

式中：下标l和m皆表示节点编号；P表示当前节点l的相邻单元所包含节点编号中，存在相邻节点m的单元数。采用文献[17]中的置大数法处理恒温边界。

2 热传导方法验证

2.1 网格适用性的验证

计算内径0.08 m，外径0.1 m的圆环热传导问题。圆环内外两侧皆为恒温边界，内侧温度T1=0 ℃，外侧温度T2=1 ℃。圆环的导热系数沿径向呈指数函数变化：

(18)

式中：k0=17 W/(m·℃)；λ为分布参数，取λ=50 m-1。

分别采用三角形单元、四边形单元、混合网格单元划分计算域，网格尺寸为0.005 m。计算得到的圆环径向温度值(见表1)与文献[18]解析解吻合良好。

表1 不同网格单元在径向上的温度值Table 1 The temperature values of different grid cells in the radial direction

图3为基于不同的网格计算得到的圆环温度分布云图，其径向方向的温度变化皆与解析解相匹配，说明CV-FVM对于线性三角形单元、双线性四边形单元和二者的混合单元皆具有良好的适用性，不存在数值振荡问题。

图3 不同网格单元温度分布云图Fig.3 Cloud map of temperature distribution in different grid cells

2.2 对二维FGM热传导问题的适用性

计算图4所示的FGM板。正方形板边长a=0.01 m。上边界B1给定对流换热边界条件，下边界B2给定热流密度边界条件，左右边界B3、B4给定恒温边界：

图4 功能梯度平板Fig.4 Functional gradient plate

(19)

计算域采用1×10-3m的四边形网格单元进行划分。

假设导热系数沿x和y方向按指数函数规律双向变化，对流换热系数沿x方向按指数规律分布，即：

k(x,y)=k0ecx+dy,h(x)=h0ecx

(20)

式中：k0为x=y=0 m处的导热系数;h0为x=0 m处的换热系数;c、d为导热系数梯度分布参数，当c=d=0 m-1时，表示材料参数均匀分布，退化为普通材料。

考虑对流换热边界温度场的解析解为[19]:

T(x,y)=Tw+

(21)

其中：

Kn=4n2π2+a2c2

Ln=1-cos(nπ)eac/2

分别计算2种条件下二维FGM板的热传导问题。

算例1下边界B2绝热q0=0 W/m2，Tf=100 ℃，Tw=0 ℃，k0=1×103W/(m·℃)，h0=1×106W/(m2·℃)，c=d=200 m-1。将计算得到的x=1×10-3m和x=5×10-3m处的温度与FEM数值解及解析解对比，如图5所示。本文计算得到的温度曲线与FEM解及解析解吻合良好。

图5 第1种情况结果对比Fig.5 Comparison chart of results in the first case

算例2下边界B2给定热流密度q0=1×108W/m2，Tf=800 ℃，Tw=300 ℃，k0=1×103W/(m·℃)，h0=1×106W/(m2·℃)，c=d=100 m-1。分别将基于CV-FVM、FEM计算得到的特征点温度与文献中的解析解对比并计算误差，见表2和表3。可以看出，相同网格尺寸下，本文CV-FVM计算得到的结果与FEM解的求解精度相当。同时，不同非均匀边界条件下，基于CV-FVM计算得到的二维FGM板的热传导结果都与解析解吻合良好，验证了本文CV-FVM对不同边界条件的适用性。

表2 第2种情况CV-FVM结果对比Table 2 Comparison of CV-FVM results in the second case

表3 第2种情况FEM结果对比Table 3 Comparison of FEM results in the second case

为验证CV-FVM与FEM的计算效率，基于Intel Core i7-10750H CPU，内存为16 GB的计算机，将算例2的网格尺寸加密至2×10-5m，分别采用CV-FVM与FEM再次计算，求解耗时及最大内存消耗见表4。可以看出目前在求解时间上CV-FVM比FEM长，但在内存消耗上要比FEM少，这是由于当前版本的软件还未在求解速度上进行优化，在后续研究中会逐步改进。

表4 CV-FVM与FEM计算效率对比Table 4 Comparison of computational efficiency between CV-FVM and FEM

2.3 CV-FVM的优势

基于2.2节的算例1，分别讨论c=d=0，-200 m-1时温度的分布。图6为c=d=0 m-1时的等温度线图。此时导热系数等于k0，材料为均质材料；换热系数等于h0，边界条件均匀。因此，计算得到的温度关于x=5×10-3m对称分布。图7为c=d=-200 m-1时的二维FGM板温度分布，由于导热系数沿x，y方向变化，且边界条件沿x方向变化，等温线不再对称分布，上边界温度梯度较大。考虑到板的两侧为T=0 ℃ 的恒温边界，上边界角点处存在剧烈的温度变化，考察此处数值计算结果，可以体现数值算法求解大梯度问题的性能。

图6 c=d=0时等温线Fig.6 Isotherm diagram when c=d=0

图7 c=d=-200 m-1时等温线Fig.7 Isotherm diagram when c=d=-200 m-1

取y=1×10-2m处的温度分布结果做进一步分析。分别对比网格尺寸为0.25×10-3、1×10-3、2×10-3m时的FEM结果和本文CV-FVM结果，并根据公式算出y=1×10-2m处的温度解析解。可以看到，随着网格尺寸的增加，FEM结果在角点附近出现了不合理的数值波动现象，而本文CV-FVM始终能给出合理的结果。图7中这种等温线汇集在一个点上的现象在物理上不合理，数学上这种位置称为奇点，这是由角点处存在2个温度而造成的[20]。在处理这种情况时，由于CV-FVM中将边界变量存放至边界网格面中心，使用这些边界条件变量时只需按照边界序号依次循环取用，若该处存在2个边界条件，则由循环顺序决定采用哪个边界。受网格精度影响，当网格尺寸较大，此处的温度阶跃变化表现为一个线性变化(见图8)。相比于FEM，CV-FVM处理此类问题时，计算所得的结果整体可信度更高，数值更加稳定。在处理复杂热传导问题时，难以提前预知大温度梯度存在的区域，采用本文CV-FVM，即使不加密该区域的网格，也能得到合理的计算结果。

图8 y=1×10-2 m处的温度曲线Fig.8 Temperature curve when y=1×10-2 m

3 结论

1) 本文发展的用于二维FGM热传导问题的CV-FVM求解方法，能够处理线性三角形单元、双线性四边形单元及混合网格，适用于不同的边界条件类型。

2) 由于采用交错网格技术考虑物性参数的空间变化，没有基于任何分布假设，该方法适用于具有任意材料物性分布的二维FGM热传导问题。

3) 在处理边界角点存在剧烈温差情况时，FEM需要通过加密网格来避免数值振荡，而本文发展的CV-FVM即使在粗糙网格情况下，也能给出合理的数值结果。因此，本文发展的CV-FVM更适用于恶劣环境下、物性参数复杂变化的FGM热传导问题，尤其对于无法提前预知大温度梯度区域的问题，本文方法具有明显优势。