二维功能梯度材料热传导格点型有限体积法研究_龚京风.pdf

1、第 44 卷第 2 期2023 年 2 月哈 尔 滨 工 程 大 学 学 报Journal of Harbin Engineering UniversityVol.44.2Feb.2023二维功能梯度材料热传导格点型有限体积法研究龚京风,徐宗著,宣领宽,江致远,郑文利(武汉科技大学 汽车与交通工程学院,湖北 武汉 430065)摘 要:为了分析二维功能梯度材料热传导问题,发展格点型有限体积法,本文利用交错网格技术处理任意空间分布的材料物性,并使用 C+语言开发数值求解程序。利用该程序分析不同网格类型下功能梯度材料圆环和不同边界条件下功能梯度材料板的热传导问题。结果表明:发展的格点型有限体积法对

2、于线性三角形单元及双线性四边形单元的适用性良好,不同边界条件下格点型有限体积法解与解析解的误差不超过0.12%,验证了求解程序的正确性。与有限元及解析解结果对比发现,当结构角点处存在剧烈温度变化且网格较为粗糙时,有限元结果出现数值振荡的现象,而格点型有限体积法结果仍然合理。分析表明:对于复杂的功能梯度材料热传导问题,格点型有限体积法比有限元方法具有更好的数值稳定性。关键词:交错网格;格点型有限体积法;二维功能梯度材料;热传导;数值振荡;C+程序语言;误差分析;CAE 软件DOI:10.11990/jheu.202105015网络出版地址:https:/ 文献标志码:A 文章编号:1006-70

3、43(2023)02-0321-08Investigation on cell-vertex finite volume method for the heat conduction of two-dimensional functionally graded materialsGONG Jingfeng,XU Zongzhu,XUAN Lingkuan,JIANG Zhiyuan,ZHENG Wenli(School of Automobile and Traffic Engineering,Wuhan University of Science and Technology,Wuhan 4

4、30065,China)Abstract:To analyze the heat conduction problem of two-dimensional functionally graded materials(FGMs),the staggered cell-vertex finite volume method(CV-FVM)is developed in this research.Staggered grid technology is ap-plied to process the material properties of arbitrary spatial distrib

5、ution,and the C+language is used to develop the numerical solution program.This program is used to analyze the heat conduction of FGM rings under different grid types and FGM plates under different boundary conditions.The results show that the developed CV-FVM method has good applicability for linea

6、r triangular elements and bilinear quadrilateral elements.The error between the CV-FVM and analytical solutions under different boundary conditions does not exceed 0.12%,which verifies the correctness of the solution procedure.Comparing the finite element and analytical results shows that when sever

7、e temperature changes occur at the corners of a structure with a coarse grid,the finite element results exhibit numerical oscillations,but the CV-FVM results are still reasonable.The analyses indicate that compared with the complicated finite element method,CV-FVM has better numerical stability for

8、complicated heat conduction analysis of FGMs.Keywords:staggered grid;cell-vertex finite volume method;functionally graded materials;heat conduction;nu-merical oscillation;C+programming language;deviation analysis;CAE software收稿日期:2021-05-01.网络出版日期:2022-10-28.基金项目:国家自然科学基金青年基金项目(51909197);湖北省自然科学基金项目

9、(2019CFB136);武汉科技大学预研基金项目(GF201911).作者简介:龚京风,女,副教授,硕士生导师;宣领宽,男,高级工程师.通信作者:宣领宽,E-mail:nsmpf_wust .功能梯度材料(functionally graded materials,FGM)被广泛应用于航空航天、核反应堆、内燃机等领域。由于 FGM 常被应用于高温高压的恶劣工作环境,FGM 结构内的温度场、应力场变化剧烈,其热力性能得到了广泛关注。刘文光等1研究了热流密度、陶瓷体积分数指数和热环境等参数对功能梯度圆柱壳热传导的影响;Fu 等2分析了在弹性基础上,考虑非线性热环境的多孔功能梯度材料圆柱壳的热声响

10、应;许新等3基于 Eluer-Bernoulli 梁理论和单向耦合的热传导理论,研究了 FGM 微梁的热弹性阻尼;梅靖等4研究了考虑厚度方向稳态温度分布的石墨烯增强功能梯度复合材料板条热弹性问题。Steinberg 等5提出为承受复杂的温度场,要求材料参数性能在 2 个甚至 3 个方向变化,因此,研究哈 尔 滨 工 程 大 学 学 报第 44 卷二维 FGM 热传导问题是有必要的。Rahul 等6在一阶剪切变形理论基础上,对二维温度变化作用下的双向功能梯度圆板进行了自由振动分析。Tang等7研究了具有横向和纵向位移耦合的双向功能梯度梁的非线性湿热动力学。这些研究都考虑多方向功能梯度材料的温度场

