含損傷熱障涂層結(jié)構(gòu)熱力耦合問(wèn)題的擴(kuò)展逐層/實(shí)體元方法研究

李頂河 麻碩

摘要：基于熱力耦合擴(kuò)展逐層/實(shí)體元方法（TM-XLWM/SE）對(duì)熱障涂層結(jié)構(gòu)傳熱過(guò)程中的力學(xué)行為進(jìn)行了研究。首先，利用熱力耦合擴(kuò)展逐層方法和三維實(shí)體元方法分別建立了熱障涂層和金屬基體的控制方程，然后，根據(jù)熱障涂層和基體連接區(qū)域的位移、溫度協(xié)調(diào)和內(nèi)力平衡條件，將熱障涂層和基體控制方程耦合為總體控制方程。依據(jù)總體控制方程，編寫(xiě)了相應(yīng)的C++程序，通過(guò)與有限元軟件Comsol建立的三維彈性模型進(jìn)行對(duì)比，驗(yàn)證了方法的正確性。最后，對(duì)含分層損傷和界面脫黏損傷的熱障涂層結(jié)構(gòu)進(jìn)行了熱力耦合研究，得到了損傷對(duì)熱力響應(yīng)的影響規(guī)律。同時(shí)，所研究的外形更接近真實(shí)的發(fā)動(dòng)機(jī)渦輪葉片，具有一定的實(shí)際應(yīng)用價(jià)值。

關(guān)鍵詞：熱障涂層結(jié)構(gòu)；動(dòng)態(tài)熱力耦合分析；分層損傷；界面脫黏；擴(kuò)展逐層

中圖分類號(hào)：V232.4文獻(xiàn)標(biāo)識(shí)碼：ADOI：10.19452/j.issn1007-5453.2021.09.002

基金項(xiàng)目：天津市教委科研計(jì)劃項(xiàng)目（2018KJ241）

渦輪葉片作為發(fā)動(dòng)機(jī)的核心部件，其承溫和承載條件極為苛刻。隨著燃?xì)鉁u輪發(fā)動(dòng)機(jī)的不斷發(fā)展，燃?xì)膺M(jìn)口溫度越來(lái)越高，僅僅依靠空氣冷卻已經(jīng)不能滿足使用要求。為了提高耐高溫、耐腐蝕能力，熱障涂層技術(shù)被越來(lái)越多地應(yīng)用于渦輪葉片結(jié)構(gòu)。

熱障涂層結(jié)構(gòu)是極其復(fù)雜的系統(tǒng)工程[1]，結(jié)構(gòu)、環(huán)境[2]、失效模式復(fù)雜，在服役過(guò)程中難以避免地會(huì)有各種損傷存在。然而目前，對(duì)于熱障涂層結(jié)構(gòu)損傷機(jī)理的研究較少，已有的研究主要是針對(duì)簡(jiǎn)單的平板模型、半圓形模型等，而非復(fù)雜的曲面結(jié)構(gòu)。近年來(lái)，有越來(lái)越多的學(xué)者關(guān)注了這一問(wèn)題，段力等[3]采用有限元固體熱傳導(dǎo)仿真模型對(duì)熱障涂層的溫度特性進(jìn)行了分析。周益春等[4]通過(guò)試驗(yàn)?zāi)M了熱障涂層的服役狀態(tài)，分別從試驗(yàn)和理論，研究了熱障涂層熱力耦合問(wèn)題的破壞機(jī)制。陳琛等[5]對(duì)橫向梯度溫度載荷作用下的熱障涂層結(jié)構(gòu)進(jìn)行了研究，并結(jié)合有限元方法分析了失效機(jī)理。徐惠彬等[6]對(duì)熱障涂層在高溫蠕變以及高溫低周疲勞作用下的損傷過(guò)程進(jìn)行了研究。董麗[7]利用有限元軟件對(duì)熱循環(huán)作用下的殘余應(yīng)力進(jìn)行了數(shù)值模擬研究。張文東等[8]基于數(shù)字圖像相關(guān)法研究了裂紋識(shí)別技術(shù)，能夠準(zhǔn)確獲得裂紋擴(kuò)展試驗(yàn)的裂紋信息。

