渦輪葉片的熱流固耦合數值方法研究及應用

摘要：針對航空發動機高壓渦輪葉片缺乏有力的熱流固耦合分析工具的現狀，構建了集流動、換熱與振動于一體的多物理場耦合數值模擬基礎框架與計算方法。依托雙重運動方法和動網格技術，使用跨求解器的弱耦合方法實現了動靜交錯環境下帶陶瓷涂層動葉的熱、流、固耦合數值分析。計算分析了簡化約束條件的渦輪葉片在多物理場環境下的振動特性，獲得了渦輪葉片重要的熱流固特性及其變化規律。結果表明：渦輪葉片在熱流固耦合環境下的振動形態為一階彎曲及扭轉振動形態的混合，其最大位移為0.209mm；動葉的變形會進一步顯著影響作用在葉片上的氣動力和熱負荷，葉尖附近的靜壓變化達9.44%，熱流密度變化達39.7%；陶瓷涂層對熱流固耦合條件下振動特性的影響極小。

關鍵詞：渦輪葉片；熱流固耦合；陶瓷涂層；振動

中圖分類號：V235.1 文獻標志碼：A

DOI：10.7652/xjtuxb202503003 文章編號：0253-987X（2025）03-0021-13

Study on Numerical Method of Thermal-Fluid-Structure Coupling for nbsp;Turbine Blades and Its Application

WANG Xiangyu， GAO Shanghong， FENG Zhenping

（School of Energy and Power Engineering， Xi’an Jiaotong University， Xi’an 710049， China）

Abstract：To address the current lack of effective thermo-fluid-structure coupling analysis tools for high-pressure turbine blades in aircraft engines， a multi-physics field coupling numerical simulation framework and computational method integrating flow， heat transfer， and vibration are developed. Leveraging the dual-motion method and dynamic mesh technology， a weak coupling approach across solvers is employed to achieve thermal-fluid-structure coupling numerical analysis of moving blades with ceramic coatings in a dynamic-steady alternating environment. The vibration characteristics of turbine blades under multi-physics field environments with simplified boundary conditions are computationally analyzed， revealing crucial thermal-fluid-structural properties and their variations. Results indicate that under the thermal-fluid-structure coupling environment， blade vibration exhibits a mix of first-order bending and torsional modes， with a maximum displacement of 0.209mm. Blade deformation further significantly affects aerodynamic forces and thermal loads acting on the blade， with static pressure changes near the blade tip reaching 9.44% and heat flux density variations of 39.7%. The impact of ceramic coatings on vibration characteristics under thermal-fluid-structure coupling conditions is minimal.

Keywords：turbine blade; thermal-fluid-structure coupling; barrier coating; vibration

半個世紀以來，渦輪進口溫度已從1600K左右提升至第四代發動機的2000K，葉片所能承受的熱負荷日益嚴峻。此外，渦輪葉片還承受來自上、下游的非定常氣動載荷^［1-3］，出現共振和高周疲勞，導致渦輪葉片的強度及振動問題日益突出。隨著燃氣溫度上升，葉片材料的導熱系數、彈性模量、屈服極限以及持久強度等均將發生顯著變化^［4］，即氣動、傳熱和振動及其設計目標是有極強聯系的，因此渦輪葉片的熱流固耦合問題得到了廣泛的關注。

然而，即使是渦輪葉片的流熱耦合（亦稱為氣熱耦合、共軛換熱）、流固耦合這樣的兩兩耦合也存在著極大的難點。流熱耦合面臨動靜干涉^［5］和強湍流條件下流體對流換熱和固體導熱之間時間尺度不匹配的問題。對于后者，流體力學和固體力學協同求解的方式難以統一。

在流熱耦合方面，Sondak等^［6］最早意識到了非定常流熱耦合的時間尺度不統一的問題。Hwang等^［7］在流體時間尺度上研究了動靜干涉條件下靜葉尾緣冷氣對下游動葉換熱的影響，表明固體的溫度幾乎保持在定值。Rahaim等^［8］使用耦合的有限體積/雙重互易邊界元方法求解了共軛耦合傳熱問題。Hernández等^［9］從較大的時間尺度（1000s）上分析了帶熱障陶瓷涂層（TBC）的葉片從冷態進入熱態的過程。Makhija等^［10］提出使用與流體域時間常數相當的物性參數與替代的勢流方程進行計算。Schmidt等^［11］對流場進行有限容積法的非定常計算，對固體域進行有限元定常計算，在流固邊界交換數據。He等^［12］量化了對流和傳導之間在時間尺度上的巨大差異，以及其對穩態和非穩態共軛解的影響，并采用時域流體解和頻域固體傳導之間的耦合混合方法來重新調整嚴重不匹配的時間尺度。上述學者的研究，使得非定常流熱耦合分析得到了較好的發展，為熱流固耦合分析獲得葉片的溫度分布提供了有力支撐。

