熱障涂層活塞熱應(yīng)力分析的格點(diǎn)型有限體積法

龔京風(fēng), 宣領(lǐng)寬, 周少偉, 明平劍

(1. 中國艦船研究設(shè)計(jì)中心, 武漢 430064; 2. 哈爾濱工程大學(xué) 動力與能源工程學(xué)院, 哈爾濱 150001)

熱障涂層活塞熱應(yīng)力分析的格點(diǎn)型有限體積法

龔京風(fēng)1, 宣領(lǐng)寬1, 周少偉1, 明平劍2

(1. 中國艦船研究設(shè)計(jì)中心, 武漢 430064; 2. 哈爾濱工程大學(xué) 動力與能源工程學(xué)院, 哈爾濱 150001)

摘要:為研究熱障涂層(thermal barrier coating, TBC)活塞熱應(yīng)力問題,將二維格點(diǎn)型有限體積法(cell vertex FVM, CV-FVM)推廣用于三維復(fù)合材料熱力性能研究．利用交錯(cuò)網(wǎng)格技術(shù),通過將物性參數(shù)和待解變量分別定義在單元中心和節(jié)點(diǎn)上,將物性參數(shù)的空間變化引入離散過程,從而避免數(shù)值不連續(xù)問題．采用CV-FVM數(shù)值模擬普通活塞熱應(yīng)力場,計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果吻合良好．?dāng)?shù)值模擬TBC活塞熱應(yīng)力場,分析活塞TBC區(qū)域的熱力性能,結(jié)果表明:活塞的最大溫度位于燃燒室與活塞頂面的交界區(qū)域;應(yīng)力最大值位于粘結(jié)層與陶瓷層的分界面附近;周向應(yīng)力對活塞的整體應(yīng)力水平起主要影響．?dāng)?shù)值結(jié)果表明CV-FVM能夠作為TBC活塞熱力性能研究的數(shù)值預(yù)測工具．

關(guān)鍵詞:交錯(cuò)網(wǎng)格技術(shù);格點(diǎn)型有限體積法;熱障涂層;活塞;熱應(yīng)力

熱障涂層(thermal barrier coating, TBC)活塞有助于改善內(nèi)燃機(jī)的燃燒性能、排放性能、抗爆震性能等，TBC最重要的作用是為活塞基體提供熱防護(hù)，減少傳入活塞基體的熱量. 而TBC應(yīng)當(dāng)采用何種形式，以及TBC對于活塞熱力性能的改變是設(shè)計(jì)階段需要研究的重要問題.

建立適用于熱障涂層活塞熱力性能研究的數(shù)值方法,作為TBC活塞設(shè)計(jì)階段的預(yù)測工具,分析TBC活塞溫度場、熱應(yīng)力場、熱變形具有重要意義．

目前,關(guān)于TBC活塞熱力性能的數(shù)值仿真相對較少,一般采用的是基于FEM的商用軟件．王素等[1-3]采用ADINA數(shù)值分析復(fù)合材料活塞的整體性能,指出使用陶瓷纖維梯度層可以明顯改變活塞溫度分布,緩和由于熱膨脹系數(shù)不匹配,在陶瓷纖維增強(qiáng)層與活塞本體交界處產(chǎn)生的應(yīng)力,但沒有詳細(xì)討論活塞表面熱應(yīng)力性能．Buyukkaya等[4-5]在活塞上表面全部涂覆TBC,采用ANSYS數(shù)值模擬其溫度場,指出普通活塞的最大溫度出現(xiàn)在燃燒室中心,而TBC活塞的最大溫度在燃燒室與活塞頂面交界處．Cerit等[6-7]在活塞頂面局部涂覆TBC,采用ANSYS數(shù)值模擬其溫度場及熱應(yīng)力場,給出了溫度、應(yīng)力在各材料交界面處的分布曲線．由于TBC活塞在基體表面有涂層材料,相對于傳統(tǒng)活塞數(shù)值分析,TBC活塞數(shù)值分析需要額外考慮如何劃分TBC區(qū)域網(wǎng)格和如何處理非均勻材料．Aboudi等[8]指出，為了獲得正確的復(fù)合材料熱應(yīng)力場,在控制方程的離散過程中需要引入物性參數(shù)的空間變化．傳統(tǒng)的FEM賦予每一個(gè)網(wǎng)格單元均勻的物性參數(shù),未考慮材料屬性的空間變化,這一處理可能導(dǎo)致計(jì)算結(jié)果存在人工數(shù)值不連續(xù)問題[9-11]．