發(fā)動(dòng)機(jī)熱障涂層結(jié)構(gòu)在非常復(fù)雜的力-熱耦合環(huán)境下工作，熱力學(xué)動(dòng)態(tài)分析是該結(jié)構(gòu)設(shè)計(jì)、優(yōu)化和制造的關(guān)鍵和瓶頸問(wèn)題，同時(shí)各種損傷形式所造成的影響也必須被考慮。而擴(kuò)展逐層方法（XLWM）非常適用于對(duì)這類問(wèn)題進(jìn)行研究。關(guān)于擴(kuò)展逐層方法，Li等[9-13]已經(jīng)做了大量的工作，XLWM能夠同時(shí)考慮多種損傷形式，計(jì)算量小，精度高。同時(shí)，其有限元控制方程中包含了上下表面的物理量自由度，這使得該方法能很方便地通過(guò)協(xié)調(diào)條件和內(nèi)力平衡條件與其他方法相耦合，建立復(fù)雜結(jié)構(gòu)的分析方法。同時(shí)，XLWM在多物理場(chǎng)耦合方面也有著大量的研究成果。Li和Fish[14]將擴(kuò)展逐層方法應(yīng)用于熱力耦合問(wèn)題的研究，對(duì)具有多個(gè)分層和橫向裂紋的復(fù)合材料層合板進(jìn)行熱力耦合動(dòng)態(tài)分析。Li和Shan等[15]進(jìn)行了復(fù)合材料層合板的穩(wěn)態(tài)熱力耦合分析。

本文將應(yīng)用考慮熱力耦合的擴(kuò)展逐層和三維實(shí)體元方法進(jìn)行研究。熱力耦合擴(kuò)展逐層方法（TM-XLWM）便于對(duì)含損傷結(jié)構(gòu)進(jìn)行分析，而熱力耦合三維實(shí)體元方法易于對(duì)復(fù)雜結(jié)構(gòu)進(jìn)行建模。結(jié)合兩種方法的分析優(yōu)勢(shì)，建立了TM-XLWM/SE方法?；赥M-XLWM/SE法，對(duì)熱障涂層結(jié)構(gòu)傳熱過(guò)程中的力學(xué)行為進(jìn)行研究，同時(shí)研究了界面脫黏損傷和分層裂紋對(duì)熱力響應(yīng)的影響。

1熱力耦合問(wèn)題的數(shù)值求解方法

1.1熱力耦合混合變分原理

2熱障涂層結(jié)構(gòu)熱力耦合分析模型的建立

2.1局部坐標(biāo)系和整體坐標(biāo)系的轉(zhuǎn)換

在第1節(jié)中介紹了擴(kuò)展逐層單元在局部坐標(biāo)系下的剛度矩陣。為了對(duì)結(jié)構(gòu)復(fù)雜的熱障涂層結(jié)構(gòu)進(jìn)行研究，需要確定總體坐標(biāo)系。如圖2所示，局部坐標(biāo)系用x， y， z表示，其中x， y軸在單元的中面內(nèi)，z軸垂直于中面?？傮w坐標(biāo)系用x′， y′， z′表示。

單元坐標(biāo)系轉(zhuǎn)換的整體思路是：將局部坐標(biāo)系下的單元?jiǎng)偠染仃嚭洼d荷矢量轉(zhuǎn)換到總體坐標(biāo)系下，獲得整體結(jié)構(gòu)的總體剛度矩陣，進(jìn)行控制方程求解，可得總體坐標(biāo)系下的位移矢量a′，將其轉(zhuǎn)換到局部坐標(biāo)系下，便可計(jì)算單元的應(yīng)力。

2.2熱障涂層結(jié)構(gòu)簡(jiǎn)化模型

目前所分析的熱障涂層結(jié)構(gòu)模型大多是理想化的板結(jié)構(gòu)、圓形結(jié)構(gòu)[17-19]。但是，實(shí)際應(yīng)用的渦輪葉片模型是非常復(fù)雜的曲面結(jié)構(gòu)，如圖3所示。葉片不同位置的曲率、厚度不同，同時(shí)熱障涂層為多層結(jié)構(gòu)，各層厚度差異很大，從而導(dǎo)致有限元分析比較困難。本節(jié)將對(duì)熱障涂層結(jié)構(gòu)進(jìn)行建模，為有限元模擬奠定基礎(chǔ)。