11、影响,研究热环境下的二维 FGM 结构性能,准确预测其温度场是前提。许杨健等8推导了复杂边界条件下的二维FGM 板的稳态温度场解析解;蒋科9在离散域采用数值法,非离散域采用解析法的混合方法,研究了FGM 结构的热传导特性。解析方法虽然能获得精确解,但难以应用于复杂工程问题,而数值模拟方法具有更广泛的适用性。Sladek 等10采用无网格局部 Petrov-Galerkin(meshless local Petrov-Galerkin,MLPG)方法分析了非均匀各向异性 FGM 三维轴对称结构的瞬态热传导问题;Ching 等11采用 MLPG方法求解了功能梯度梁的瞬态热弹性变形。MLPG方法计算

12、量较大,处理结构复杂网格量大的热机械问题存在一定困难。仝国军等12采用有限元方法(finite element methed,FEM)研究了非均匀温度场下变物性二维 FGM 板的瞬态热传导分布问题。针对 FGM 热传导问题,Gong 等13-14提出非结构时域有限体积法(finite volume methed,FVM),该方法将四边形单元处理为双线性单元,有效避免由于将四边形单元作为线性单元处理的数值振荡问题;随后进一步发展该方法,提出交错格点型有限体积法(cell-vetex finite volume methed,CV-FVM),用于研究活塞、功能梯度多孔材料等热传导问题15。当前,基

13、于 FEM 的数值计算方法在固体热传导分析领域占据主导地位。考虑到 FVM 在流体力学里的广泛应用及其在固体力学中的成功应用,为了发展统一的多物理场数值求解方法,拟基于 FVM 进一步探究其在 FGM 热传导问题中的应用。为了研究二维 FGM 材料热传导问题,本文进一步发展交错 CV-FVM,将材料物性作为空间坐标的函数,考虑全域的物性参数变化。基于线性三角形单元和双线性四边形单元,建立二维 FGM 热传导数值求解方法,并验证其正确性,讨论其适用性;同时,在处理复杂 FGM 热传导问题,考虑大温度梯度区域难以预知情况时,探究 CV-FVM 在不加密该区域网格情况下的数值稳定性。本文通过研究求解

14、二维FGM 问题的 CV-FVM 解法,为开发拥有自主知识产权的 CAE 分析软件奠定基础。1 热传导数值方法1.1 控制方程 在控制体 M 上建立导热方程,控制体体积为 V,边界面积为 S。根据能量守恒定律建立热传导方程:V(cT)tdV=S-qndS+VQdV(1)式中:、c、T 分别是密度、比热容和温度;q 为热流密度;n 为控制体边界面的外法线矢量;Q 为热源在单位体积内产生的热量。由傅里叶定律可知-q=kT,则可将式(1)写成:V(cT)tdV=SkTndS+VQdV(2)式中 k 为导热系数,考虑各项同性问题,k 为标量。1.2 数值离散1.2.1 控制体的建立 对于二维问题,CV

15、-FVM 依次将围绕节点的相邻单元中心、边中点连接,构造控制体。采用三角形单元或四边形单元离散计算域。图 1 为局部网格示意图,空心点表示单元中心或边中点,实心点为单元节点。图 1 二维 CV-FVM 离散单元示意Fig.1Schematic 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 内的控制体边界面,节点 27 皆为节点1 的相邻节点。采用交错

16、网格的思想,将待解温度定义在网格节点,并假设温度在控制体内均匀分布;将材料参数、体积热源定义在单元中心,假设其在网格单元内均匀分布。由于材料物性的空间分布,在控制体内,物性参数非均匀,从而将物性参数的空间分布引入离散方程。1.2.2 导热方程的离散 为了方便方程的离散,将式(2)改写为:VcTtdV=SkTxnx+Tyny()dS+VQdV(3)式中 nx、ny分别表示控制体边界面的法向矢量在 x、223第 2 期龚京风,等:二维功能梯度材料热传导格点型有限体积法研究y 方向的分量。方程(3)右端第 1 项借用 FEM 的思想,采用形函数对该项进行离散:SkTxnx+Tyny()dS=ni=1

