李顶河 麻硕

摘要：基于热力耦合扩展逐层/实体元方法（TM-XLWM/SE）对热障涂层结构传热过程中的力学行为进行了研究。首先，利用热力耦合扩展逐层方法和三维实体元方法分别建立了热障涂层和金属基体的控制方程，然后，根据热障涂层和基体连接区域的位移、温度协调和内力平衡条件，将热障涂层和基体控制方程耦合为总体控制方程。依据总体控制方程，编写了相应的C++程序，通过与有限元软件Comsol建立的三维弹性模型进行对比，验证了方法的正确性。最后，对含分层损伤和界面脱黏损伤的热障涂层结构进行了热力耦合研究，得到了损伤对热力响应的影响规律。同时，所研究的外形更接近真实的发动机涡轮叶片，具有一定的实际应用价值。

关键词：热障涂层结构；动态热力耦合分析；分层损伤；界面脱黏；扩展逐层

中图分类号：V232.4文献标识码：ADOI：10.19452/j.issn1007-5453.2021.09.002

基金项目：天津市教委科研计划项目（2018KJ241）

涡轮叶片作为发动机的核心部件，其承温和承载条件极为苛刻。随着燃气涡轮发动机的不断发展，燃气进口温度越来越高，仅仅依靠空气冷却已经不能满足使用要求。为了提高耐高温、耐腐蚀能力，热障涂层技术被越来越多地应用于涡轮叶片结构。

热障涂层结构是极其复杂的系统工程[1]，结构、环境[2]、失效模式复杂，在服役过程中难以避免地会有各种损伤存在。然而目前，对于热障涂层结构损伤机理的研究较少，已有的研究主要是针对简单的平板模型、半圆形模型等，而非复杂的曲面结构。近年来，有越来越多的学者关注了这一问题，段力等[3]采用有限元固体热传导仿真模型对热障涂层的温度特性进行了分析。周益春等[4]通过试验模拟了热障涂层的服役状态，分别从试验和理论，研究了热障涂层热力耦合问题的破坏机制。陈琛等[5]对横向梯度温度载荷作用下的热障涂层结构进行了研究，并结合有限元方法分析了失效机理。徐惠彬等[6]对热障涂层在高温蠕变以及高温低周疲劳作用下的损伤过程进行了研究。董丽[7]利用有限元软件对热循环作用下的残余应力进行了数值模拟研究。张文东等[8]基于数字图像相关法研究了裂纹识别技术，能够准确获得裂纹扩展试验的裂纹信息。

发动机热障涂层结构在非常复杂的力-热耦合环境下工作，热力学动态分析是该结构设计、优化和制造的关键和瓶颈问题，同时各种损伤形式所造成的影响也必须被考虑。而扩展逐层方法（XLWM）非常适用于对这类问题进行研究。关于扩展逐层方法，Li等[9-13]已经做了大量的工作，XLWM能够同时考虑多种损伤形式，计算量小，精度高。同时，其有限元控制方程中包含了上下表面的物理量自由度，这使得该方法能很方便地通过协调条件和内力平衡条件与其他方法相耦合，建立复杂结构的分析方法。同时，XLWM在多物理场耦合方面也有着大量的研究成果。Li和Fish[14]将扩展逐层方法应用于热力耦合问题的研究，对具有多个分层和横向裂纹的复合材料层合板进行热力耦合动态分析。Li和Shan等[15]进行了复合材料层合板的稳态热力耦合分析。

本文将应用考虑热力耦合的扩展逐层和三维实体元方法进行研究。热力耦合扩展逐层方法（TM-XLWM）便于对含损伤结构进行分析，而热力耦合三维实体元方法易于对复杂结构进行建模。结合两种方法的分析优势，建立了TM-XLWM/SE方法。基于TM-XLWM/SE法，对热障涂层结构传热过程中的力学行为进行研究，同时研究了界面脱黏损伤和分层裂纹对热力响应的影响。

1热力耦合问题的数值求解方法

1.1热力耦合混合变分原理

2热障涂层结构热力耦合分析模型的建立

2.1局部坐标系和整体坐标系的转换