本節(jié)對(duì)葉片的幾何模型進(jìn)行了一定程度的簡(jiǎn)化[20]，并參考文獻(xiàn)[20]～文獻(xiàn)[26]進(jìn)行了熱障涂層結(jié)構(gòu)建模。忽略了冷卻通道的具體形狀與數(shù)目，假設(shè)成只有單個(gè)冷卻通道。即使做此簡(jiǎn)化，也不會(huì)對(duì)溫度場(chǎng)、應(yīng)力場(chǎng)等人們所關(guān)注的結(jié)果產(chǎn)生太大的誤差。模型中陶瓷層用TBC表示，厚度為ht；渦輪葉片基底用SUB表示，厚度為hs。簡(jiǎn)化后的模型和尺寸如圖4所示。

2.3功能梯度材料參數(shù)模型

熱障涂層結(jié)構(gòu)多為雙層結(jié)構(gòu)[27]，這種結(jié)構(gòu)的局限性是會(huì)引起較大的殘余應(yīng)力，服役過(guò)程中也會(huì)有較大的熱應(yīng)力產(chǎn)生，這會(huì)大大縮短結(jié)構(gòu)的使用年限。那么如果能連續(xù)地或近連續(xù)地改變成分比例和組織結(jié)構(gòu)，使之形成成分梯度、結(jié)構(gòu)梯度或功能梯度，便可以消除宏觀界面，緩解殘余應(yīng)力和熱應(yīng)力的影響。由此產(chǎn)生了功能梯度熱障涂層結(jié)構(gòu)[28]，如圖5所示。

2.5程序?qū)崿F(xiàn)

為了便于理解所建立的分析模型，本節(jié)給出了相應(yīng)的程序流程圖，如圖6所示。首先，分別讀取依據(jù)熱力耦合擴(kuò)展逐層方法和三維實(shí)體元方法建模的熱障涂層和金屬基體輸入文件。隨后處理它們的邊界條件，并確定兩種方法的稀疏矩陣維度。依據(jù)式（8）和式（11）分別進(jìn)行兩結(jié)構(gòu)剛度矩陣、質(zhì)量矩陣和阻尼矩陣的拼裝，并進(jìn)行坐標(biāo)系的轉(zhuǎn)換。根據(jù)式（14）和式（15），分別確定兩結(jié)構(gòu)的接觸和非接觸自由度，并進(jìn)行行列變換。根據(jù)式（17），拼裝整體結(jié)構(gòu)的TMXLWM/SE方法控制方程。根據(jù)式（37），計(jì)算K矩陣和Ft +？t，求解矢量Qt +？t。最后，進(jìn)行時(shí)間積分。當(dāng)滿足所設(shè)定的步數(shù)要求，則停止運(yùn)算。

3數(shù)值算例

3.1收斂性分析及正確性驗(yàn)證

對(duì)于熱障涂層結(jié)構(gòu)，陶瓷層（TBC）采用氧化鋯[30]，基體（SUB）采用GH4169鎳基合金[31]。TBC的熱力學(xué)參數(shù)為：彈性模量E = 34GPa，泊松比γ=0.12，密度ρ= 5600kg/m3，熱膨脹系數(shù)α= 8.657×10-6K-1，導(dǎo)熱系數(shù)κ= 2.09W/（m？K），比熱容c = 460J/（kg？K）。SUB的熱力學(xué)參數(shù)為：彈性模量E = 160GPa，泊松比γ=0.35，密度ρ= 8240kg/m3，熱膨脹系數(shù)α= 17×10-6K-1，導(dǎo)熱系數(shù)κ= 23.6W/（m？K），比熱容c = 615J/（kg？K）。在本算例中，每層均為各向同性，邊界條件為底面固支。整個(gè)結(jié)構(gòu)初始溫度為0，當(dāng)最外層溫度突然增加為800K時(shí)，研究其傳熱過(guò)程中產(chǎn)生的一系列熱力響應(yīng)。