在涉及流體力學和固體力學協同求解時，從耦合機理上講，研究固體和流體采用的坐標系不同，在固體的求解中習慣采用Lagrange坐標系，著眼于質點；而在流體的求解中多使用Euler坐標系，著眼于空間點。因此，如何將兩種坐標系協調統一以完成求解耦合問題，是流固耦合以及熱流固耦合中不可回避的一個難點。周盛^［13］以葉片上70%～75%展向處葉片截面為研究對象，將固體力學部分簡化為平移扭轉彈性模型，分析了二維流固耦合問題。過去數十年間，隨著商業軟件的發展，通過求解器之間的數據交換，使得流體和固體之間以弱耦合的方式完成數值協同求解成為可能。Davison等^［14］詳細研究了流體分析、熱傳導分析和應力分析的耦合分析方法。Dawes等^［15］以氣冷渦輪葉片和機匣的壽命預測為例，進行了氣動、熱傳導和應力等控制方程的求解以及耦合分析。Collard等^［16］研究了動靜干涉以及葉片振動情況。全金樓等^［17］、周迪等^［18］均使用計算流體動力學（CFD）和計算結構動力學（CSD）協同求解的方法分析了葉片顫振問題。在計算準確度和速度得以保證的前提下，使用跨求解器的弱耦合是一種協調統一行之有效的工程求解方法。

部分研究者嘗試將原本雙向強耦合的問題朝著單向弱耦合的方向簡化。Rademaker等^［19］研究了中型商用發動機高壓渦輪機匣的熱負荷及冷卻通道設計。謝永慧等^［20］采用熱流固耦合方法獲得了葉片溫度分布和應力應變分布的具體信息。黃朝暉等^［21］研究了液體火箭發動機渦輪泵的葉片裂紋。湯文章^［22］對帶熱障涂層的渦輪葉片的溫度場和應力場進行數值模擬。除了商業軟件外，也有不少研究者嘗試自己開發耦合計算程序，其中Guo等^［23］基于有限差分法開發了渦輪熱流固耦合仿真計算平臺，李濤^［24］開發了一套氣熱彈耦合計算程序。盡管單向弱耦合的數值求解如今在部分領域能得到令人滿意的解^［25-27］，但大多數所謂的耦合研究僅僅只是將其中一個物理場作為邊界條件單向輸入從而進行孤立的數值求解，固體葉片發生的形變量以及溫度重新分布并沒有進一步返回到流場，從而無法形成熱、流、固三者耦合的閉環，以至于流場在微小變形下發生的改變對后續帶來的影響并不可知。陳剛等^［28］明確指出，耦合必定是兩者的相互作用。所謂單向耦合這種單方面作用的形式是否能被稱之為耦合受到了質疑。兩場耦合并沒有得到良好充分的發展，反而類似于MPCCI、Ansys system coupling等用于數據交換和傳遞的商業軟件開始逐步出現且升級。

然而，熱、流、固耦合在流熱耦合依然存在不少問題的情況下，還引入了另外一個物理場，使問題進一步復雜化。原本流熱耦合已經存在時間尺度不匹配問題，葉片的往復振動又將引入一個新的時間尺度，如何平衡上述時間尺度，成了一個新的難點。為了得到準確的傳熱求解，通常流體域棱柱網格的第一層節點需要深入黏性底層，但葉片的振動使得保持合理的精細化動網格質量和穩定求解變得尤為困難。而在諸多單向耦合的數值分析中，并沒有直面雙向耦合過程中協同求解的核心問題。

針對上述問題及其復雜物理現象，本文應用作者獨立開發的邊界層動網格技術和雙重運動方法^［29-30］，構建了渦輪葉片熱流固耦合數值計算方法，計算評估了求解器的數據傳遞誤差，并應用于GE-E3發動機第一級高壓渦輪葉片，實現了工程意義上的熱流固耦合計算與分析，進一步給出了多物理場耦合的情況下帶陶瓷涂層的該渦輪葉片的振動和換熱特性。

1 熱流固耦合計算方法

1.1 耦合計算基礎框架

在熱流固三者耦合的計算框架中，固體內的導熱是在有限元求解器中完成，即流體與固體之間的共軛換熱屬于弱耦合。相比于強耦合，弱耦合的缺點是耦合傳熱需要在兩個求解器之間迭代，增大了計算工作量。優點是在計算熱應力前，溫度數據是以面的形式而不是體的形式進行傳遞和插值，減小了插值的計算量和誤差，也避免了有限容積法求解中動網格需要同時兼顧流體與固體兩個區域的網格質量的束縛。

熱流固耦合的計算需要同時使用有限容積求解軟件Fluent、有限元求解軟件Abaqus和數據交換軟件MPCCI。具體是利用Fluent求解流體域，利用Abaqus求解固體域，兩者通過MPCCI聯立起來并交換數據。圖1（a）給出了熱流固耦合基本框架。在一個時間步內，Fluent和Abaqus分別完成內迭代之后再交換數據，如此最多循環5次進入下一個時間步。圖1（b）為時間步迭代方法，展示了時間步、動靜干涉下的掃掠周期以及整體數值計算過程的構成關系。為了獲得求解精度與求解速度的平衡，在計算中各個參數取值通常為：單求解器迭代步數nsolid=10，nfluid=40，耦合迭代步數k=5，時間步數m=25，周期數l=20，按照總迭代步數計算式 （nsolid + nfluid）×k×m×l，則整個熱流固耦合計算總共需要計算迭代至多125000步。