17、kimj=1TijSiNijxnxdS+mj=1TijSiNijynydS()(4)式中:n 表示当前节点的相邻单元总数;m 表示第 i个相邻单元的节点总数;Si表示当前节点在第 i 个相邻单元内的控制体边界面积;Si(Nij/x)nxdS 以及Si(Nij/y)nydS 分别表示第 i 个单元内第 j 个节点上的形函数对 x、y 的偏导与积分路径的乘积,该项只与网格有关,具体计算方法参考文献16,对于不同类型的网格单元,形函数不同;下标 i 表示第i 个单元中心的变量值,下标 ij 表示在第 i 个单元内第 j 个节点上的变量值。由于体积源项定义在单元中心并假设在单元内均匀分布,因此方程(3

18、)右端第 2 项可离散为:VQdV=ni=1QiVim(5)方程(3)左端为一阶时间项,采用向后差分的方式离散:VcTtdV=ni=1iciVim()Tt-Tt-tt()(6)式中:t 为当前时刻;t 为时间步长。对于边界上的节点,相应的导热方程离散需要考虑恒温边界 SD,热流密度边界 SN和换热边界 SR的影响:SD:T=TB(7)SN:-kTn=qB(8)SR:-kTn=hB(T-Tf)(9)式中:TB为边界温度;qB为热流密度;hB为对流换热系数;Tf为环境温度。对于边界 SD上的节点,采用式直接给定温度。对于边界 SN和边界 SR上的节点,将式(8)和式(9)代入式(4)得:SkTxn

19、x+Tyny()dS=ni=1kimj=1TijSiNijxnxdS+mj=1TijSiNijynydS()-SNqBdS-SRhB(T-Tf)dS(10)将方程的最终离散形式整理为:ni=1iciVim()Tt-Tt-tt()=ni=1kimj=1TtijSiNijxnxdS+mj=1TtijSiNijynydS()-nNi=1qBiABi-nNi=1hBi(Tt-Tfi)ABi-ni=1QiVim(11)式中:ABi表示与节点相邻的第 i 个边界网格面的面积矢量;nN表示与当前节点相邻的位于 SN上的边界网格面的个数;nR表示与当前节点相邻的位于 SR上的边界网格面的个数。本文只考虑各项同

20、性且无内热源的稳态问题,因此,式(11)可以简化为:ni=1kimj=1TijSiNijxnxdS+mj=1TijSiNijynydS()-nNi=1qBiABi-nNi=1hBi(Tt-Tfi)ABi=0(12)1.3 程序设计 为实现上述数值算法,采用 C+语言编写程序,图 2 为程序流程图。图 2 程序流程Fig.2 Program flow chart为了考虑 FGM 热传导问题,在物性参数初始化时,根据其变化规律,基于单元中心坐标计算网格单元内的物性参数。为了考虑非均匀温度边界条件,在初始化时,根据边界面单元中心坐标计算非均匀边界参数。同时,计算程序不采用迭代法求解,而采用直接法求解

21、。323哈 尔 滨 工 程 大 学 学 报第 44 卷将离散后的控制方程整理为式的形式:AX=B(13)式 中:A 表 示 系 数 矩 阵,与S(N/x)nxdS、S(N/y)nydS、k 以及边界条件有关;B 表示已知列向量,与边界条件、时间项、源项有关;X 为待解变量即温度。A=11121e21222ee1e2ee|ee(14)B=b1b2beT(15)X=T1T2TeT(16)式中:e 表示节点总数;A 中的元素 lm为:lm=0,m 不是相邻节点时Pi=1kiSiNimxnxdS+SiNimynydS(),m 为相邻节点时Pi=1kiSiNimxnxdS+SiNimynydS()-nN

22、i=1hBABi,m 为相邻节点且为换热边界节点|(17)式中:下标 l 和 m 皆表示节点编号;P 表示当前节点l 的相邻单元所包含节点编号中,存在相邻节点 m的单元数。采用文献17中的置大数法处理恒温边界。2 热传导方法验证2.1 网格适用性的验证 计算内径 0.08 m,外径 0.1 m 的圆环热传导问题。圆环内外两侧皆为恒温边界,内侧温度 T1=0,外侧温度 T2=1。圆环的导热系数沿径向呈指数函数变化:k=k0ex2+y2-R1()(18)式中:k0=17 W/(m );为分布参数,取=50 m-1。分别采用三角形单元、四边形单元、混合网格单元划分计算域,网格尺寸为 0.005 m。

