龔京風(fēng)，張文平，宣領(lǐng)寬，明平劍，韓春旭

（1.武漢理工大學(xué)能源與動(dòng)力工程學(xué)院，湖北武漢430063；2.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，黑龍江哈爾濱150001）

功能梯度熱障涂層熱應(yīng)力的有限體積法研究

龔京風(fēng)1，2，張文平2，宣領(lǐng)寬2，明平劍2，韓春旭2

（1.武漢理工大學(xué)能源與動(dòng)力工程學(xué)院，湖北武漢430063；2.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，黑龍江哈爾濱150001）

為了研究功能梯度熱障涂層（TBC）熱應(yīng)力問(wèn)題，改進(jìn)了交錯(cuò)格點(diǎn)型有限體積方法（SCV-FVM）。運(yùn)用交錯(cuò)網(wǎng)格技術(shù)，待解變量定義在單元節(jié)點(diǎn)上，物性參數(shù)定義在單元中心上。計(jì)算功能梯度圓盤(pán)的應(yīng)力場(chǎng)，結(jié)果與理論解吻合良好，驗(yàn)證方法的正確性。通過(guò)與有限元結(jié)果及高階功能梯度理論（HOTFGM）結(jié)果對(duì)比發(fā)現(xiàn)，SCV-FVM可以有效避免物性參數(shù)變化引起的人工數(shù)值不連續(xù)。研究功能梯度TBC熱應(yīng)力問(wèn)題，分析涂層厚度變化對(duì)熱應(yīng)力場(chǎng)的影響，通過(guò)與ANSYS計(jì)算結(jié)果對(duì)比驗(yàn)證數(shù)值模擬的正確性。

熱障涂層；功能梯度；有限體積法；交錯(cuò)網(wǎng)格；應(yīng)力場(chǎng)

在內(nèi)燃機(jī)部件表面噴涂熱障涂層（TBC）是發(fā)展低散熱發(fā)動(dòng)機(jī)的重要手段。TBC通常由粘結(jié)底層和陶瓷頂層組成，由于材料屬性的突變，在分層界面上產(chǎn)生應(yīng)力集中。采用功能梯度材料（FGM），使材料屬性逐漸變化，可以有效緩解應(yīng)力集中問(wèn)題，同時(shí)增強(qiáng)材料粘結(jié)強(qiáng)度。準(zhǔn)確預(yù)測(cè)功能梯度TBC熱力性能是設(shè)計(jì)優(yōu)化涂層結(jié)構(gòu)和材料分布的前提。

目前，關(guān)于功能梯度TBC熱力性能的研究主要基于試驗(yàn)［1-5］和數(shù)值模擬［6-11］。其中，數(shù)值模擬主要利用有限元商用軟件，如ADINA［8］、ANSYS［7，10-11］。由于有限元方法（FEM）需要計(jì)算并存儲(chǔ)大型剛度矩陣，計(jì)算量和內(nèi)存消耗較大［12］。另外，TBC屬于非均勻材料，對(duì)物性參數(shù)的不合理處理可能產(chǎn)生人工數(shù)值不連續(xù)問(wèn)題［13-15］。為了克服這一問(wèn)題，Santare和Lambros［16］將梯度單元引入傳統(tǒng)FEM中，而Aboudi等人［17］提出高階功能梯度理論（HOTFGM）。2種方法及其發(fā)展形式均能夠有效減輕不連續(xù)問(wèn)題，但不能完全避免［10，13-20］。

本文旨在發(fā)展一種能有效預(yù)測(cè)功能梯度TBC熱力性能的數(shù)值方法。格點(diǎn)型有限體積法（CVFVM）兼有FVM的守恒性和FEM的靈活性，是FEM的一種有效補(bǔ)充方法。文獻(xiàn)［12］將交錯(cuò)網(wǎng)格技術(shù)引入CV-FVM提出SCV-FVM，采用時(shí)間推進(jìn)法解離散的平衡方程，因此該方法僅適用于動(dòng)力學(xué)問(wèn)題。本文進(jìn)一步發(fā)展SCV-FVM用于研究穩(wěn)態(tài)熱應(yīng)力問(wèn)題，利用隱格式耦合法求解離散的熱彈性方程。通過(guò)數(shù)值算例驗(yàn)證發(fā)展的SCV-FVM預(yù)測(cè)熱應(yīng)力場(chǎng)的準(zhǔn)確性和對(duì)人工不連續(xù)問(wèn)題的有效性，進(jìn)而研究功能梯度TBC熱應(yīng)力問(wèn)題。

1 熱應(yīng)力控制方程