龔京風(fēng)等[12-15]發(fā)展了一種適用于二維復(fù)合材料熱應(yīng)力研究的格點(diǎn)型有限體積方法(cell vertex FVM, CV-FVM),數(shù)值結(jié)果表明,CV-FVM能夠有效避免物性參數(shù)引起的數(shù)值不連續(xù)問題．本文進(jìn)一步將CV-FVM推廣應(yīng)用于三維問題,通過數(shù)值算例驗(yàn)證計(jì)算方法的正確性,并基于CV-FVM分析TBC活塞熱力性能．

1基本方程

考慮穩(wěn)態(tài)情況下各向同性線彈性復(fù)合材料的熱應(yīng)力問題,材料的物性參數(shù)隨空間變化．基于控制體Ω建立熱傳導(dǎo)和熱彈性方程,Ω體積為V,邊界面為S．1.1熱傳導(dǎo)方程

熱傳導(dǎo)控制方程為

(1)

式中:ρ、c、T、k分別為密度、比熱容、溫度、熱傳導(dǎo)系數(shù); nα為微元邊界S的單位外法線矢量分量,α=x, y, z; qVT為熱源在單位時(shí)間單位體積內(nèi)產(chǎn)生的熱量,本文不涉及熱源,后文中將省略相應(yīng)內(nèi)容．考慮3種邊界條件:

(2)

(3)

(4)

式中:TB為Dirichlet邊界SD上給定的溫度;qB為Neumann邊界SN上給定的法向熱流量;hB為Robin邊界SR上給定的對流換熱系數(shù),T∞為環(huán)境溫度．下標(biāo)B代表邊界．

1.2熱彈性方程

熱彈性控制方程為

式中: uα為位移矢量的分量,σαβ為應(yīng)力張量的分量,fα為體積力矢量的分量．本文不涉及體積力,后文中將省略相應(yīng)內(nèi)容．

利用線彈性體的本構(gòu)方程及幾何方程,熱彈性方程可寫為

(5)

式中: Γ為熱彈性系數(shù); G、λ為拉姆系數(shù); a為線膨脹系數(shù);Tr為參考溫度,當(dāng)T=Tr時(shí)結(jié)構(gòu)的熱應(yīng)變?yōu)?;δαβ為克羅尼克爾符號,當(dāng)α=β時(shí)δαβ=1,當(dāng)α ≠ β時(shí)δαβ=0．

考慮邊界條件:

(6)

式中: uαB為邊界SD上給定的位移;σnB為邊界SN上給定的法向應(yīng)力．

2數(shù)值離散

對于本文考慮的活塞熱應(yīng)力問題,數(shù)值方法應(yīng)能夠處理三維非結(jié)構(gòu)化網(wǎng)格．圍繞網(wǎng)格節(jié)點(diǎn)依次連接相鄰單元中心和面中心建立控制體,見圖1．空心圓點(diǎn)代表網(wǎng)格節(jié)點(diǎn),實(shí)心圓點(diǎn)代表體中心及面中心．用粗實(shí)線及虛線圍成網(wǎng)格,細(xì)實(shí)線及點(diǎn)劃線圍成控制體．關(guān)于二維問題的數(shù)值離散參考文獻(xiàn)[12-15]．