首先，進(jìn)行了網(wǎng)格的收斂性分析，選擇時(shí)間步？t = 0.005s，以及4種不同的網(wǎng)格劃分方式，如圖7所示。研究任一點(diǎn)A（20mm，28.38mm，4.5mm）物理量隨時(shí)間的變化情況，如圖8所示。從圖8中可以看出，隨著網(wǎng)格的不斷加密，曲線逐漸重合，結(jié)果不斷收斂，第三、四劃分方式的物理量變化極小，則第三種網(wǎng)格劃分方式已經(jīng)收斂。取第三種網(wǎng)格劃分方式，選擇4種不同的時(shí)間步，物理量變化曲線如圖9所示，在時(shí)間步??？t = 0.005s時(shí)，計(jì)算結(jié)果已經(jīng)收斂。在之后的數(shù)值算例中，取第三種網(wǎng)格劃分方式，時(shí)間步？t = 0.005s進(jìn)行計(jì)算。

依據(jù)所選取的網(wǎng)格劃分方式和時(shí)間步，在有限元軟件Comsol中進(jìn)行了相同工況的算例計(jì)算，網(wǎng)格劃分方式和時(shí)間步相一致，用以驗(yàn)證方法的正確性，對(duì)比情況在圖10中給出。從圖中可以看出，變化曲線基本重合，證明了應(yīng)用TM-XLWM/SE方法所建立的模型是正確的。

3.2功能梯度熱障涂層結(jié)構(gòu)算例計(jì)算

為緩解雙層結(jié)構(gòu)殘余應(yīng)力的影響，熱障涂層常選用功能梯度材料。本算例中，熱障涂層頂層選用ZrO2隔熱材料，ZrO2的熱力學(xué)參數(shù)為：彈性模量E = 34GPa，泊松比？= 0.12，密度ρ= 5600kg/m3，熱膨脹系數(shù)α= 8.657×10-6K-1，導(dǎo)熱系數(shù)κ= 2.09W/（m？K），比熱容c = 460J/（kg？K）。底層選擇NiCoCrAlY材料，NiCoCrAlY的熱力學(xué)參數(shù)為：彈性模量E = 214.5GPa，泊松比γ=0.3，密度ρ= 7320kg/m3，熱膨脹系數(shù)α= 11.6×10-6K-1，導(dǎo)熱系數(shù)κ= 16.1W/（m？K），比熱容c = 501J/（kg？K）。選取了功能梯度指數(shù)具有代表性的三個(gè)值n（0，0.5，1.0）進(jìn)行了算例計(jì)算，工況和上節(jié)相同。A點(diǎn)位移、溫度隨時(shí)間的變化曲線如圖11所示，由于功能梯度指數(shù)的不斷變化，材料在復(fù)合體中的成分比例不斷變化，對(duì)物理量的變化速度、達(dá)到穩(wěn)定的時(shí)間會(huì)產(chǎn)生較大影響。從圖中可以看出，隨著功能梯度指數(shù)的增大，位移、溫度的變化速度會(huì)逐漸加快，會(huì)更快趨于穩(wěn)定。

3.3含損傷功能梯度熱障涂層結(jié)構(gòu)算例計(jì)算

熱障涂層結(jié)構(gòu)在服役過(guò)程中長(zhǎng)期受到熱循環(huán)的作用，會(huì)導(dǎo)致各種損傷的出現(xiàn)，進(jìn)而導(dǎo)致熱障涂層體系的失效。在本算例中，模擬熱障涂層在厚度方向上層間裂紋對(duì)位移、溫度分布的影響。整個(gè)結(jié)構(gòu)的工況和之前相一致，功能梯度指數(shù)取0.5。分兩種情況進(jìn)行研究：（1）結(jié)構(gòu)不含任何損傷；（2）在涂層厚度方向的中間位置處有分層損傷，如圖12所示。熱障涂層結(jié)構(gòu)上一點(diǎn)B（20mm，14.2mm，2.5mm）位移、溫度隨時(shí)間的變化如圖13所示，B點(diǎn)在損傷區(qū)域的投影為區(qū)域的中間位置。從圖12中可以看出，層間裂紋的存在影響了溫度的正常傳導(dǎo)，溫度達(dá)到穩(wěn)定的時(shí)間會(huì)被延長(zhǎng)。同時(shí)，裂紋的存在對(duì)B點(diǎn)的位移變化產(chǎn)生顯著影響，位移變化小于不含裂紋的情況。