考慮二維各向同性線彈性材料穩(wěn)態(tài)熱應(yīng)力問(wèn)題，基于計(jì)算域內(nèi)任意控制體微元S給出熱應(yīng)力基本方程及邊界條件。

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

不考慮內(nèi)熱源，穩(wěn)態(tài)熱傳導(dǎo)方程為

式中：k和T分別為熱傳導(dǎo)系數(shù)和溫度。考慮3種邊界條件

式中：TB為邊界LD上給定的溫度（Dirichlet邊界）；qB為邊界LN上給定的法向熱流密度（Neumann邊界）；hB和T∞分別為邊界LR上給定的對(duì)流換熱系數(shù)和環(huán)境溫度（Robin邊界）。nα（α＝x，y）為邊界單位外法線矢量n的分量。

1.2 平衡方程

忽略體積力，根據(jù)牛頓第二定律得到平衡方程：

式中：σ為應(yīng)力張量。本構(gòu)方程和幾何方程為

式中：σαβ，εαβ和uα分別為應(yīng)力張量σ、應(yīng)變張量ε、位移矢量u的分量。a為熱膨脹系數(shù)，Tr為參考溫度。G和λ為拉姆系數(shù)，Γ為熱系數(shù)，根據(jù)不同條件由楊氏模量E和泊松比μ計(jì)算得到：

平面應(yīng)變

平面應(yīng)力

考慮2種邊界條件

式中：uαB為邊界LD上給定的位移（Dirichlet邊界）；tαB為邊界LN上給定的法向應(yīng)力（Neumann邊界）。簡(jiǎn)支邊界可由上述2種邊界條件組合得到。

2 數(shù)值方法

本文基于CV-FVM發(fā)展適用于功能梯度TBC熱應(yīng)力問(wèn)題的數(shù)值方法。圍繞單元節(jié)點(diǎn)依次連接相鄰單元中心及邊長(zhǎng)中點(diǎn)建立控制體，如圖1所示。單元由實(shí)線圍成，其節(jié)點(diǎn)用空心圓點(diǎn)表示?？刂企w由虛線圍成，其中心及邊長(zhǎng)中點(diǎn)用實(shí)心圓點(diǎn)表示。

圖1 網(wǎng)格單元及控制體示意圖Fig.1 Sketch map of grids and control volumes

采用交錯(cuò)網(wǎng)格技術(shù)，從而將物性參數(shù)的空間變化引入離散過(guò)程。待解變量定義在單元節(jié)點(diǎn)上，并假設(shè)變量在控制體內(nèi)均勻分布。物性參數(shù)定義在單元中心并假設(shè)均勻分布，則物性參數(shù)在控制體內(nèi)是變化的，如圖1所示。

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

采用高斯公式將方程寫(xiě)為式中：L為控制體的邊。采用FEM的思想，利用型函數(shù)離散方程

式中：nc為當(dāng)前節(jié)點(diǎn)周?chē)膯卧倲?shù)，ncni為第i個(gè)單元內(nèi)的節(jié)點(diǎn)個(gè)數(shù)。三角形單元內(nèi)，ncni＝3；四邊形單元內(nèi)，ncni＝4。N為型函數(shù)。下標(biāo)i代表第i個(gè)單元中心的變量值，下標(biāo)ij代表第i個(gè)單元內(nèi)第j個(gè)節(jié)點(diǎn)上的變量值。

對(duì)于邊界上的節(jié)點(diǎn)，相應(yīng)的控制方程離散需要考慮邊界條件的影響。對(duì)于Dirichlet邊界上的節(jié)點(diǎn)，其溫度直接根據(jù)式給定。對(duì)于Neumann邊界和Robin邊界上的節(jié)點(diǎn)，將式（3）和（4）代入方程（13）得

2.2 平衡方程的離散

將式（6）和（7）代入方程（5）得利用高斯公式及型函數(shù)，方程（15）可離散寫(xiě)為

式中：Liα為積分線Li在α方向上的投影長(zhǎng)度。

對(duì)于Dirichlet邊界上的節(jié)點(diǎn)，其位移直接根據(jù)式（10）給定。對(duì)于Neumann邊界上的節(jié)點(diǎn)，平衡方程的積分利用式（11）得

離散方程（14）和（17）中的型函數(shù)積分及積分線長(zhǎng)度的計(jì)算可參考文獻(xiàn)［12］。

2.3 數(shù)值求解