1.2 雙重運動方法

為了使得形變傳輸精確，本文未使用常規的滑移網格模型來實現動靜葉之間的相對運動。其原因在于，在滑移網格模型下，動葉葉柵計算域網格會在每一個時間步進行一次移位，其數量級為10mm。在每一次流固之間進行數據交換時，葉片形變使流體計算域產生的網格位移，其數量級為1×10^－3mm。最終的網格節點位移為兩者之和，形變位移在計算過程中極易被忽略。因此，本文提出了一種新的計算方案，即雙重運動方法，來模擬動靜葉之間的相互交錯。所謂雙重運動方法，就是靜葉同時使用動坐標系模型（frame motion）和滑移網格模型（slide mesh motion），而動葉使用動坐標系模型。雖然靜葉和動葉同時采用動坐標系模型，但反向的滑移網格模型使得靜葉相對于絕對坐標系依然是靜止狀態。而動葉網格的網格不滑移，其旋轉運動由動坐標系模型實現，最終計算域的網格只有來自于結構變形產生的位移。

雙重運動系統下流動規律的正確性已經通過理論推導得到了證實。流體的真實速度u、流體在動坐標系下的速度ur、流體在雙重運動系統下的合成速度uDMS之間有如下關系

ur=u－uframe（1）

uDMS=ur－umesh（2）

式中：下標frame表示動坐標系的相關變量；mesh表示動網格的相關變量。

特別地，在雙重運動系統中，動坐標系角速度、動網格角速度數值相同而方向相反，即

ωmesh=－ωframe（3）

將上述關系代入至旋轉坐標系下的相對動量方程

ρDurDtr+ρ（2ωframe×ur+ωframe×ωframe×r）=

－Δp+Δ·τ+ρ·f（4）

則能得到如下表達式

ρDuDMSDt=－Δp+Δ·τ+ρ·f（5）

式中：ρ是流體密度；p是壓力；τ為剪切力張量；f為徹體力。上式與絕對靜止坐標系下的動量方程完全一致，除此之外的連續性方程和能量方程同理，從而證明了雙重運動系統和絕對靜止坐標系下流動的等價性。

總之，雙重運動方法在模擬葉片真實轉動和相互干擾進行流場數值計算時，在流固數據傳遞中確保了動葉微小變形的數據真實性和精確度。

2 計算模型與邊界條件

2.1 金屬基、陶瓷涂層及流體域

本文以GE-E3高壓渦輪第一級靜動葉排為研究對象，對動靜交錯環境下帶陶瓷涂層的渦輪葉片進行熱流固耦合研究。對于靜葉模型，將氣膜孔排簡化為槽縫處理，即在葉片吸力面前部、壓力面中后部繪制模擬氣膜孔的槽縫，在尾緣處繪制劈縫，三者均含有氣膜通道，以降低網格量和計算量。動靜葉數分別取76和38，即一個靜葉和兩個動葉形成對應關系。根據GE-E3的原始設計，動葉內部帶有冷卻通道，動葉表面繪制0.1mm厚的熱障陶瓷涂層，葉尖部分采用凹槽。其中，陶瓷涂層通過向固體側按當地法向方向偏置曲面得到，以保證葉片外部幾何完全相同。鑒于目前國內的陶瓷涂層噴涂技術和工藝，葉片的葉尖及其凹槽沒有陶瓷涂層覆蓋。圖2（a）展示了熱流固耦合研究的渦輪級幾何模型，圖2（b）所示為葉片內冷通道、金屬基及表面的陶瓷涂層模型。

2.2 邊界條件及參數

靜葉的壓力面、吸力面氣膜孔槽縫和尾緣劈縫均為流量進口邊界，流量和總溫分別為0.0336kg/s和960K、0.0332kg/s和980K以及0.0224kg/s和987K。冷氣設置為垂直邊界面進入。表1展示了主流進出口邊界條件。

通過Abaqus求解固體導熱，涂層和金屬基形成裝配體，兩者接觸面設定綁定約束，葉根施加鉗定約束。使用動力-溫度-位移耦合的顯式求解器，設定額定轉速為1322.6rad/s時非慣性參考系下離心力形成的預應力場為初始應力場，設定1143K為初始應力溫度場。陶瓷涂層表面、陶瓷葉尖側面、金屬基葉尖面分別賦予來自流體的氣動載荷。動葉內部通道依據其氣熱耦合計算，分別設定恒定的內冷換熱系數，以模擬冷氣在蛇形通道中帶走內壁上熱量的過程。K405高溫合金葉片的彈性模量、剪切模量、泊松比和線膨脹系數均隨溫度變化而變化。陶瓷涂層簡化為一層計算，共彈性模量取87GPa，泊松比為0.12，密度為6150kg/m3，線膨脹率為12.16×10^－6K^－1。