在第1节中介绍了扩展逐层单元在局部坐标系下的刚度矩阵。为了对结构复杂的热障涂层结构进行研究，需要确定总体坐标系。如图2所示，局部坐标系用x， y， z表示，其中x， y轴在单元的中面内，z轴垂直于中面。总体坐标系用x′， y′， z′表示。

单元坐标系转换的整体思路是：将局部坐标系下的单元刚度矩阵和载荷矢量转换到总体坐标系下，获得整体结构的总体刚度矩阵，进行控制方程求解，可得总体坐标系下的位移矢量a′，将其转换到局部坐标系下，便可计算单元的应力。

2.2热障涂层结构简化模型

目前所分析的热障涂层结构模型大多是理想化的板结构、圆形结构[17-19]。但是，实际应用的涡轮叶片模型是非常复杂的曲面结构，如图3所示。叶片不同位置的曲率、厚度不同，同时热障涂层为多层结构，各层厚度差异很大，从而导致有限元分析比较困难。本节将对热障涂层结构进行建模，为有限元模拟奠定基础。

本节对叶片的几何模型进行了一定程度的簡化[20]，并参考文献[20]～文献[26]进行了热障涂层结构建模。忽略了冷却通道的具体形状与数目，假设成只有单个冷却通道。即使做此简化，也不会对温度场、应力场等人们所关注的结果产生太大的误差。模型中陶瓷层用TBC表示，厚度为ht；涡轮叶片基底用SUB表示，厚度为hs。简化后的模型和尺寸如图4所示。

2.3功能梯度材料参数模型

热障涂层结构多为双层结构[27]，这种结构的局限性是会引起较大的残余应力，服役过程中也会有较大的热应力产生，这会大大缩短结构的使用年限。那么如果能连续地或近连续地改变成分比例和组织结构，使之形成成分梯度、结构梯度或功能梯度，便可以消除宏观界面，缓解残余应力和热应力的影响。由此产生了功能梯度热障涂层结构[28]，如图5所示。

2.5程序实现

为了便于理解所建立的分析模型，本节给出了相应的程序流程图，如图6所示。首先，分别读取依据热力耦合扩展逐层方法和三维实体元方法建模的热障涂层和金属基体输入文件。随后处理它们的边界条件，并确定两种方法的稀疏矩阵维度。依据式（8）和式（11）分别进行两结构刚度矩阵、质量矩阵和阻尼矩阵的拼装，并进行坐标系的转换。根据式（14）和式（15），分别确定两结构的接触和非接触自由度，并进行行列变换。根据式（17），拼装整体结构的TMXLWM/SE方法控制方程。根据式（37），计算K矩阵和Ft +？t，求解矢量Qt +？t。最后，进行时间积分。当满足所设定的步数要求，则停止运算。

3数值算例

3.1收敛性分析及正确性验证

對于热障涂层结构，陶瓷层（TBC）采用氧化锆[30]，基体（SUB）采用GH4169镍基合金[31]。TBC的热力学参数为：弹性模量E = 34GPa，泊松比γ=0.12，密度ρ= 5600kg/m3，热膨胀系数α= 8.657×10-6K-1，导热系数κ= 2.09W/（m？K），比热容c = 460J/（kg？K）。SUB的热力学参数为：弹性模量E = 160GPa，泊松比γ=0.35，密度ρ= 8240kg/m3，热膨胀系数α= 17×10-6K-1，导热系数κ= 23.6W/（m？K），比热容c = 615J/（kg？K）。在本算例中，每层均为各向同性，边界条件为底面固支。整个结构初始温度为0，当最外层温度突然增加为800K时，研究其传热过程中产生的一系列热力响应。

首先，进行了网格的收敛性分析，选择时间步？t = 0.005s，以及4种不同的网格划分方式，如图7所示。研究任一点A（20mm，28.38mm，4.5mm）物理量随时间的变化情况，如图8所示。从图8中可以看出，随着网格的不断加密，曲线逐渐重合，结果不断收敛，第三、四划分方式的物理量变化极小，则第三种网格划分方式已经收敛。取第三种网格划分方式，选择4种不同的时间步，物理量变化曲线如图9所示，在时间步取？t = 0.005s时，计算结果已经收敛。在之后的数值算例中，取第三种网格划分方式，时间步？t = 0.005s进行计算。

