二维功能梯度材料热传导格点型有限体积法研究_龚京风.pdf
《二维功能梯度材料热传导格点型有限体积法研究_龚京风.pdf》由会员分享,可在线阅读,更多相关《二维功能梯度材料热传导格点型有限体积法研究_龚京风.pdf(8页珍藏版)》请在文库网上搜索。
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 二维 功能 梯度 材料 热传导 格点型 有限 体积 研究 龚京风