3 熱流固耦合數值方法評估

有關熱流固耦合實驗數據的報道，在公開文獻中非常少。為了保證熱流固耦合計算結果的可靠性，必須同時確保流熱耦合的準確性和結構力學求解的準確性。在結構變形和氣動載荷耦合方面，數值求解準確性借助跨聲速機翼顫振的實驗數據得到了驗證^［31］。本文在流熱耦合的基礎上進行網格無關性驗證和湍流模型驗證，以提高流場和溫度場計算的準確性。借助某航發全尺寸、全氣膜中溫中壓渦輪動葉中葉展的溫度分布^［32］進行數值驗證，如圖3所示，實驗點的數字為測溫熱電偶的編號。雖然計算結果在壓力面尾緣區域差距較大，但k-ω模型與實驗結果吻合最好，數值計算結果的最大誤差不足5%，均方根誤差為1.2%。網格無關性亦使用該實驗數據進行評價，網格的擬合外推誤差為2.064%。

本文的流體域使用四面體網格，由Ansys Meshing繪制，在壁面上繪制棱柱邊界層網格，邊界層第一層網格厚度為0.0016mm，保證葉片表面的近壁面y+在1附近。流體域動葉部分網格數為571萬，靜葉部分網格數為384萬。對于固體域，通過掃掠的方法繪制一層四棱柱陶瓷涂層網格，使用三棱錐繪制葉片金屬基網格。葉片金屬基的網格數為33萬，陶瓷涂層的網格數為1.4萬。圖4所示為整體網格和局部細節。流體域的計算采用k-ω湍流模型。

3.1 數據傳遞誤差評估

計算過程中監視了數據傳遞的殘差，如圖5所示。除了最開始時間步的殘差較高外，剩下時間步內的殘差均有下降，其中流體和固體之間節點傳遞的位移、氣動力誤差均方根不超過0.1%，溫度數據的均方根殘差保持在0.001%以下。

3.2 動網格質量評估

計算過程中發現，在對固體施加至流體的位移增量進行積分后，其位移總量達到了3×10^－4m數量級。這一位移量為陶瓷涂層厚度的幾倍，在真實情況下振動幅度肉眼可見；在流體計算域中這一變形量也已大幅度超過了邊界層棱柱網格的總厚度。為了保證流體網格質量，使用動網格技術，即在劃分網格時將棱柱網格和四面體網格分離到兩個計算域，賦予邊界層棱柱網格和壁面相同的運動。在整個計算過程中，最差網格質量為0.89，經過一個葉片振動的往復運動之后，最差網格質量有所好轉并穩定在0.86附近，網格平均質量保持0.24不變。計算過程中并沒有出現因為邊界運動產生負網格而計算失敗或網格質量過低，導致單個有限體積內計算精度過低而浮點溢出的情況。圖6為計算過程中保存葉片振動至總體位移最大處時的cas文件所觀察到的網格細節。

從圖6可見，葉片附近的網格依然保持著良好的形態，沒有出現可觀的大面積畸變。由此可見，本文發展的細分到黏性底層的低y+值邊界層網格的高容錯率動網格技術，在實際的數值計算中是切實有效且可行的。這種方法解決了諸多動網格技術在數值工程中應用的難點，提升了計算穩定性，避免了計算中途出錯帶來的資源和時間的浪費。

4 結果分析

4.1 熱流固耦合場下葉片應力應變分析

圖7為葉片隱藏涂層后其金屬基部分在兩個掃掠周期內不同時刻的von Mises等效應力與變形，為了更明顯觀察葉片的振動形態，將位移整體放大了20倍。可以看出，在葉片的3個位置處出現了較大的應力，分別是葉尖肩部中弦區域、葉根尾緣和內部通道葉根吸力面。

為了研究傳熱和熱應力在整個熱流固耦合系統中的作用機理和重要性，設置了流固對照組。流固對照組在數據交換中去掉熱流密度和溫度的傳輸及插值，在固體域直接賦予恒定溫度1143K，原本的熱流固耦合計算簡化為流固耦合計算。與圖7相似，圖8為對照組兩個掃掠周期內不同時刻葉片的von Mises等效應力。

除了上述帶有形變的應力云圖外，圖9和圖10分別給出了熱流固耦合計算及流固耦合對照組的葉根前后緣及高應力區域的應力變化情況。

在對照組中，除離心力外其主要應力來源于氣動載荷，較大應力區域在葉根前后緣、尾緣和吸力面。在除去第一個周期內的振動后，其最大值分別為108.30、67.29MPa，比考慮傳熱的數值結果要小一個數量級。