依据所选取的网格划分方式和时间步，在有限元软件Comsol中进行了相同工况的算例计算，网格划分方式和时间步相一致，用以验证方法的正确性，对比情况在图10中给出。从图中可以看出，变化曲线基本重合，证明了应用TM-XLWM/SE方法所建立的模型是正确的。

3.2功能梯度热障涂层结构算例计算

为缓解双层结构残余应力的影响，热障涂层常选用功能梯度材料。本算例中，热障涂层顶层选用ZrO2隔热材料，ZrO2的热力学参数为：弹性模量E = 34GPa，泊松比？= 0.12，密度ρ= 5600kg/m3，热膨胀系数α= 8.657×10-6K-1，导热系数κ= 2.09W/（m？K），比热容c = 460J/（kg？K）。底层选择NiCoCrAlY材料，NiCoCrAlY的热力学参数为：弹性模量E = 214.5GPa，泊松比γ=0.3，密度ρ= 7320kg/m3，热膨胀系数α= 11.6×10-6K-1，导热系数κ= 16.1W/（m？K），比热容c = 501J/（kg？K）。选取了功能梯度指数具有代表性的三个值n（0，0.5，1.0）进行了算例计算，工况和上节相同。A点位移、温度随时间的变化曲线如图11所示，由于功能梯度指数的不断变化，材料在复合体中的成分比例不断变化，对物理量的变化速度、达到稳定的时间会产生较大影响。从图中可以看出，随着功能梯度指数的增大，位移、温度的变化速度会逐渐加快，会更快趋于稳定。

3.3含损伤功能梯度热障涂层结构算例计算

热障涂层结构在服役过程中长期受到热循环的作用，会导致各种损伤的出现，进而导致热障涂层体系的失效。在本算例中，模拟热障涂层在厚度方向上层间裂纹对位移、温度分布的影响。整个结构的工况和之前相一致，功能梯度指数取0.5。分两种情况进行研究：（1）结构不含任何损伤；（2）在涂层厚度方向的中间位置处有分层损伤，如图12所示。热障涂层结构上一点B（20mm，14.2mm，2.5mm）位移、温度随时间的变化如图13所示，B点在损伤区域的投影为区域的中间位置。从图12中可以看出，层间裂纹的存在影响了温度的正常传导，温度达到稳定的时间会被延长。同时，裂纹的存在对B点的位移变化产生显著影响，位移变化小于不含裂纹的情况。

涂层的部分剥落，即热障涂层和金属基体间脱黏区域的出现也是热障涂层体系失效的一个重要原因。关于这一问题，分两种情况进行研究：（1）结构不含损伤；（2）存在含脱黏区域，脱黏区域的大小和上一算例相同，位置在交界面处，如图14所示。B点（20mm，14.2mm，2.5mm）位移、温度随时间的变化情况如图15所示。由于损伤区域的大小相同，位置相近，各物理量随时间的变化规律是相类似的。

最后，对同时含有层间裂纹、界面脱黏的热障涂层结构进行了分析。在本算例中考虑了三种不同的情况：（1）结构不含任何损伤；（2）结构只含分层损伤，损伤区域和之前相同；（3）结构同时含有两种损伤，损伤区域的大小、位置不变。B点（20mm，14.2mm，2.5mm）位移、温度随时间的变化情况如图16所示。将情况1和情况3进行比较可以看出，裂纹的存在阻隔了温度的正常传导，增加了B点位移、温度达到稳定的时间，位移变化的速率和幅值会小于无损情况。比较情况2和情况3可以看出，u1、u2、T随时间的变化情况基本一致。而u3随时间的变化出现了差异。对于B点而言，u3为与表面相垂直方向的位移。在温度传导的起始阶段，分层损伤的影响占主导，u3随时间的变化曲线是重合的。随着温度传导的进行，由于界面脱黏的存在，二者出现了差异。当达到稳态时，情况3的u3值要略小于情況2。

4结论

本文基于热力耦合扩展逐层方法和三维实体元方法，建立了热障涂层结构热力耦合分析模型。通过数值算例验证了方法的正确性，同时分析了不同功能梯度指数对物理量分布的影响，也对含层间裂纹、界面脱黏，以及同时含有两种损伤的模型进行了分析。得到了相应规律：（1）对于功能梯度热障涂层结构，随着功能梯度指数的增大，位移、温度的变化速度会逐渐加快，会更快趋于稳定。（2）对于分别含分层、脱黏损伤的热障涂层结构，由于裂纹的存在，基体上投影于损伤区域内的点位移、温度达到稳定的时间被大大延长，位移的变化会小于不含裂纹的情况。