23、计算得到的圆环径向温度值(见表 1)与文献18解析解吻合良好。表 1 不同网格单元在径向上的温度值Table 1The temperature values of different grid cells in the radial direction径向距离/m解析解/三角形单元/四边形单元/混合单元/0.08000000.083 0.260 0.257 0.257 0.249 0.087 0.476 0.469 0.468 0.466 0.0900.649 0.643 0.643 0.649 0.093 0.788 0.786 0.786 0.784 0.097 0.905 0.904 0

24、.904 0.901 0.1111 1 图 3 为基于不同的网格计算得到的圆环温度分布云图,其径向方向的温度变化皆与解析解相匹配,说明 CV-FVM 对于线性三角形单元、双线性四边形单元和二者的混合单元皆具有良好的适用性,不存在数值振荡问题。2.2 对二维 FGM 热传导问题的适用性 计算图 4 所示的 FGM 板。正方形板边长 a=0.01 m。上边界 B1给定对流换热边界条件,下边界B2给定热流密度边界条件,左右边界 B3、B4给定恒温边界:B1:-k(x,y)T(x,y)y=h(x,y)(T(x,y)-Tf)B2:-k(x,y)T(x,y)y=q0B3:T(x,y)=TwB4:T(x,y

25、)=Tw|(19)计算域采用 110-3 m 的四边形网格单元进行划分。假设导热系数沿 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+n=18h0(Tf-Tw)nKnLne-cx/2sinnxaMnh0e(-d+n)b/2+-d+nd+ne(-d-n)b/2()+Nn|(2

26、1)其中:n=c2+d2+4n22a2Kn=4n22+a2c2Ln=1-cos(n)eac/2423第 2 期龚京风,等:二维功能梯度材料热传导格点型有限体积法研究Mn=e(-d+n)y/2+-d+nd+ne(-d-n)y/2Nn=k0edb-d+n2(e(-d+n)b/2-e(-d-n)b/2)图 3 不同网格单元温度分布云图Fig.3Cloud map of temperature distribution in differ-ent grid cells 分别计算 2 种条件下二维 FGM 板的热传导问题。算例 1 下边界 B2绝热 q0=0 W/m2,Tf=100,Tw=0,k0=11

27、03 W/(m),h0=1106 W/(m2),c=d=200 m-1。将计算得到的 x=110-3 m 和 x=510-3 m 处的温度与 FEM 数值解及解析解对比,如图 5 所示。本文计算得到的温度曲线与 FEM 解及解析解吻合良好。图 4 功能梯度平板Fig.4 Functional gradient plate图 5 第 1 种情况结果对比Fig.5 Comparison chart of results in the first case算例 2 下 边 界 B2给 定 热 流 密 度 q0=1108 W/m2,Tf=800,Tw=300,k0=1 103 W/(m ),h0=1

28、106 W/(m2),c=d=100 m-1。分别将基于 CV-FVM、FEM 计算得到的特征点温度与文献中的解析解对比并计算误差,见表 2 和表 3。可以看出,相同网格尺寸下,本文CV-FVM 计算得到的结果与 FEM 解的求解精度相当。同时,不同非均匀边界条件下,基于 CV-FVM计算得到的二维 FGM 板的热传导结果都与解析解吻合良好,验证了本文 CV-FVM 对不同边界条件的适用性。为验证 CV-FVM 与 FEM 的计算效率,基于 Intel Core i7-10750H CPU,内存为 16 GB 的计算机,将算例 2 的网格尺寸加密至 210-5 m,分别采用 CV-FVM 与

29、FEM 再次计算,求解耗时及最大内存消耗见表 4。可以看出目前在求解时间上 CV-FVM 比FEM 长,但在内存消耗上要比 FEM 少,这是由于当前版本的软件还未在求解速度上进行优化,在后续研究中会逐步改进。523哈 尔 滨 工 程 大 学 学 报第 44 卷表 2 第 2 种情况 CV-FVM 结果对比Table 2 Comparison of CV-FVM results in the second case坐标/(m,m)CV-FVM 解/解析解/误差/%(510-3,0)530.66530.380.05(510-3,110-3)479.55479.380.04(510-3,210-3)

30、445.97445.970.00(510-3,310-3)427.56427.72-0.04(510-3,410-3)422.27422.56-0.07(510-3,510-3)428.77429.14-0.09(510-3,610-3)446.48446.88-0.09(510-3,710-3)475.38475.79-0.09(510-3,810-3)515.63516.05-0.08(510-3,910-3)566.82567.30-0.08(510-3,110-2)627.04627.73-0.11表 3 第 2 种情况 FEM 结果对比Table 3 Comparison of FE