這是由于當加入流熱耦合中的換熱，其應力應變的主要來源從單一的氣動力轉變成氣動載荷和熱載荷的綜合作用。對于葉片葉肩部分（固體），在泄漏流和刮削渦的影響下，流體對固體形成了一個強換熱效果，葉肩一側有熱障涂層，頂面和凹槽內側由于直接與高溫燃氣接觸，底面通過金屬熱傳導傳出熱量，在這樣的環境下其溫度梯度較大，形成了較大的熱應力，其時均值為752.1MPa。對于葉片葉根部分，內部通道溫度低于無預應力初始溫度，而葉片表面溫度高于無預應力初始溫度，形成了葉型平面內的溫度差異；受端壁附近的馬蹄渦的影響，其外部的熱負荷與中葉展區域的熱負荷存在差異，形成了葉高方向上的溫度差異，兩者一起形成了葉根部分的熱應力。

需要特別指出的是：真實情況下渦輪葉片同時承受氣動載荷和熱載荷，兩者無法剝離獨立，而本文計算中通過關閉數據傳遞的方式進行對照，考慮氣動載荷的獨立作用才呈現出了數量級的差異。其次，渦輪動葉還有一個非常重要的載荷來源即離心力，當考慮離心力時，在葉根處呈現出的應力亦在百兆帕的量級。然而在轉速恒定時，離心力形成的應力分布也是恒定的。本文主要是分析動靜干涉下非定常氣動載荷和熱載荷的影響，根據結構力學中的線性疊加原理，在分析上述載荷時都統一約去了離心力，以便能更好地分析兩種載荷對于渦輪葉片的作用。

圖11展示了兩個掃掠周期內不同時刻下葉片的變形， 葉片各節點上的矢量箭頭長度和顏色代表當地的位移，葉片整體變形形態同樣放大了20倍。依照同樣的視角和放大倍率，繪制了如圖12所示的對照組兩個掃掠周期內不同時刻葉片的總變形。在熱流固耦合仿真結果中，振動穩定后葉片的最大位移為2.09×10^－4m。在流固耦合對照組中，計算一開始，由于氣動力直接施加在葉片上，葉片在兩個掃掠周期內達到最大位移為8.33×10^－5m，在流體形成的阻尼耗散消除初場的影響后，其最大位移為4.95×10^－5m。

對于考慮傳熱的熱流固耦合計算，其最大位移比流固耦合對照組高出了一個數量級。如果將熱流固耦合的計算結果與對照組的結果畫在同一張圖里，則由于其差異太大無法觀察出其波動情況。因此，對于葉尖處，分別繪制了如圖13和圖14所示的兩組前緣和尾緣點以及剛心的位移。在熱流固耦合研究中，前緣點和尾緣點的位移相近。而在對照組中，氣動載荷產生的負縱向力矩使得前緣點和尾緣點的位移均值并不相同，前緣點比尾緣點低0.025mm。

綜上，在耦合場中是否考慮傳熱，不僅會導致應力結果的數值和分布的差異，而且在變形、整體位移上也存在差異。相比于流固耦合計算，盡管熱流固耦合計算時固體域的平均溫度和流固耦合計算時賦予的溫度相等，但在前者計算時，葉片由于外部燃氣施加熱負荷，內部通道帶走熱量形成葉片外部溫度高內部溫度低的情況，從而葉片的材料模量隨溫度變化，形成了葉片外部模量低、內部模量高的情況。同時，由于壓力面溫度高于吸力面，其熱膨脹量高于吸力面，形成了一個與氣動載荷相同的使葉片向吸力面彎曲的趨勢。因此，考慮固體內溫度分布和熱應力的熱流固耦合數值結果，在整體位移和變形上要高于流固耦合的數值結果。

4.2 熱流固耦合場下葉片振動形態分析

對于熱流固耦合和流固耦合兩組位移分別做快速傅里葉變換，其結果如圖15和圖16所示。當頻率為4500Hz和7314Hz時，熱流固耦合存在兩個峰值，分別與葉片在流固耦合對照組定常氣動力和離心力作用下的一階、二階模態接近。熱流固耦合所形成的葉片各處模量的頻率和對照組有略微的差別。葉片在耦合情況下呈現彎曲和扭轉兩種形式組合振動，在葉片扭轉變形并不引起剛心的位移的前提下，為了分離兩種振動形態單獨研究，同樣也提取并分析了葉尖剛心的位移，如圖中黑色線所示。

隨時間變化，剛心的位移表現了葉片的彎曲振動情況，在對照組中，剛心的位移隨時間變化存在兩個峰值，第二個峰值與動靜干涉即動葉被掃掠過的頻率7999Hz吻合，而在熱流固耦合的條件下僅在葉片一階彎曲處存在一個峰值。因此，熱流固耦合中熱應力是影響其主要變形的因素，氣動載荷對其形變的貢獻次之，而在對照組得到的振動對氣動載荷這一主要因素的響應更為明顯。對照組前后緣振動的第二個峰值的頻率為7339Hz，與動靜干涉的頻率更為接近，即受交變氣動載荷的驅動，其在第二個峰值的分量更大。