涂層的部分剝落，即熱障涂層和金屬基體間脫黏區(qū)域的出現(xiàn)也是熱障涂層體系失效的一個(gè)重要原因。關(guān)于這一問(wèn)題，分兩種情況進(jìn)行研究：（1）結(jié)構(gòu)不含損傷；（2）存在含脫黏區(qū)域，脫黏區(qū)域的大小和上一算例相同，位置在交界面處，如圖14所示。B點(diǎn)（20mm，14.2mm，2.5mm）位移、溫度隨時(shí)間的變化情況如圖15所示。由于損傷區(qū)域的大小相同，位置相近，各物理量隨時(shí)間的變化規(guī)律是相類似的。

最后，對(duì)同時(shí)含有層間裂紋、界面脫黏的熱障涂層結(jié)構(gòu)進(jìn)行了分析。在本算例中考慮了三種不同的情況：（1）結(jié)構(gòu)不含任何損傷；（2）結(jié)構(gòu)只含分層損傷，損傷區(qū)域和之前相同；（3）結(jié)構(gòu)同時(shí)含有兩種損傷，損傷區(qū)域的大小、位置不變。B點(diǎn)（20mm，14.2mm，2.5mm）位移、溫度隨時(shí)間的變化情況如圖16所示。將情況1和情況3進(jìn)行比較可以看出，裂紋的存在阻隔了溫度的正常傳導(dǎo)，增加了B點(diǎn)位移、溫度達(dá)到穩(wěn)定的時(shí)間，位移變化的速率和幅值會(huì)小于無(wú)損情況。比較情況2和情況3可以看出，u1、u2、T隨時(shí)間的變化情況基本一致。而u3隨時(shí)間的變化出現(xiàn)了差異。對(duì)于B點(diǎn)而言，u3為與表面相垂直方向的位移。在溫度傳導(dǎo)的起始階段，分層損傷的影響占主導(dǎo)，u3隨時(shí)間的變化曲線是重合的。隨著溫度傳導(dǎo)的進(jìn)行，由于界面脫黏的存在，二者出現(xiàn)了差異。當(dāng)達(dá)到穩(wěn)態(tài)時(shí)，情況3的u3值要略小于情況2。

4結(jié)論

本文基于熱力耦合擴(kuò)展逐層方法和三維實(shí)體元方法，建立了熱障涂層結(jié)構(gòu)熱力耦合分析模型。通過(guò)數(shù)值算例驗(yàn)證了方法的正確性，同時(shí)分析了不同功能梯度指數(shù)對(duì)物理量分布的影響，也對(duì)含層間裂紋、界面脫黏，以及同時(shí)含有兩種損傷的模型進(jìn)行了分析。得到了相應(yīng)規(guī)律：（1）對(duì)于功能梯度熱障涂層結(jié)構(gòu)，隨著功能梯度指數(shù)的增大，位移、溫度的變化速度會(huì)逐漸加快，會(huì)更快趨于穩(wěn)定。（2）對(duì)于分別含分層、脫黏損傷的熱障涂層結(jié)構(gòu)，由于裂紋的存在，基體上投影于損傷區(qū)域內(nèi)的點(diǎn)位移、溫度達(dá)到穩(wěn)定的時(shí)間被大大延長(zhǎng)，位移的變化會(huì)小于不含裂紋的情況。

與之前的研究成果相比，所建立的TM-XLWM/SE方法的優(yōu)點(diǎn)包括：（1）同時(shí)應(yīng)用兩種方法進(jìn)行建模，不需要增添任何幾何假設(shè)，便可對(duì)復(fù)雜結(jié)構(gòu)進(jìn)行含損傷的熱力耦合分析，同時(shí)也保障了分析的精度；（2）考慮了幾種典型的損傷形式，可對(duì)同時(shí)具有多種損傷的熱障涂層結(jié)構(gòu)進(jìn)行分析，這是其他方法很難實(shí)現(xiàn)的；（3）本文提出的損傷分析模型可用于結(jié)構(gòu)復(fù)雜的熱障涂層結(jié)構(gòu)的熱力分析、設(shè)計(jì)和優(yōu)化，而不是簡(jiǎn)單的板結(jié)構(gòu)或圓形結(jié)構(gòu)。