与之前的研究成果相比，所建立的TM-XLWM/SE方法的优点包括：（1）同时应用两种方法进行建模，不需要增添任何几何假设，便可对复杂结构进行含损伤的热力耦合分析，同时也保障了分析的精度；（2）考虑了几种典型的损伤形式，可对同时具有多种损伤的热障涂层结构进行分析，这是其他方法很难实现的；（3）本文提出的损伤分析模型可用于结构复杂的热障涂层结构的热力分析、设计和优化，而不是简单的板结构或圆形结构。

参考文献

[1]李应红，魏悦广，周益春.航空发动机涡轮叶片热障涂层[J].湘潭大学学报（自然科学版），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]邹学锋，潘凯，燕群，等.多场耦合环境下高超声速飞行器结构动强度问题综述[J].航空科学技术，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]段力，高均超，汪瑞军，等.航空发动机叶片表面热障涂层温度分布的仿真分析[J].上海交通大学学报， 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]周益春，桥田俊之，段祝平.热障涂层的热-力耦合破坏机制[C]//“力学2000”学术大会.北京：中国力学学会， 2000. Zhou Yichun， Hashida Toshiyuki， Duan Zhuping. Thermalmechanical coupling failure mechanism of thermal barrier coatings[C]//“Mechanics2000”Proceedings of Academic Conferences. Beijing：Chinese Society of Mechanics， 2000.（in Chinese）

[5]陈琛，郭洪波，宫声凯.横向梯度温度场下热障涂层的失效分析[J].中国腐蚀与防护学报， 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]徐惠彬，宫声凯，陈立强，等.热、力耦合作用下热障涂层的失效机制[J].北京航空航天大學学报， 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热障涂层残余应力有限元模拟[D].大连：大连理工大学， 2012. Dong Li. Finite-element analysis of the residual stress in ZrO2thermal barrier coatings[D]. Dalian： Dalian University of Technology，2012.（in Chinese）

[8]张文东，卓轶，董登科，等.一种裂纹识别方法的研究及试验验证[J].航空科学技术，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，Fish 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]郑成洲.低速冲击对热障涂层动态响应的影响[J].机械工程与自动化， 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]杨雪松.基于内聚力单元热障涂层体系失效过程有限元模拟分析[D].湘潭：湘潭大学， 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溫度场及应力场模拟分析[D].天津：中国民航大学， 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].湘潭：湘潭大学， 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]毛卫国.热—力联合作用下热障涂层界面破坏分析[D].湘潭：湘潭大学， 2006. Mao Weiguo. Analysis of interface failure of thermal barrier ceramiccoatingunderthermo-mechanicalloadings[D]. Xiangtan： Xiangtan University，2006.（in Chinese）

[22]张治彪.基于真实TGO界面形貌的热障涂层热应力及界面失效有限元分析[D].湘潭：湘潭大学， 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]王明军.有限元方法分析冷却气膜孔对涡轮叶片TBCs温度场和应力场的影响[D].湘潭：湘潭大学， 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]李晓军.考虑冷却通道的热障涂层涡轮叶片应力场的有限元模拟[D].湘潭：湘潭大学， 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]董平.航空发动机气冷涡轮叶片的气热耦合数值模拟研究[D].哈尔滨：哈尔滨工业大学， 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]廖涛，董一巍，张赛涛，等.基于虚拟测量的涡轮叶片气膜孔误差分析方法[J].航空科学技术，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]郭洪波，宫声凯，徐惠彬.梯度热障涂层的设计[J].航空学报，

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中国功能材料科技与产业高层论坛， 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].天津：中国民航大学，2019. Yang Xiao. Extended layerwise method for dynamic fracture of functionally graded materials[D]. Tianjin： Civil Aviation University of China， 2019.（in Chinese）

[30]单武奎.多层结构多物理场断裂问题扩展逐层法研究[D].天津：中国民航大学，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镍基高温合金脉冲激光焊工艺参数优化[J].工具技术，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