4.3 熱流固耦合場下葉尖位移對流動和換熱的影響分析

固體葉片造成的網格運動傳遞至流體邊界時將產生網格運動和變形，其中葉尖最大位移在0.3mm。在許多葉片振動的簡化研究中，通常只對固體葉片施加氣動載荷而不考慮固體的變形對流場的影響。通過三維熱流固耦合的計算結果，可說明微小變形在耦合計算中的重要性。

圖17提取了不同時刻90%葉高截面處葉片表面的壓力分布。通過保證動靜葉均處于相同的相位，來排除動靜葉不同的相對位置對動葉表面氣動力的影響。可以明顯看出，不同時刻葉片表面上靜壓的差異，即葉片微小的變形能帶來氣動上可觀的變化。圖18所示為熱流固耦合數值結果所得的不同時刻下90%葉高截面處葉片表面的熱流密度分布，可見除了氣動上的差異外，微小的變形也帶來了明顯的傳熱差異。總之，以往不考慮固體的變形對流場影響的簡化研究與真實情況是存在差異的，尤其是在考慮葉片換熱和熱應力的耦合計算中，流體的邊界運動和計算域變形是極其重要的。

4.4 熱障涂層對葉片振動特性影響

熱障陶瓷涂層會影響葉片整體耦合振動，主要原因有：①涂層阻隔了一部分熱負荷，降低了葉片的整體溫度，金屬材料的力學特性隨著溫度變化從而其振動特性發生改變；②涂層的密度、模量和熱膨脹系數與高溫合金不盡相同，涂層會影響葉片整體的質量和剛度從而改變其振動特性。上述兩者都會進一步影響流場的變化，如此循環造成整個熱流固耦合特性的變化；③陶瓷涂層的模量比金屬基低，其密度也低于金屬基，根據一維無阻尼自由振動公式，其彈性模量和質量的改變會引起固有頻率變化。為此，本小節研究在三維復雜熱流固耦合環境下，陶瓷涂層彈性力學特性對葉片振動的影響規律。

本節同樣設置一個對照組，即把陶瓷涂層的彈性力學參數設置為高溫合金的力學參數，其傳熱參數依舊保持不變。對照組的應力分布如圖19所示。可見，相比對照組，真實涂層參數條件下高溫金屬的模量大于陶瓷涂層，高應力大部分是由金屬基承擔，但對照組的陶瓷涂層也會呈現較大的應力。

如圖20所示，對照組與真實涂層的熱流固耦合的振動特性十分接近，如果繪制形變云圖或者位移線圖，兩者幾乎看不出任何區別。經過對位移數據進行快速傅里葉變換，振動的兩個峰值頻率分別為4523Hz和7327Hz，比真實頻率分別高出0.51%和0.18%，可見其差異十分微小。

因此，在葉片表面噴涂陶瓷涂層，其核心的影響機理是陶瓷涂層改變了葉片的耦合換熱特性和金屬基的溫度分布，形成了不一樣的熱應力分布，從而引起耦合條件下振動特性的改變。由于陶瓷涂層本身非常薄，單從彈性力學特性方面看，其對葉片熱流固耦合條件下振動特性的影響極小。

5 結 論

本文以GE-E3高壓渦輪第一級靜葉和動葉為研究對象，應用動網格方法和雙重運動方法，構建了熱流固耦合數值模擬基本框架，并對動靜交錯環境下帶陶瓷涂層的渦輪葉片進行了熱流固耦合數值分析與研究，主要結論如下。

（1）在熱流固耦合計算中，數據的跨求解器傳遞流體和固體之間節點傳遞的位移、氣動力誤差均方根不超過0.1%，溫度誤差均方根不超過0.001%。動葉在流固耦合數值計算中只考慮氣動載荷及變形，其最大位移為4.95×10^－5m。在熱流固耦合中，同時考慮氣動載荷熱負荷及變形，其最大位移為2.09×10^－4m。兩種耦合下的振動形態均為一階彎曲及扭轉振動形態的混合，其響應頻率與其固有頻率接近。除了總變形存在著數量級的差異以外，熱流固耦合的應力分布和單流固耦合也不同，其應力的主要來源是熱應力，葉根前緣應力為496.9MPa，亦比流固耦合高一個數量級。

（2）動葉在氣動載荷和熱載荷下的變形會進一步影響作用在葉片上的氣動力和熱負荷。由于其變形，葉尖附近的靜壓變化達87kPa（相較其當地靜壓占比9.44%），熱流密度變化達4.6×105W/m2（相較其當地熱流密度占比39.7%），即在考慮葉片換熱和熱應力的耦合計算中，葉片振動產生的微小邊界運動與計算域變形是極其重要的。

（3）陶瓷熱障涂層對葉片振動影響的核心機理是涂層改變了葉片的耦合換熱特性和金屬基的溫度分布，形成了不一樣的熱應力分布，從而引起其耦合條件下振動特性的改變，單從彈性力學特性方面看，這對葉片熱流固耦合條件下振動特性的影響極小。