參考文獻(xiàn)

[1]李應(yīng)紅，魏悅廣，周益春.航空發(fā)動(dòng)機(jī)渦輪葉片熱障涂層[J].湘潭大學(xué)學(xué)報(bào)（自然科學(xué)版），2020，42（3）：3. Li Yinghong， Wei Yuegaung， Zhou Yichun. Thermal barrier coatings on aero-engine turbine blades[J].Journal of Xiangtan University （Natural Science Edition）， 2020， 42（3）： 3.（in Chinese）

[2]鄒學(xué)鋒，潘凱，燕群，等.多場(chǎng)耦合環(huán)境下高超聲速飛行器結(jié)構(gòu)動(dòng)強(qiáng)度問(wèn)題綜述[J].航空科學(xué)技術(shù)，2020，31（12）：3-15. Zou Xuefeng，Pan Kai，Yan Qun，et al. Overview of dynamic strength of hypersonic vehicle structure in multi-field coupling environment[J].Aeronautical Science & Technology，2020，31（12）：3-15.（in Chinese）

[3]段力，高均超，汪瑞軍，等.航空發(fā)動(dòng)機(jī)葉片表面熱障涂層溫度分布的仿真分析[J].上海交通大學(xué)學(xué)報(bào)， 2017， 51（8）： 915-920. Duan Li， Gao Junchao， Wang Ruijun， et al. Simulation analysis of temperature distribution of turbine blades by thermal buffer coating for aero engine[J]. Journal of Shanghai Jiaotong University， 2017， 51（8）： 915-920.（in Chinese）

[4]周益春，橋田俊之，段祝平.熱障涂層的熱-力耦合破壞機(jī)制[C]//“力學(xué)2000”學(xué)術(shù)大會(huì).北京：中國(guó)力學(xué)學(xué)會(huì)， 2000. Zhou Yichun， Hashida Toshiyuki， Duan Zhuping. Thermalmechanical coupling failure mechanism of thermal barrier coatings[C]//“Mechanics2000”P(pán)roceedings of Academic Conferences. Beijing：Chinese Society of Mechanics， 2000.（in Chinese）

[5]陳琛，郭洪波，宮聲凱.橫向梯度溫度場(chǎng)下熱障涂層的失效分析[J].中國(guó)腐蝕與防護(hù)學(xué)報(bào)， 2013， 33（5）： 400-406. Chen Chen， Guo Hongbo， Gong Shengkai. Failure analysis of thermal barrier coating being subjected to lateral thermal gradient on surface[J]. Journal of Chinese Society for Corrosion and Protection，2013，33（5）： 400-406.（in Chinese）

[6]徐惠彬，宮聲凱，陳立強(qiáng)，等.熱、力耦合作用下熱障涂層的失效機(jī)制[J].北京航空航天大學(xué)學(xué)報(bào)， 2004 （10）： 919-924. Xu Huibin， Gong Shengkai，Chen Liqiang， et al. Failure process of thermal barrier coatings under thermal and mechanical loading[J]. Journal of Beijing University of Aeronautics andAstronautics， 2004 （10）：919-924.（in Chinese）

[7]董麗. ZrO2熱障涂層殘余應(yīng)力有限元模擬[D].大連：大連理工大學(xué)， 2012. Dong Li. Finite-element analysis of the residual stress in ZrO2thermal barrier coatings[D]. Dalian： Dalian University of Technology，2012.（in Chinese）

[8]張文東，卓軼，董登科，等.一種裂紋識(shí)別方法的研究及試驗(yàn)驗(yàn)證[J].航空科學(xué)技術(shù)，2020，31（5）：81-88. Zhang Wendong，Zhuo Yi，Dong Dengke，et al. Research on the method of crack identification and experimental verification[J].Aeronautical Science & Technology，2020，31（5）：81-88.（in Chinese）