本文采用交錯(cuò)網(wǎng)格技術(shù)模擬TBC活塞物性參數(shù)的空間變化．在網(wǎng)格節(jié)點(diǎn)上定義待解變量(溫度和位移),并假設(shè)在控制體內(nèi)均勻分布;將已知量(如物性參數(shù))定義在單元中心,并假設(shè)在網(wǎng)格單元內(nèi)均勻分布．考慮控制體與網(wǎng)格單元的關(guān)系,可知控制體內(nèi)物性參數(shù)是變化的,從而在離散過程中考慮空間變化的物性參數(shù)．

圖1　三維網(wǎng)格示意

2.1熱傳導(dǎo)方程的離散

利用向后差分格式離散熱傳導(dǎo)方程(1)左端的一階時(shí)間項(xiàng),得

式中: t為當(dāng)前時(shí)刻,Δt為時(shí)間步長,Vi為當(dāng)前節(jié)點(diǎn)周圍第i個(gè)單元的體積,nc為節(jié)點(diǎn)周圍單元的總數(shù),ncni為第i個(gè)單元的節(jié)點(diǎn)數(shù)．下標(biāo)i代表第i個(gè)單元中心的變量值．

參考FEM的思想,利用型函數(shù)離散方程(1)右端第一項(xiàng)

式中:N為型函數(shù),下標(biāo)ij代表第i個(gè)單元中第j個(gè)節(jié)點(diǎn)上的變量值．

考慮邊界條件,對于邊界SD上的節(jié)點(diǎn),其溫度直接根據(jù)式(2)給定．對于邊界SN和SR上的節(jié)點(diǎn),將式(3)和(4)代入,得到熱傳導(dǎo)方程的最終離散形式:

式中:ABi為與節(jié)點(diǎn)相鄰的第i個(gè)邊界單元面的面積矢量．nN為與當(dāng)前節(jié)點(diǎn)相鄰的位于SN上的面單元數(shù),nR為與當(dāng)前節(jié)點(diǎn)相鄰的位于SR上的面單元數(shù)．

2.2熱彈性方程的離散

利用中心差分格式離散熱彈性方程(5)左端的二階時(shí)間項(xiàng),得

參考熱傳導(dǎo)方程空間項(xiàng)的離散方法,離散熱彈性方程(5)右端第一項(xiàng)和第二項(xiàng)

考慮邊界條件,對于SD上的節(jié)點(diǎn),位移給定為uαB．對于SN上的點(diǎn),將式(6)代入得到熱彈性方程的最終離散形式

2.3數(shù)值方法的實(shí)施

采用Fortran語言編程實(shí)現(xiàn)本文的數(shù)值方法．文獻(xiàn)[16-18]采用分離解法求解不同方向的位移方程,即在求解某一方向的位移時(shí),相關(guān)的其它位移采用上一次的計(jì)算結(jié)果,迭代求解各方向位移,直到獲得收斂解．區(qū)別于上述文獻(xiàn),本文采用耦合解法同時(shí)求解3個(gè)方向位移,從而1次計(jì)算即可獲得位移場．不同網(wǎng)格單元的型函數(shù)表達(dá)式見文獻(xiàn)[18]．

3活塞計(jì)算模型

3.1活塞幾何模型及網(wǎng)格模型

模擬活塞熱應(yīng)力問題,考慮到活塞幾何形狀的對稱性,取如圖2所示的1/4模型,活塞直徑為170mm．活塞模型1與模型2的區(qū)別在于,活塞模型1上表面具有避閥坑,沒有布置TBC,通過冷卻油道冷卻活塞頭部;而活塞模型2上表面被簡化了,沒有避閥坑,布置有TBC,沒有冷卻油道,利用TBC的隔熱作用保護(hù)活塞頭部．采用四面體單元和棱柱單元劃分計(jì)算域,見圖2．

(a) 活塞1網(wǎng)格模型　　　　　 (b) 活塞2網(wǎng)格模型

3.2邊界條件及物性參數(shù)