真實的動葉通過樅樹型葉根固定在輪盤上，若要計算得到葉片的真實振動結果，未來需要考慮葉片端壁的幾何構型、葉根與輪盤的接觸，以及輪盤和轉子軸構成的整個系統。此外，使用熱流固耦合研究可以同時準確地得到渦輪葉片上各處的溫度和應力，從而預測相應時間內的蠕變量，未來應開展基于渦輪葉片熱流固耦合的蠕變及疲勞壽命研究。

參考文獻：

［1］陳矛章. 風扇/壓氣機技術發展和對今后工作的建議［J］. 航空動力學報， 2002， 17（1）： 1-15.

CHEN Maozhang. Development of fan/compressor techniques and suggestions on further researches［J］. Journal of Aerospace Power，2002， 17（1）： 1-15.

［2］張明明， 李紹斌， 侯安平， 等. 葉輪機械葉片顫振研究的進展與評述［J］. 力學進展， 2011， 41（1）： 26-38.

ZHANG Mingming， LI Shaobin， HOU Anping， et al. A review of the research on blade flutter in turbomachinery［J］. Advances in Mechanics， 2011， 41（1）： 26-38.

［3］吳錦濤， 王珺， 徐自力， 等. 高轉速部分進氣渦輪盤氣流力及葉片振動響應研究［J］. 西安交通大學學報， 2022， 56（7）： 108-117.

WU Jintao， WANG Jun， XU Zili， et al. Airflow force and vibration response for high speed partial-admission turbine disk［J］. Journal of Xi’an Jiaotong University，2022， 56（7）： 108-117.

［4］《中國航空材料手冊》委員會. 中國航空材料手冊［M］. 北京： 中國標準出版社， 2002.

［5］GILES M B. Stator/rotor interaction in a transonic turbine［J］. Journal of Propulsion and Power， 1990， 6（5）： 621-627.

［6］SONDAK D L， DORNEY D J. Simulation of coupled unsteady flow and heat conduction in turbine stage［J］. Journal of Propulsion and Power， 2000， 16（6）： 1141-1148.

［7］HWANG S， SON C， SEO D， et al. Comparative study on steady and unsteady conjugate heat transfer analysis of a high pressure turbine blade［J］. Applied Thermal Engineering， 2016， 99： 765-775.

［8］RAHAIM C P， KASSAB A J， CAVALLERI R J. Coupled dual reciprocity boundary element/finite volume method for transient conjugate heat transfer［J］. Journal of Thermophysics and Heat Transfer， 2000， 14（1）： 27-38.

［9］HERNNDEZ R A， MAZUR C Z， BAUTISTA P E A， et al. Unsteady 3-D conjugated heat transfer simulation of a thermal barrier coated gas turbine bucket［C］//ASME Turbo Expo 2008： Power for Land， Sea， and Air. New York， USA： ASME， 2009： 2413-2423.

［10］MAKHIJA D S， BERAN P S. Time scale effects in topology optimization of the interior material distribution of a body subject to transient conjugate heat transfer［C］//AIAA Scitech 2019 Forum. Reston， VA， USA： AIAA， 2019： AIAA 2019-1468.

［11］SCHMIDT M， STARKE C. Comparison of steady and unsteady coupled heat-transfer simulations of a high-pressure turbine blade［C］//ASME Turbo Expo 2015： Turbine Technical Conference and Exposition. New York， USA： ASME， 2015： V05AT10A016.

［12］HE L， OLDFIELD M L G. Unsteady conjugate heat transfer modeling［J］. Journal of Turbomachinery， 2011， 133（3）： 031022.

［13］周盛. 葉輪機氣動彈性力學引論［M］. 北京： 國防工業出版社， 1989： 65-74.

［14］DAVISON J B， FERGUSON S W， MENDONA F G， et al. Towards an automated simulation process in combined thermal， flow and stress in turbine blade cooling analyses［C］//ASME Turbo Expo 2008： Power for Land， Sea， and Air. New York， USA： ASME， 2008： 2641-2648.

［15］DAWES W N， KELLAR W P， HARVEY S A. Towards cooled turbine preliminary life prediction via concurrent aerodynamic， thermal and material stress simulations on conjugate meshes［C］//ASME Turbo Expo 2010： Power for Land， Sea， and Air. New York， USA： ASME， 2010： 711-722.

［16］COLLARD J E III， CIZMAS P G A. Blade-forced vibration effects on turbomachinery rotor-stator interaction［C］//37th Joint Propulsion Conference and Exhibit. Reston， VA， USA： AIAA， 2001： AIAA 2001-3472.

［17］全金樓， 張偉偉， 蘇丹， 等. 基于CFD/CSD時域耦合方法的多通道葉柵顫振分析［J］. 航空學報， 2013， 34（9）： 2019-2028.