[9]Li D H，Liu Y，Zhang X. An extended Layerwise method for composite laminated beams with multiple delaminations and matrix cracks[J]. International Journal for Numerical Methods in Engineering，2015，101（6）：407-434.

[10]Li D H. Delamination and transverse crack growth prediction for laminated composite plates and shells[J]. Computers and Structures，2016，177：39-55.

[11]Li D H，Zhang X，Sze K Y，et al. Extended layerwise method for laminated composite plates with multiple delaminations and transverse cracks[J]. Computational Mechanics，2016，58（4）：657-679.

[12]Li D H，Zhang F，Xu J. Incompatible extended layerwise method for laminated composite shells[J]. International Journal of Mechanical Sciences，2016，119：243-252.

[13]Li D H，Zhang F. Full extended layerwise method for the simulation of laminated composite plates and shells[J]. Computers & Structures，2017，187：101-113.

[14]Li D H，F(xiàn)ish J. Thermomechanical extended layerwise method for laminated composite plates with multiple delaminations and transverse cracks[J]. Composite Structures，2018，185：665-683.

[15]Li D H，Shan W，Zhang F. Steady-state thermomechanical analysis of composite laminated plate with damage based on extended layerwise method[J]. Archive of Applied Mechanics，2020，90（2）：415-435.

[16]Benjeddou A，Andrianarison O. A heat mixed variational theoremforthermoelasticmultilayeredcomposites[J]. Computers & Structures，2006，84（19/20）：1247-1255.

[17]鄭成洲.低速?zèng)_擊對(duì)熱障涂層動(dòng)態(tài)響應(yīng)的影響[J].機(jī)械工程與自動(dòng)化， 2020（5）： 116-118. Zheng Chengzhou. Effect of low-velocity impact on dynamic responseofthermalbarriercoatings[J].Mechanical Engineering &Automation， 2020（5）：116-118.（in Chinese）

[18]楊雪松.基于內(nèi)聚力單元熱障涂層體系失效過(guò)程有限元模擬分析[D].湘潭：湘潭大學(xué)， 2013. Yang Xuesong. Finite element analysis of failure process of thermal barrier coatings by using cohesive element[D]. Xiangtan： Xiangtan University，2013.（in Chinese）

[19]秦令超.激光熔覆CoNiCrAlY溫度場(chǎng)及應(yīng)力場(chǎng)模擬分析[D].天津：中國(guó)民航大學(xué)， 2019. Qin Lingchao. Simulation analysis of temperature field and stress field of laser cladding CoNiCrAlY[D]. Tianjin： Civil Aivation University of China，2019.（in Chinese）

[20]劉奇星.熱障涂層渦輪葉片失效的有限元模擬[D].湘潭：湘潭大學(xué)， 2012. Liu Qixing. Numerical simulation of failure of turbine blade with thermal barrier coatings based on finite element method[D]. Xiangtan： Xiangtan University，2012.（in Chinese）

[21]毛衛(wèi)國(guó).熱—力聯(lián)合作用下熱障涂層界面破壞分析[D].湘潭：湘潭大學(xué)， 2006. Mao Weiguo. Analysis of interface failure of thermal barrier ceramiccoatingunderthermo-mechanicalloadings[D]. Xiangtan： Xiangtan University，2006.（in Chinese）

[22]張治彪.基于真實(shí)TGO界面形貌的熱障涂層熱應(yīng)力及界面失效有限元分析[D].湘潭：湘潭大學(xué)， 2016. Zhang Zhibiao. Analysis of thermal stress and the crack propagation within real TGO interface in thermal barrier coatings by finite element modeling[D]. Xiangtan： Xiangtan University，2016.（in Chinese）

[23]王明軍.有限元方法分析冷卻氣膜孔對(duì)渦輪葉片TBCs溫度場(chǎng)和應(yīng)力場(chǎng)的影響[D].湘潭：湘潭大學(xué)， 2017. Wang Mingjun. Analysis of film cooling holes influence on the temperature and stress field of turbine blade with TBCs by finite element modeling[D]. Xiangtan： Xiangtan University， 2017.（in Chinese）