以活塞2為例給出計(jì)算邊界條件．活塞x=0及z=0(見圖3)平面為對稱面,設(shè)為絕熱邊界,給定簡支位移約束．為了避免活塞的剛體平移,活塞底面給定簡支約束邊界條件．參考文獻(xiàn)[19]中的活塞表面對流換熱系數(shù)和環(huán)境溫度(由技術(shù)參數(shù)和經(jīng)驗(yàn)公式得到),取平均值做穩(wěn)態(tài)計(jì)算．活塞熱邊界條件見圖3和表1．活塞物性參數(shù)見表2．

圖3　活塞邊界

邊界對流換熱系數(shù)/(W·m-2·℃-1)環(huán)境溫度/℃活塞上表面(頂面及燃燒室)600687 火力岸50427上壁1400427 第一道環(huán)內(nèi)側(cè)80307下壁450227 第一環(huán)岸100127上壁1985127 第二道環(huán)內(nèi)側(cè)100127下壁3400127 第二環(huán)岸200127上壁1435112 第三道環(huán)內(nèi)側(cè)200112下壁1435112 冷卻油道460076.8 裙部80076.8 活塞銷孔50076.8上部60076.8 活塞內(nèi)腔中部50076.8下部40076.8

表2　活塞物性參數(shù)

4活塞熱應(yīng)力分析

4.1普通活塞熱應(yīng)力

為驗(yàn)證本文數(shù)值方法的正確性,數(shù)值模擬普通活塞1的熱應(yīng)力場,活塞材料為硅鋁合金．圖4為不同方法計(jì)算得到的45°平面上位移ux云圖．圖5為不同方法計(jì)算的45°平面r=0.072m處的應(yīng)力σxx曲線,其中r=0.072m位于冷卻油腔和活塞環(huán)槽之間,熱應(yīng)力場較為復(fù)雜．通過對比可以看出,本文發(fā)展的三維CV-FVM計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果一致,說明本文三維CV-FVM可用于活塞熱應(yīng)力數(shù)值模擬．

(a) ANSYS計(jì)算結(jié)果　　　　 (b) CV-FVM計(jì)算結(jié)果

圖5　活塞1的應(yīng)力σxx曲線對比(45°平面上r=0.072 m處)4.2　TBC活塞熱應(yīng)力分析

在活塞頂面涂覆1 mm的TBC．活塞基體為鋁硅合金ZL 109,粘結(jié)材料為NiCrAl(厚度為0.2 mm),陶瓷材料為MgZrO3(厚度為0.8 mm),材料物性參數(shù)見表2．圖6給出了活塞的表面熱流矢量及溫度分布云圖．活塞2表面的TBC承受了大部分的熱載荷,熱流量集中于TBC區(qū)域,活塞上表面溫度明顯高于活塞基體溫度．

圖6　活塞2的表面熱流矢量及溫度云圖

圖7為活塞2內(nèi)θ=45°平面上的熱應(yīng)力場云圖．由圖可見,活塞的最大溫度T位于燃燒室與活塞頂面的交界區(qū)域(A區(qū)域),主要原因是A區(qū)為活塞上表面轉(zhuǎn)折位置,相較于其它位置熱阻更大,熱量傳遞更少,因而溫度越高．軸向應(yīng)力σyy、徑向應(yīng)力σr和周向應(yīng)力σθ幅值相當(dāng),其最大值位于粘結(jié)層與陶瓷層的分界面附近,主要原因是陶瓷層阻隔了大部分的熱量,相對而言傳遞到粘結(jié)層、金屬基體的熱量較小,故陶瓷層和粘結(jié)層之間的溫度梯度比粘結(jié)層和金屬基體間的溫度梯度大得多,陶瓷層和粘結(jié)層間存在更大的熱應(yīng)力．由于TBC內(nèi)物性參數(shù)的突變,TBC區(qū)域存在明顯的應(yīng)力集中現(xiàn)象,容易發(fā)生破壞．

圖7　活塞2的θ=45°平面上熱應(yīng)力場云圖

分析TBC活塞2上表面θ=45°處的溫度及應(yīng)力曲線,見圖8．