將僅與網(wǎng)格幾何形狀有關(guān)的型函數(shù)積分及積分線長(zhǎng)度作為幾何常數(shù)處理，計(jì)算一次并存儲(chǔ)，減少計(jì)算量。離散的平衡方程（17）中涉及2個(gè)方向的位移，采用耦合解法同時(shí)求解，避免迭代，從而一次計(jì)算獲得整個(gè)熱應(yīng)力場(chǎng)。采用Fortran語(yǔ)言編程實(shí)現(xiàn)數(shù)值模擬。

3 FGM應(yīng)力問(wèn)題

適用于功能梯度TBC熱應(yīng)力問(wèn)題的數(shù)值方法，必須能夠準(zhǔn)確計(jì)算FGM應(yīng)力場(chǎng)。為了評(píng)價(jià)SCVFVM對(duì)FGM應(yīng)力問(wèn)題的適用性，計(jì)算文獻(xiàn)［15］中的平面應(yīng)變FGM圓盤(pán)應(yīng)力場(chǎng)，幾何模型見(jiàn)圖2。圓盤(pán)內(nèi)徑ri＝0.025 4 m，受到法向壓力pi＝68.95 MPa。圓盤(pán)外徑ro＝0.050 8 m，受到法向壓力po＝6.89 MPa。圓盤(pán)泊松比為常數(shù)μ＝0.22，楊氏模量沿徑向變化E＝E0r2，其中E0＝300 GPa，r為半徑。

圖2 FGM圓盤(pán)示意圖Fig.2 Sketch map of the FGM disk

選取圓盤(pán)的1／4作為計(jì)算域，網(wǎng)格劃分與文獻(xiàn)［15］相同。圖3為不同方法得到的計(jì)算結(jié)果，其中Q4代表4節(jié)點(diǎn)四邊形單元，Q8代表8節(jié)點(diǎn)四邊形單元。通過(guò)比較可以看出SCV-FVM結(jié)果與理論解吻合良好，且能夠有效避免應(yīng)力場(chǎng)的不連續(xù)，而HOTFGM理論及傳統(tǒng)FEM不能。另外，SCV-FVM結(jié)果與理論解誤差在邊界附近相對(duì)較大，這是由插值誤差引起的。SCV-FVM以位移為求解變量，單元中心的應(yīng)力根據(jù)式（6）和（7）計(jì)算，然后插值得到節(jié)點(diǎn)應(yīng)力，從而引入誤差。

圖3 FGM圓盤(pán)計(jì)算結(jié)果Fig.3 Results of the FGM disk

4 功能梯度TBC熱應(yīng)力問(wèn)題

將本文發(fā)展的SCV-FVM用于研究含F(xiàn)GM層的平面應(yīng)變TBC熱力性能，計(jì)算模型見(jiàn)圖4，材料屬性見(jiàn)表1。其中第3層（L3）為FGM層，其物性參數(shù)變化規(guī)律為

式中：P代表任意物性參數(shù)，Pa、Pb分別為陶瓷頂層L4和粘結(jié)層L2的物性參數(shù)。

表1 功能梯度TBC材料屬性Table 1 Material properties of the functionally graded TBC

圖4 功能梯度TBC計(jì)算模型示意圖Fig.4 Sketch map of the functionally graded TBC

TBC下表面溫度Td＝25℃，上表面溫度沿x方向變化Tu＝1 325cos20（|x-1|）+25℃，左右表面絕熱。參考溫度Tr＝0℃。TBC下表面簡(jiǎn)支，其余自由。采用均勻四邊形單元?jiǎng)澐钟?jì)算域，網(wǎng)格總數(shù)為5 000。

當(dāng)金屬底層L1厚度為Δy1＝0.2 m，L2厚度為Δy2＝0.1 m，L3厚度為Δy3＝0.5 m，L4厚度為Δy4＝0.2 m時(shí)，熱應(yīng)力場(chǎng)云圖如圖5所示。由于L1和L2交界面處物性參數(shù)突變，存在應(yīng)力集中現(xiàn)象（圖5（c））。FGM層L3起到過(guò)渡作用，L2、L3及L3、L4交界面處沒(méi)有應(yīng)力集中。SCV-FVM計(jì)算結(jié)果沒(méi)有出現(xiàn)應(yīng)力不連續(xù)現(xiàn)象。比較SCV-FVM與ANSYS在y＝0.5 m處的計(jì)算結(jié)果（見(jiàn)圖6），兩者吻合良好，驗(yàn)證本文方法的正確性。

圖5 TBC熱應(yīng)力場(chǎng)（Δy1＝0.2 m，Δy2＝0.1 m，Δy3＝0.5 m，Δy4＝0.2 m）Fig.5 The thermoelastic field of the TBC（Δy1＝0.2 m，Δy2＝0.1 m，Δy3＝0.5 m，Δy4＝0.2 m）