31、M results in the second case坐标/(m,m)FEM 解/解析解/误差/%(510-3,0)530.79530.380.08(510-3,110-3)479.67479.380.06(510-3,210-3)445.09445.970.03(510-3,310-3)427.66427.72-0.01(510-3,410-3)422.37422.56-0.04(510-3,510-3)428.87429.14-0.06(510-3,610-3)446.56446.88-0.07(510-3,710-3)475.45475.79-0.07(510-3,810-3)515.

32、69516.05-0.07(510-3,910-3)566.86567.30-0.08(510-3,110-2)627.07627.73-0.10表 4 CV-FVM 与 FEM 计算效率对比Table 4Comparison of computational efficiency between CV-FVM and FEM数值方法求解时间/s最大内存消耗/GBCV-FVM71.5FEM31.92.3 CV-FVM 的优势 基于 2.2 节的算例 1,分别讨论 c=d=0,-200 m-1时温度的分布。图 6 为 c=d=0 m-1时的等温度线图。此时导热系数等于 k0,材料为均质材料;换热

33、系数等于 h0,边界条件均匀。因此,计算得到的温度关于 x=510-3 m 对称分布。图 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取

34、 y=110-2 m 处的温度分布结果做进一步分析。分别对比网格尺寸为 0.2510-3、110-3、210-3 m 时的 FEM 结果和本文 CV-FVM 结果,并根据公式算出 y=110-2 m 处的温度解析解。可以看到,随着网格尺寸的增加,FEM 结果在角点附近出现了不合理的数值波动现象,而本文 CV-FVM 始终能给出合理的结果。图 7 中这种等温线汇集在一个点上的现象在物理上不合理,数学上这种位置称为奇点,这是由角点处存在 2 个温度而造成的20。在处理这种情况时,由于 CV-FVM 中将边界变量存放至边界网格面中心,使用这些边界条件变量时只需按照边界序号依次循环取用,若该处存在 2

35、 个边界条件,则由循环顺序决定采用哪个边界。受网格精度影响,当网格尺寸较大,此处的温度阶跃变化表现为一个线性变化(见图 8)。相比于 FEM,CV-FVM处理此类问题时,计算所得的结果整体可信度更高,数值更加稳定。在处理复杂热传导问题时,难以提前预知大温度梯度存在的区域,采用本文 CV-FVM,即使不加密该区域的网格,也能得到合理的计算结果。623第 2 期龚京风,等:二维功能梯度材料热传导格点型有限体积法研究图 8 y=110-2 m 处的温度曲线Fig.8 Temperature curve when y=110-2 m3 结论 1)本文发展的用于二维 FGM 热传导问题的CV-FVM 求

36、解方法,能够处理线性三角形单元、双线性四边形单元及混合网格,适用于不同的边界条件类型。2)由于采用交错网格技术考虑物性参数的空间变化,没有基于任何分布假设,该方法适用于具有任意材料物性分布的二维 FGM 热传导问题。3)在处理边界角点存在剧烈温差情况时,FEM需要通过加密网格来避免数值振荡,而本文发展的CV-FVM 即使在粗糙网格情况下,也能给出合理的数值结果。因此,本文发展的 CV-FVM 更适用于恶劣环境下、物性参数复杂变化的 FGM 热传导问题,尤其对于无法提前预知大温度梯度区域的问题,本文方法具有明显优势。参考文献:1 刘文光,丰霞瑶,姚婉.功能梯度圆柱壳的热应力与热传导分析J.航空材

37、料学报,2019,39(6):81-89.LIU Wenguang,FENG Xiayao,YAO Wan.Analysis on thermal stress and heat conduction of functionally graded cylindrical shell J.Journal of aeronautical materials,2019,39(6):81-89.2 FU Tao,WU Xing,XIAO Zhengming,et al.Ther-moacoustic response of porous FGM cylindrical shell sur-round

38、by elastic foundation subjected to nonlinear thermal loading J.Thin-walled structures,2020,156:106996.3 许新,李世荣.功能梯度材料微梁的热弹性阻尼研究J.力学学报,2017,49(2):308-316.XU Xin,LI Shirong.Analysis of thermoelastic damping for functionally graded material micro-beam J.Chinese journal of theoretical and applied mechan

39、ics,2017,49(2):308-316.4 梅靖,陈鼎,杨博.石墨烯增强功能梯度材料板条的热弹性响应J.应用力学学报,2020,37(5):1857-1863,2309.MEI Jing,CHEN Ding,YANG Bo.Thermo-elastic re-sponse of functionally graded graphene reinforced stripsJ.Chinese journal of applied mechanics,2020,37(5):1857-1863,2309.5 STEINBERG M A.Materials for aerospaceJ.Scien