(a) 溫度曲線

(b) 應(yīng)力曲線

燃燒室頂面包含3個(gè)區(qū)域:區(qū)域1、區(qū)域2位于燃燒室,區(qū)域3是活塞頂面．燃燒室邊緣區(qū)域用斜線標(biāo)出．溫度最小值位于區(qū)域2,最大值位于B點(diǎn)．參考溫度曲線,可以看到在燃燒室邊緣區(qū)域溫度梯度最大,應(yīng)力在該區(qū)域變化幅值相對其他區(qū)域大得多．與應(yīng)力σr和σyy相比,應(yīng)力σθ沿徑向變化平緩且幅值最大．Von Mises應(yīng)力σvm可代表活塞的綜合應(yīng)力情況,其變化趨勢與應(yīng)力σθ基本一致,說明應(yīng)力σθ對活塞的整體應(yīng)力水平起主要影響．

5結(jié)論

1) 本文基于CV-FVM求解熱傳導(dǎo)方程和熱彈性方程,將適用于二維復(fù)合材料熱應(yīng)力問題的CV-FVM推廣應(yīng)用于三維熱應(yīng)力問題．采用CV-FVM分析普通活塞的熱應(yīng)力問題,計(jì)算結(jié)果與商用軟件ANSYS計(jì)算結(jié)果吻合良好,表明該方法的正確性．

2) 基于CV-FVM進(jìn)一步分析TBC活塞熱應(yīng)力場．通過分析活塞熱應(yīng)力場可知,陶瓷層能夠阻隔活塞頂面的大部分熱量,使得粘結(jié)層和陶瓷層間存在較大溫度梯度,引起應(yīng)力集中,應(yīng)力最大值出現(xiàn)在粘結(jié)層和陶瓷層分界面區(qū)域．

3) 通過分析TBC活塞熱應(yīng)力曲線可知,溫度最大值出現(xiàn)在燃燒室與活塞頂面的交界區(qū)域,活塞的應(yīng)力水平主要受周向應(yīng)力σθ影響．

4) TBC活塞熱應(yīng)力場沒有出現(xiàn)數(shù)值不連續(xù)現(xiàn)象,CV-FVM能夠作為TBC活塞熱力性能研究的數(shù)值預(yù)測工具．

參考文獻(xiàn)

[1] 王素, 倪春陽, 朱心雄. 功能梯度材料活塞三維溫度場分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2005, 31(12): 1299-1302.

[2] 王素, 倪春陽, 朱心雄. 一種功能梯度材料活塞的材料梯度方程選擇[J]. 機(jī)械科學(xué)與技術(shù), 2006, 25(2): 133-138.

[3] 王素, 倪春陽, 朱玉明, 等. 陶瓷纖維梯度增強(qiáng)活塞的梯度方程研究[J]. 航空學(xué)報(bào), 2007, 28(1): 234-239.

[4] BUYUKKAYA E, CERIT M. Thermal analysis of a ceramic coating diesel engine piston using 3-D finite element method[J]. Surface and Coatings Technology, 2007, 202(2): 398-402.

[5] BUYUKKAYA E. Thermal analysis of functionally graded coating AlSi alloy and steel pistons[J]. Surface and Coatings Technology, 2008, 202(16): 3856-3865.

[6] CERIT M. Thermo mechanical analysis of a partially ceramic coated piston used in an SI engine[J]. Surface and Coatings Technology, 2011, 205(11): 3499-3505.

[7] CERIT M, AYHAN V, PARLAK A, et al. Thermal analysis of a partially ceramic coated piston: Effect on cold start HC emission in a spark ignition engine[J]. Applied Thermal Engineering, 2011, 31(2/3): 336-341.

[8] ABOUDI J, PINDERA M J, ARNOLD S M. Higher-order theory for functionally graded materials[J]. Composites, Part B, 1999, 33: 777-832.