保持Δy2＝0.1 m，Δy4＝0.2 m，Δy1+Δy3＝0.7 m，改變?chǔ)3研究FGM層厚度對(duì)TBC熱應(yīng)力場(chǎng)的影響，見(jiàn)圖6。隨著Δy3的增加，位移ux及溫度T增大。合理選擇FGM層厚度，能夠有效減小應(yīng)力幅值，使應(yīng)力分布更光滑，如Δy3＝0.5 m和Δy3＝0.4 m。不合理的FGM層厚度不但不能改善應(yīng)力分布，甚至使應(yīng)力幅值比無(wú)FGM層時(shí)大得多，如Δy3＝0.3 m和Δy3＝0.2 m。

圖6 TBC中y＝0.5 m處的熱應(yīng)力曲線Fig.6 The curves of variables at y＝0.5 m in the TBC

4 結(jié)束語(yǔ)

本文發(fā)展SCV-FVM，采用交錯(cuò)網(wǎng)格技術(shù)，使得SCV-FVM適用于不均勻材料穩(wěn)態(tài)熱應(yīng)力分析。計(jì)算FGM圓盤(pán)應(yīng)力場(chǎng)，結(jié)果與理論解吻合良好，通過(guò)與文獻(xiàn)對(duì)比發(fā)現(xiàn)SCV-FGM能夠有效避免不合理的應(yīng)力不連續(xù)現(xiàn)象?；赟CV-FVM進(jìn)一步研究功能梯度TBC熱應(yīng)力問(wèn)題，分析涂層厚度對(duì)TBC熱力性能的影響。計(jì)算得到的熱應(yīng)力場(chǎng)不存人工數(shù)值不連續(xù)現(xiàn)象，結(jié)果與ANSYS吻合良好，表明SCV-FVM可以用作TBC設(shè)計(jì)階段的數(shù)值預(yù)測(cè)工具，作為FEM的有效補(bǔ)充方法。

［1］胡春枝，肖金生.陶瓷／金屬梯度熱障涂層的制備及性能評(píng)價(jià)［J］.海軍工程大學(xué)學(xué)報(bào)，2006，18（6）：75-78.HU Chunzhi，XIAO Jinsheng.Fabrication and characterization of ceramic／metal gradient thermal barrier coating［J］.Journal of Naval University of Engineering，2006，18（6）：75-78.

［2］原紅星.等離子噴涂法制備鋁基體梯度陶瓷涂層的組織及性能研究［D］.西安：西安理工大學(xué)，2009：21-45.YUAN Hongxing.Structure and performance of gradient ceramic coatings on aluminum substrate prepared by plasma spray［D］.Xi’an：Xi’an University of Technology，2009：21-45.

［3］HEJWOWSKI T，WERONSKI A.The effect of thermal barrier coatings on diesel engine performance［J］.Vacuum，2002，65：427-432.

［4］PARLAK A，SAHIN B，YASAR H.Performance optimisation of an irreversible dual cycle with respect to pressure ratio and temperature ratio-experimental results of a ceramic coated IDI diesel engine［J］.Energy Conversion and Management，2004，45：1219-1232.

［5］PARLAK A，YASAR H，ELDOGAN O.The effect of thermal barrier coating on a turbo-charged diesel engine performance and exergy potential of the exhaust gas［J］.Energy Conversion and Management，2005，46：489-499.

［6］劉杰，肖金生，覃峰，等.陶瓷／金屬梯度熱障涂層圓筒的傳熱與熱應(yīng)力有限差分分析［J］.武漢理工大學(xué)學(xué)報(bào)：交通科學(xué)與工程版，2002，26（3）：379-382.LIU Jie，XIAO Jinsheng，QIN Feng，et al.Finite difference analysis for heat transfer and thermal stress of ceramic／metal gradient thermal barrier coating cylinder［J］.Journal of Wuhan University of Technology：Transportation Science and Engineering，2002，26（3）：379-382.

［7］劉春曉.催化／熱障功能梯度涂層的建模、仿真與實(shí)驗(yàn)研究［D］.武漢：武漢理工大學(xué)，2003：30-47.LIU Chunxiao.Modeling simulation and experimental study of catalysis／thermal barrier coating functionally gradient coatings［D］.Wuhan：Wuhan University of Technology，2003：30-47.

