熱障涂層活塞熱應力高階控制體有限元法研究

劉 琦， 杜雁霞， 齊文亮， 劉 磊， 肖光明

(1. 中國空氣動力研究與發展中心空氣動力學國家重點實驗室， 綿陽 621000； 2. 哈爾濱工程大學動力與能源工程學院， 哈爾濱 150001； 3. 中國航空工業集團公司西安航空計算技術研究所， 西安 710065)

1 引 言

熱障涂層(thermal barrier coating， TBC)活塞是在傳統活塞頂端涂覆陶瓷等離子體的一種新型活塞. 陶瓷材料導熱系數小，可有效減少活塞頂端熱量向活塞內部的傳遞，降低活塞基體溫度的同時提高活塞工作溫度，改善內燃機燃燒性能、燃油消耗和排放. TBC熱力失效問題是TBC活塞研究中的一個重要問題. TBC材料的組成形式以及TBC與基體材料熱匹配特性研究已成為TBC活塞設計評估階段的一個重要方面[1-2]. 建立適用于TBC活塞熱學及力學特性研究的數值方法, 掌握TBC活塞服役過程中的熱力分布規律，對實現TBC活塞結構優化設計具有重要意義.

目前，針對TBC活塞熱力特性的數值研究還相對較少，主要是借助商業FEM軟件進行開展. 王素等[3]基于ADINA有限元分析軟件對陶瓷纖維梯度增強活塞的熱力分布進行了計算，發現陶瓷纖維梯度層可以有效改善纖維與基體層間因熱匹配而產生的應力集中問題，但未考慮粘結層對熱力分布的影響. 牛小強等[4]基于ANSYS分析了層和型TBC對活塞的溫度、應力場的影響，指出TBC可以有效降低活塞頭部和環槽區域的工作溫度，活塞基體頂面、喉口區域及邊緣棱角處熱應力集中明顯. Buyukkaya等[5-6]采用ANASYS研究了TBC對鋁合金活塞和鋼制活塞溫度場的影響，結果表明普通活塞最大溫度出現在燃燒室中心，而TBC活塞的最大溫度出現在燃燒室與活塞頂面交界處，TBC可使鋁合金活塞溫度下降48%，使鋼制活塞溫度下降35%. Hejwowski[7]采用試驗方法對比研究了不同厚度的TBC對柴油機活塞傳熱、疲勞性能的影響，指出黏結層厚度為0.15 mm、陶瓷層厚度為0.35 mm時較為合適. 龔京風等[8-9]基于低階四面體及三棱柱單元發展了一種適用于TBC活塞熱應力計算的三維控制體有限元法(Control volume finite element method, CV-FEM)，該方法可以有效避免物性參數引起的數值不連續問題. 基于CV-FEM，龔京風[8]對比研究了層和型和功能型TBC活塞熱應力問題，發現相比于層和型TBC活塞，應用功能型TBC后粘結層附近應力集中問題并未得到明顯改善.

劉琦等[10-13]基于高階二維及三維單元建立了高階CV-FEM，結果表明高階CV-FEM既可以有效避免物性參數引起的數值不連續問題，又能在不增加求解自由度下提高計算精度. 本文將該方法應用到TBC活塞的熱應力問題分析中，研究不同材料排列形式下的TBC對活塞熱力特性的影響.

2 控制方程

本文基于任意控制體Ω建立熱傳導方程和平衡方程，研究穩態下各向同性線彈性復合材料熱應力問題. 其中，控制體Ω的體積為V，控制體的邊界面為S.

2.1 熱傳導方程

無內熱源積分格式穩態熱傳導方程為：

(1)

其中，n為界面外法矢量;k為熱傳導系數;T為待解溫度.

考慮三種換熱邊界條件：

T=TB, 邊界SD

(2)

(3)

(4)

其中,α(=x，y，z)為邊界S的單位外法矢量n的分量；Dirichlet邊界SD上的給定溫度；qB為Neumann邊界SN上的給定的流進或流出的法向熱流；hB及T∞分別為Robin邊界SR上的對流換熱系數及環境溫度.