[9] SANTARE M H, LAMBROS J. Use of graded finite elements to model the behavior of nonhomogeneous materials[J]. ASME Journal of Applied Mechanics, 2000, 67(4): 819-822.

[10]KIM J H, PAULINO G H. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials[J]. ASME Journal of Applied Mechanics, 2002, 69(4): 502-514.

[11]CAVALCANTE M A A, MARQUES S P C, PINDERA M J. Parametric formulation of the finite-volume theory for functionally graded materials-Part II: numerical results [J]. Journal of Applied Mechanics, 2007, 74(5): 946-957.[12]GONG J F, XUAN L K, MING P J, et al. Thermoelastic analysis of functionally graded solids using a staggered finite volume method[J]. Composite Structures, 2013, 104: 134-143.

[13]GONG J F, XUAN L K, MING P J, et al. An unstructured finite-volume method for transient heat conduction analysis of multilayer functionally graded materials with mixed grids[J]. Numerical Heat Transfer: Part B, 2013, 63(3): 222-247.

[14]GONG J F, XUAN L K, MING P J, et al. Application of the staggered cell-vertex finite volume method to thermoelastic analysis in heterogeneous materials[J]. Journal of Thermal Stresses, 2014, 37(4): 506-531.

[15]龔京風(fēng), 明平劍, 宣領(lǐng)寬, 等. 基于有限體積法求解FGM動態(tài)響應(yīng)及固有特性[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2013, 41(5): 28-33.

[16]TAYLOR G A, BAILEY C, CROSS M. Solution of the elastic/visco-plastic constitutive equations: a finite volume approach[J]. Applied Mathematical Modelling, 1995, 19(12): 746-760.

[17]BAILEY C, CROSS M. A finite volume procedure to solve elastic solid mechanics problems in three dimensions on an unstructured mesh [J]. International Journal for Numerical Methods in Engineering, 1995, 38(10): 1757-1776.

[18]TAYLOR G A, BAILEY C, CROSS M. A vertex-based finite volume method applied to non-linear material problems in computational solid mechanics[J]. International Journal for Numerical Methods in Engineering, 2003, 56(4): 507-529.

[19]劉友. 船用柴油機(jī)活塞瞬態(tài)溫度場測試與分析研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2013.

(編輯楊波)

Cell vertex FVM for thermoelastic analysis of the piston with thermal barrier coating

GONG Jingfeng1, XUAN Lingkuan1, ZHOU Shaowei1, MING Pingjian2

(1. China Ship Development and Design Center, Wuhan 430064, China;2. College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

Abstract:To analyze the performance of the piston with thermal barrier coating (TBC), two dimensional cell vertex FVM (CV-FVM) has been extended to the research of the thermoelastic performance of three-dimensional composite structures. The staggered grid technique was adopted to incorporate property variation into the course of solution so as to avoid numerical discontinuity. The unknown variable was defined at the cell vertex, while the material property was defined at the cell center. The developed CV-FVM was used to simulate the thermoelastic fields of the ordinary piston and the predicted results agreed with the ANSYS results. Then, the CV-FVM was used to simulate the thermoelastic fields of the TBC piston and the performance were analyzed. The maximum temperature is in the area between the combustion chamber and the top surface, and the maximum stress is around the interface between the ceramic layer and the cohesive layer. The circumferential stress contributes the main part of the thermoelastic stress in the TBC area. The successful application demonstrates that the CV-FVM can be used as a predictive tool for the thermoelastic analysis of the TBC piston.

Keywords:staggered grid technique; cell vertex FVM; thermal barrier coating; piston; thermoelastic

doi：10.11918/j.issn.0367-6234.2016.07.012

收稿日期:2015-05-03

基金項(xiàng)目:國家自然科學(xué)基金(51206031)

作者簡介:龔京風(fēng)(1986—),女,博士,工程師

通信作者:宣領(lǐng)寬, xuanlingkuan@163.com

中圖分類號:TK422

文獻(xiàn)標(biāo)志碼:A

文章編號:0367-6234(2016)07-0076-06