40、tific American,1986,255(4):66-72.6 SAINI R,LAL R.Axisymmetric vibrations of temperature-dependent functionally graded moderately thick circular plates with two-dimensional material and temperature distri-butionJ.Engineering with computers,2022,38(1):437-452.7 TANG Ye,DING Qian.Nonlinear vibration an

41、alysis of a bi-directional functionally graded beam under hygro-thermal loadsJ.Composite structures,2019,225:111076.8 许杨健,白雪,刘硕.复杂边界条件下二维功能梯度结构温度场解析解J.燕山大学学报,2018,42(6):525-531.XU Yangjian,BAI Xue,LIU Shuo.Analytical solution of temperature field of two dimensional functionally graded structure unde

42、r complex boundary conditionsJ.Journal of Yanshan University,2018,42(6):525-531.9 蒋科.基于混合数值法的功能梯度结构热传导研究D.西安:西安工业大学,2018:15-22.JIANG Ke.Heat conduction investigation of the functionally graded structure based on hybrid numerical method D.Xian:Xian Technological University,2018:15-22.10 SLADEK J,SLAD

43、EK V,HELLMICH C,et al.Heat con-duction analysis of 3-D axisymmetric and anisotropic FGM bodies by meshless local petrov-Galerkin method J.Computational mechanics,2007,39(3):323-333.11 CHING H K,YEN S C.Transient thermoelastic deforma-tions of 2-D functionally graded beams under nonuniformly convecti

44、ve heat supply J.Composite structures,2006,73(4):381-393.12 仝国军,许杨健.非均匀温度场下变物性二维功能梯度材料板的瞬态热应力分析J.机械科学与技术,723哈 尔 滨 工 程 大 学 学 报第 44 卷2018,37(6):980-984.TONG Guojun,XU Yangjian.Transient thermal stress a-nalysis of variable properties for two-dimensional function-ally graded material plates under non-u

45、niform temperature field J.Mechanical science and technology for aero-space engineering,2018,37(6):980-984.13 GONG Jingfeng,XUAN Lingkuan,MING Pingjian,et al.An unstructured finite-volume method for transient heat conduction analysis of multilayer functionally graded materials with mixed grids J.Num

46、erical heat transfer,part B:fundamentals,2013,63(3):222-247.14 GONG Jingfeng,XUAN Lingkuan,YING Baosheng,et al.Thermoelastic analysis of functionally graded porous materials with temperature-dependent properties by a stag-gered finite volume method J.Composite structures,2019,224:111071.15 龚京风.热障涂层活

47、塞热应力的有限体积方法研究与应用D.哈尔滨:哈尔滨工程大学,2014:34-38.GONG Jingfeng.Research and application of the finite volume method for the thermal barrier coated piston D.Harbin:Harbin Engineering University,2014:34-38.16 宣领宽.结构振动声辐射的有限体积法研究与应用D.哈尔滨:哈尔滨工程大学,2014:30-35.XUAN Lingkuan.Research and application of the FVM in t

48、he structural vibration and sound radiation D.Harbin:Harbin Engineering University,2014:30-35.17 孔祥谦.热应力有限单元法分析M.上海:上海交通大学出版社,1999.18 CHAROENSUK J,VESSAKOSOL P.A high order con-trol volume finite element procedure for transient heat con-duction analysis of functionally graded materials J.Heat and ma

49、ss transfer,2010,46(11/12):1261-1276.19 许杨健,李祥瑞,杨秋足,等.对流换热边界下 2D-FGM 板稳态温度场精确解J.北京工业大学学报,2018,44(10):1284-1290.XU Yangjian,LI Xiangrui,YANG Qiuzu,et al.Accu-rate solution of steady state temperature field of 2D-FGM plate under convective heat transfer boundary J.Jour-nal of Beijing University of Tec

50、hnology,2018,44(10):1284-1290.20 陶文铨.传热学M.5 版.北京:高等教育出版社,2019.本文引用格式:龚京风,徐宗著,宣领宽,等.二维功能梯度材料热传导格点型有限体积法研究J.哈尔滨工程大学学报,2023,44(2):321-328.GONG Jingfeng,XU Zongzhu,XUAN Lingkuan,et al.Investigation on Cell-vertex finite volume method for the heat conduction of two-dimensional func-tionally graded materi