2.2 熱彈性方程

無體積力穩態熱彈性方程為：

∮Sσ·ndS=0

(5)

其中σ為Cauchy應力張量，對于各向同性線彈性體，其本構關系為：

σ=2με+λtr(ε)I-Γ(T-Tr)I

(6)

其中,ε表示Cauchy應變張量;Tr表示參考溫度;I為單位張量;Γ為熱彈性系數;μ、λ表示拉梅系數.平面應變假設下，Γ、μ、λ的表達式為：

(7)

其中,E表示楊氏模量;ν表示泊松比;a表示熱膨脹系數.

Cauchy應變張量表達式為：

(8)

考慮兩種力學邊界條件：

uα=uαB, 邊界SD

(9)

σαβnβ=σnB, 邊界SN

(10)

其中,uαB為Dirichlet邊界上給定的位移；σαB為Neumann邊界上給定的力載荷.

3 數值方法

采用高階10點四面體單元劃分計算域.依次連接大四面體1-2-3-4，小四面體1-5-8-7、2-5-9-6、3-6-7-10、4-10-8-9單元中心和面中點，形成以單元角點(n1-n4)和單元邊中點(n5-n6)共10個子控制體(見圖1).利用交錯網格技術，將材料屬性及待解溫度T、待解位移uα分別存儲在單元中心及單元節點.

圖1 10點四面體單元子控制

3.1 熱傳導方程離散

借鑒FEM，單元內任意變量φ采用形函數進行插值近似

(11)

其中,Ni為形函數;ns為單元內的節點數.對于10點四面體單元ns=10. 10點四面體單元的形函數表達式可見文獻[10].

將式(11)引入式(1)，對其進行離散得

(12)

其中,下角標i表示第i個單元，ij表示第i個單元的第j個節點. 式(12)中形函數Nij對全局坐標系xα的導數將其映射到標準單元進行計算，對其面積分采用中點積分公式進行近似計算，具體計算流程可見文獻[10].

引入邊界條件式(12)～式(4)，得熱傳導方程最終離散格式：

(13)

其中，ABi為與節點相鄰的第i個單元的邊界面面積矢量；nN為與節點相鄰的位于SN上的邊界面個數；nR為與節點相鄰的位于SR上的邊界面個數.

將式(6)和式(8)帶入式(5)，采用與式(12)相同的離散方法，式(5)可離散為：

(14)

引入邊界條件式(9)和式(10)，得熱彈性方程最終離散格式：

(15)

4 數值求解

采用Fortran語言實現本文計算方法，開始計算之前，我們將與網格信息相關的形函數導數、積分點位置、形函數導數積分作為幾何常數進行一次性存儲，減少計算量．對與臨近單元存在Dirichlet邊界的節點，我們采用置大數法將對其導熱系數矩陣或剛度矩陣相應的行進行處理.離散后的熱傳導及熱彈性方程,我們采用直接求解法進行求解.

5 數值算例

5.1 圓環熱應力問題算例驗證

為驗證發展高階CV-FEM對復雜換熱邊界條件下熱應力問題求解的正確性，我們首先驗證受非均勻熱流加熱的圓筒熱應力問題(見圖2(a))，并將計算結果與ANASYS結果進行對比.計算模型和網格模型如圖2(b)所示. 圓筒采用10 068個10點四面體單元進行劃分，包含16 478個節點. CV-FEM與ANSYS采用相同的網格進行數值模擬. 圓筒的上表面受不均勻熱源加熱，熱源強度沿徑向的變化規律如式(16)所示. 圓筒外表面(Rout=0.5 m)、圓筒內表面(Rin=0.1 m)以及圓筒底面(z=0)施加對流換熱邊界條件.對流換熱系數分別為hout=hin=530 W·m-2·℃-1,hdown=180 W·m-2·℃-1,環境溫度Tout=Tin=30 ℃，Tdown=60 ℃.圓筒下表面施加絕熱邊界. 材料的彈性模量E=70 GPa，泊松比ν=0.33，導熱系數k=156 W·m-1·℃-1, 線性膨脹系數a=2.3×10-5℃-1.