［8］王素，倪春陽(yáng)，朱心雄.功能梯度材料活塞三維溫度場(chǎng)分析［J］.北京航空航天大學(xué)學(xué)報(bào)，2005，31（12）：1299-1302.WANG Su，NI Chunyang，ZHU Xinxiong.Temperature field analysis for heterogeneous material piston［J］.Journal of Beijing University of Aeronautics and Astronautics，2005，31（12）：1299-1302.

［9］杜雙松.梯度熱障涂層活塞的穩(wěn)態(tài)熱分析研究［D］.合肥：合肥工業(yè)大學(xué)，2008：29-54.DU Shuangsong.Steady-state thermal analysis of functional gradient materials as thermal barrier coating of engine piston［D］.Hefei：Hefei University of Technology，2008：29-54.

［10］BANSAL Y，PINDERA M J.Efficient reformulation of the thermoelastic higher-order theory for functionally graded materials［J］.Journal of Thermal Stresses，2003，26：1055-1092.

［11］BUYUKKAYA E.Thermal analysis of functionally graded coating AlSi alloy and steel pistons［J］.Surface and Coatings Technology，2008，202：3856-3865.

［12］龔京風(fēng)，明平劍，宣領(lǐng)寬，等.基于有限體積法求解FGM動(dòng)態(tài)響應(yīng)及固有特性［J］.華中科技大型學(xué)報(bào)：自然科學(xué)版，2013，41（5）：25-30.GONG Jingfeng，MING Pingjian，XUAN Lingkuan，et al.Solving dynamic performance and natural characteristics of FGM by finite volume method［J］.Journal of Huazhong U-niversity of Science＆Technology：Natural Science Edition，2013，41（5）：25-30.

［13］KIM J H，PAULINO G H.Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials［J］.ASME Journal of Applied Mechanics，2002，69：502-514.

［14］CHAROENSUK J，VESSAKOSOL P.Numerical solutions for functionally graded solids under thermal and mechanical loads using a high-order control volume finite element method［J］.Applied Thermal Engineering，2011，31：213-227.

［15］CAVALCANTE M A A，MARQUES S P C，PINDERA M J.Computational aspects of the parametric finite-volume theory for functionally graded materials［J］.Computational Materials Science，2008，44：422-438.

［16］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：819-822.

［17］ABOUDI J，PINDERA M J，ARNOLD S M.Higher-order theory for functionally graded materials［J］.Composites Part B，1999，30：777-832.

［18］ZHONG Y，BANSAL Y，PINDERA M J.Efficient reformulation of the thermal higher-order theory for FGMs with locally variable thermal conductivity［J］.International Journal of Computational Engineering Science，2004，5：795-831.

［19］CAVALCANTE M A A，MARQUES S P C，PINDERA M J.Parametric formulation of the finite-volume theory for functionally graded materials—Part I：analysis［J］.ASME Journal of Applied Mechanics，2007，74：935-945.

［20］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］.ASME Journal of Applied Mechanics，2007，74：946-957.

Research on the finite volume method for the thermoelastic stress analysis of thermal barrier coatings

GONG Jingfeng1，2，ZHANG Wenping2，XUAN Lingkuan2，MING Pingjian2，HAN Chunxun2
（1.School of Energy and Power Engineering，Wuhan University of Technology，Wuhan 430063，China；2.College of Power and Energy Engineering，Harbin Engineering University，Harbin 150001，China）

In order to analyze the thermoelastic performance of a thermal barrier coating with a functionally graded material（FGM）layer，the staggered grid finite volume method（SCV-FVM）has been developed.The staggered grid technique was adopted by defining variables at the vertex and material properties at the cell center.SCV-FVM was applied to predict the stress of an FGM disk.The results agree well with the theoretical solution which validates its accuracy.The comparisons between the results from SCV-FVM，the conventional finite element method and the higher order theory functionally graded method（HOTFGM）show that SCV-FVM can avoid the stress discontinuity caused by the change in the matter parameters.The thermoelastic performance of functionally graded TBC was investigated.The obtained results are consistent with the ANSYS results which validates the simulations.The effects of the thickness have been discussed.

thermal barrier coating；functionally graded；finite volume method；staggered grid；stress field

10.3969／j.issn.1006-7043.201306012

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006-7043.201306012.html

O343.1

A

1006-7043（2014）06-0713-06

2013-06-03.網(wǎng)絡(luò)出版時(shí)間：2014-05-14 15：49：09.

中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金資助項(xiàng)目（HEUCF130302）.

龔京風(fēng)（1986-），女，博士研究生；張文平（1956-），男，教授，博士生導(dǎo)師.

張文平，E-mail：zhangwenping＠hrbeu.edu.cn.