QUAN Jinlou， ZHANG Weiwei， SU Dan， et al. Flutter analysis of turbomachinery cascades based on coupled CFD/CSD method［J］. Acta Aeronautica et Astronautica Sinica， 2013， 34（9）： 2019-2028.

［18］周迪， 陸志良， 郭同慶， 等. 基于CFD/CSD耦合的葉輪機葉片失速顫振計算［J］. 航空學報， 2015， 36（4）： 1076-1085.

ZHOU Di， LU Zhiliang， GUO Tongqing， et al. Stall flutter computation of turbomachinery blade based on a CFD/CSD coupling method［J］. Acta Aeronautica et Astronautica Sinica， 2015， 36（4）： 1076-1085.

［19］RADEMAKER E R， HULS R A， SOEMARWOTO B I， et al. Modeling approach to calculate redistributions of HPT-shroud cooling channels minimizing thermal stresses including some turbine blade tip effects［C］//ASME 2013 Turbine Blade Tip Symposium. New York， USA： ASME， 2013： V001T05A003.

［20］謝永慧， 鄧實， 張荻， 等. 基于熱流固耦合分析的重型燃氣輪機透平高壓葉片壽命研究［J］. 熱力透平， 2011， 40（3）： 151-158.

XIE Yonghui， DENG Shi， ZHANG Di， et al. Study on life prediction for HP turbine blade of heavy-duty gas turbine based on analysis of multi-factor coupling physical fields［J］. Thermal Turbine， 2011， 40（3）： 151-158.

［21］黃朝暉， 袁奇， 張弘斌， 等. 某型火箭發動機渦輪轉子流熱固耦合強度及疲勞壽命分析［J］. 西安交通大學學報， 2022， 56（8）： 73-84.

HUANG Chaohui， YUAN Qi， ZHANG Hongbin， et al. Analysis of fluid-thermal-solid coupling strength and fatigue life of a certain rocket engine turbine rotor［J］. Journal of Xi’an Jiaotong University， 2022， 56（8）： 73-84.

［22］湯文章. 基于流固耦合熱障涂層渦輪葉片應力場的數值模擬［D］. 湘潭： 湘潭大學， 2015.

［23］GUO Zhaoyuan， WANG Qiang， DONG Ping， et al. The key techniques for thermal-flow-elastic coupling numerical simulation platform in turbines［C］//ASME 2008 International Mechanical Engineering Congress and Exposition. New York， USA： ASME， 2008： 1077-1083.

［24］李濤. 基于非結構網格的氣冷渦輪氣熱彈耦合數值計算［D］. 哈爾濱： 哈爾濱工業大學， 2014.

［25］李立州， 王婧超， 呂震宙， 等. 渦輪葉片流固耦合協作優化方法［J］. 推進技術， 2008， 29（5）： 604-608.

LI Lizhou， WANG Jingchao， L Zhenzhou， et al. Collaborative optimization method for fluid-solid coupling problem of turbine blades［J］. Journal of Propulsion Technology， 2008， 29（5）： 604-608.

［26］XU Ning， LIU Zhansheng， PENG He. Transient thermal stress analysis of the typical structures in turbine rotor based on thermomechanical analysis［C］//ASME Turbo Expo 2014： Turbine Technical Conference and Exposition. New York， USA： ASME， 2014： V05CT20A006.

［27］艾延廷， 劉仲堯， 王志. 起動階段渦輪導向葉片熱應力及冷卻作用研究［J］. 機械設計與制造， 2015（1）： 97-100.

AI Yanting， LIU Zhongyao， WANG Zhi. Research of thermal stress and cooling effect on a turbine guide vane during the start-up period［J］. Machinery Design amp;Manufacture， 2015（1）： 97-100.

［28］陳剛， 李躍明. 非定常流場降階模型及其應用研究進展與展望［J］. 力學進展， 2011， 41（6）： 686-701.

CHEN Gang， LI Yueming. Advances and prospects of the reduced order model for unsteady flow and its application［J］. Advances in Mechanics， 2011， 41（6）： 686-701.

［29］汪翔宇， 豐鎮平. 渦輪葉片振動對其換熱影響的數值研究［C］//中國工程熱物理年會2017年學術年會論文. 2017： 20172146.

［30］西安交通大學. 葉輪機械流動仿真計算方法及系統： CN202311167825.8［P］. 2023-12-12.

［31］WANG X Y， FENG Z P. Study on vibration characteristics of rotor blades and damping scheme with film cooling holes under rotor-stator interaction［C］//The 15th International Symposium on Ultrasonic Applications and Advanced Acoustic Technologies. Oxford， UK： ISUAAAT， 2018： ISUAAAT15-085.

［32］楊星， 汪翔宇， 豐鎮平. 渦輪葉片綜合冷卻有效度評價指標研究［J］. 工程熱物理學報， 2020， 41（7）： 1642-1648.

YANG Xing， WANG Xiangyu， FENG Zhenping. A study on scaling method of turbine blade metal temperatures［J］. Journal of Engineering Thermophysics， 2020， 41（7）： 1642-1648.

（編輯 亢列梅）