(16)

(a) 計算模型以及邊界條件

(b) 網格模型

圖3 位移ur沿圓筒徑向分布規律

圖4 應力σr沿圓筒徑向分布規律

圖3和圖4為圓筒中面徑向位移和徑向應力沿r方向變化曲線. 從圖3和圖4可以看出，基于相同高階單元網格模型，CV-FEM計算結果與FEM計算結果吻合良好. 該算例驗證了高階CV-FEM對三維復雜換熱下熱應力問題求解的正確性.

文獻[8]研究了層合型和功能型TBC活塞熱應力問題，對比研究了不同分布形式下涂層界面處應力分布，發現相比于層合型TBC活塞，功能梯度型TBC活塞熱應力集中現象并未得到改善，各涂層交界面處仍存在較大的應力集中問題. 本節對該問題做進一步研究，優化功能型TBC活塞涂層內應力集中問題. TBC活塞計算模型如圖5(a)所示，活塞直徑為170 mm，采用1/4模型進行分析，見圖5(b). 活塞上表面鋪設5層涂層，每層厚度為0.2 mm. 考慮層合型涂層和功能梯度型涂層兩種情況(見圖6)，采用六種排布形式進行數值模擬(見表1). 幾何模型采用66 935個10點四面體單元進行網格劃分(見圖7). 活塞基體材料為鋁硅合金，陶瓷材料為MgZrO3，粘結材料為NiCrAl，材料的物性參數見表2. 活塞x=0和y=0平面施加對稱邊界和絕熱邊界條件，底部施加簡支邊界條件，活塞其他位置換熱邊界如表3所示，參考溫度T0=27 ℃.

(a)活塞模型

(b) 1/4活塞模型及邊界條件

(a)層合型涂層

(b)梯度型涂層

(a)計算網格 (b)局部放大圖

表1 TBC活塞材料分布

表2 活塞物性參數

表3邊界條件[8]

梯度型涂層物性參數采用Mori-Tanaka方法計算[14]：

(17)

式中，下角標C表示陶瓷MgZrO3；下角標B表示粘結劑NiCrAl；VC表示過渡層內陶瓷體積分數.指數型表示為VC=(d/4)m，其中m表示材料參數，d表示無量綱距離，其定義見圖6(b). 楊氏模量E和泊松比υ可采用彈性系數關系式計算. 計算結果取值位置如圖8所示.

E=9Kμ/(3K+μ)

ν=(3K-2μ)/(6K+2μ)

(19)

(a)活塞徑向應力σr和周向應力σθ定義

(b) 數據取值位置

圖9給出了無涂層活塞、TBC活塞4和TBC活塞5的結果對比. 從圖9(a)可以看出，涂覆熱障涂層后的活塞表面溫度明顯增大，相比于無涂層活塞最高溫度，活塞4最高溫度增加了130.57 ℃，提升39.32%，活塞5最高溫度增加了81.21 ℃，提升14.8%. 活塞4最高溫度高于活塞5是因為陶瓷的導熱系數要明顯小于粘結劑，由此導致其隔熱效果更好. 活塞4和活塞5活塞頭表面溫度的提高增加了活塞的工作溫度，降低了燃燒室內熱量向活塞的傳遞. 從圖9(b～c)可以看出，涂覆熱障涂層后的活塞最大徑向位移減小，同時涂層內出現了明顯的應力集中現象. 為了更加清晰地描述活塞內應力分布情況，圖10～圖14給出了無涂層活塞、活塞4和活塞5在θ=45°平面內計算結果. 從中可以看出, 無論是層合型TBC亦或功能型TBC，粘結層附近均出現明顯應力集中現象，在過渡層(第2層～第4層)內引入功能型TBC并未改善應力集中問題. 此外，我們嘗試改變層合型TBC陶瓷/粘結劑組合份數和功能型TBC過渡層內材料的指數系數mT，雖然應力幅值有所改變, 但對應力集中問題的改善不大，在粘結層附近仍存在較明顯的應力集中. 層合型TBC和功能型TBC活塞內應力幅值見圖15(a～b).