[24]李曉軍.考慮冷卻通道的熱障涂層渦輪葉片應(yīng)力場(chǎng)的有限元模擬[D].湘潭：湘潭大學(xué)， 2014. Li Xiaojun. Finite element simulation on thermal stress of a turbine blade with thermal barrier coatings considering cooling channels[D].Xiangtan：XiangtanUniversity，2014.（in Chinese）

[25]董平.航空發(fā)動(dòng)機(jī)氣冷渦輪葉片的氣熱耦合數(shù)值模擬研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)， 2009. Dong Ping. Research on conjugate heat transfer simulation of aero turbine engine air-cooled vane [D]. Harbin： Harbin Institute of Technology， 2009.（in Chinese）

[26]廖濤，董一巍，張賽濤，等.基于虛擬測(cè)量的渦輪葉片氣膜孔誤差分析方法[J].航空科學(xué)技術(shù)，2021，32（4）：50-59. Liao Tao，Dong Yiwei，Zhang Saitao，et al. Error analysis method of turbine blade film cooling hole based on virtual measurement[J]. Aeronautical Science & Technology，2021，32（4）：50-59.（in Chinese）

[27]郭洪波，宮聲凱，徐惠彬.梯度熱障涂層的設(shè)計(jì)[J].航空學(xué)報(bào)，

2002（5）： 467-472. Guo Hongbo， Gong Shengkai， Xu Huibin. Design of gradient thermal barrier coatings[J]. Chinese Journal of Aeronautics， 2002（5）： 467-472.（in Chinese）

[28]唐逾，陳德茂，劉慶賓.功能梯度熱障涂層材料[C]//2011中國(guó)功能材料科技與產(chǎn)業(yè)高層論壇， 2011. Tang Yu， Chen Demao， Liu Qingbin. Functional gradients materials for thermal barrier coatings[C]//2011 China Functional Materials Technology and Industry Forum， 2011.（in Chinese）

[29]楊曉.功能梯度材料動(dòng)態(tài)斷裂問(wèn)題的擴(kuò)展逐層法研究[D].天津：中國(guó)民航大學(xué)，2019. Yang Xiao. Extended layerwise method for dynamic fracture of functionally graded materials[D]. Tianjin： Civil Aviation University of China， 2019.（in Chinese）

[30]單武奎.多層結(jié)構(gòu)多物理場(chǎng)斷裂問(wèn)題擴(kuò)展逐層法研究[D].天津：中國(guó)民航大學(xué)，2020. Shan Wukui. Research on extended layerwise method of multiphysics fracture problem for multilayer structures[D]. Tianjin： CivilAviation University of China， 2020.（in Chinese）

[31]吳承隆，尹浩，黃澤涵. GH4169鎳基高溫合金脈沖激光焊工藝參數(shù)優(yōu)化[J].工具技術(shù)，2020，54（10）：38-42 Wu Chenglong，Yin Hao，Huang Zehan. Process parameter optimization of GH4169 pulse laser welding[J].Tool Engineering，2020，54（10）：38-42（in Chinese）

Study on Thermomechanical Extended-layewise/Solid-elements Method for

Thermomechanical Problems of Thermal Barrier Coatings Structure with Damage

Li Dinghe，Ma Shuo

Civil Aviation University of China，Tianjin 300300，China

Abstract： Based on the thermomechanical extended-layerwise/solid-element method （TM-XLWM/SE）， the mechanical behavior of thermal barrier coatings structure is studied during heat transfer. Firstly， the governing equations of the thermal barrier coating and the metal substrate are established by using the thermomechanical extended-layerwise method and the solid-element method， and then the governing equations of the thermal barrier coating and the metal substrate are coupled into the final governing equations according to the conditions of displacement， temperature coordination and internal force balance in the connecting region. According to the final governing equation， the corresponding C++ code is developed. By comparison with the three-dimensional elastic model developed by Comsol， the correctness of TM-XLWM/SE to analyze thermal barrier coating structure is verified. Finally， thermal barrier coatings structures with delamination and interface debonding are studied， and the corresponding influence law is obtained. The structure studied is more similar to the real engine turbine blade， which has certain practical application value.

Key Words：thermal barrier coatings structure；dynamic thermomechanical analysis；delamination damage； interface debonding； XLWM