無涂層活塞 TBC活塞(Case4) TBC活塞(Case5)考慮到活塞4和活塞5的應力集中出現在粘結層內，該問題應是由粘結層與基體材料參數不匹配所致. 因此，在層合型TBC及功能型TBC的粘結層(第5層)內引入FGM，采用材料屬性由基體材料到粘結劑材料過渡的方式緩解粘結層附近應力集中問題. 功能型TBC粘結層內基體體積分數VT和粘結劑體積分數VB一致，VT=VB=0.5，粘結層內材料系數mB=0.3. 圖15(c)給出了不同過渡層材料系數mT活塞6的最大應力分布. 從圖中明顯看出，相比于僅在過渡層(第2層～第4層)內引入FGM計算結果(見圖15(b))，不同系數mT下粘結層(第5層)在引入FGM后的應力幅值明顯降低，最大應力σθ、σr和σz分別下降了22.9%(316.62→243.93)、24.13%(272.13→206.45)和24.8%(270.27→203.17). 過渡層內FGM系數mT對應力幅值影響明顯，隨著mT的增加，應力幅值增大. 同時，在層合型TBC活塞3和活塞4中的粘結層引入FGM后應力幅值降低，應力集中問題明顯改善. 相比于粘結層不引入FGM，活塞3應力σθ、σr和σz分別下降了32.83%(303.39→203.77)和32.31%(261.04→176.7)、30.62%(259.07→179.72); 活塞4應力σθ、σr和σz下降了21.34%(294.02→231.27)、22.68%(256.66→198.43)和23.61%(256.4→195.85)(見圖15(d)). 粘結層引入FGM活塞3應力分布見圖15. 圖15(d)給出了活塞6過渡層內材料系數mT=0.1，粘結層材料系數mB=0.3，0.5和0.7時的最大應力計算結果. 從中可以看出，mB=0.3時最大應力幅值最小. 圖17給出了位置A和位置B處TBC活塞沿涂層厚度方向的應力分布. 從中可以看出，在層合型及功能型TBC的粘結層內引入FGM后，涂層內粘結層交界面處應力幅值降低，應力集中問題得到明顯改善.

(a)T/℃

(b)ur/mm

(c) σo/MPa

(d) σr/MPa

(e) σz/MPa

(a) T/℃

(b) σθ/MPa

(c)σr/MPa

(d)σz/MPa

(a) T/℃

(b) σθ/MPa

(c) σr/MPa

(d) σz/MPa

(a) T/℃

(b) σθ/MPa

(c) σr/MPa

(d) σz/MPa

(a) T/℃

(b) σθ/MPa

(c) σr/MPa

(d) σz/MPa

(a) T/℃

(b) σθ/MPa

(c) σr/MPa

(d) σz/MPa

(a) 不同陶瓷厚度

(b) FGM過渡層

(c) FGM (過渡層+粘結層)

(d) FGM (過渡層+粘結層)及4/3層陶瓷+FGM粘結層

(a) T/℃

(b) σθ/MPa

(c) σr/MPa

(d) σz/MPa

(a) A點σθ

(b) A點σr

(c) B點σθ

(d) B點σr

(e) B點σz

6 結 論

本文基于高階10點四面體CV-FEM求解熱傳導和熱彈性方程，將適用于三維復合材料傳熱問題的高階CV-FEM推廣到熱彈性問題求解. 我們將CV-FEM計算結果與ANSYS結果進行了對比驗證，對層合型TBC活塞、功能梯度型TBC活塞熱應力集中問題進行研究，發現當粘結層為均勻材料時，改變功能型TBC材料系數或層合型TBC陶瓷/粘結劑的分布形式, 并不能避免粘結層附近出現的應力集中現象，而在粘結層內引入材料屬性由粘結劑到基體過渡的FGM可以有效改善應力集中問題，僅在1層或2層粘結層內引入FGM應力幅值即可下降20%～30%.