梁賢燁,弭光寶,李培杰,黃 旭,曹春曉

（1.中國航發(fā)北京航空材料研究院 鈦合金研究所, 北京 100095；2.清華大學 新材料國際研發(fā)中心, 北京 100084；3.中國航發(fā)先進鈦合金重點實驗室, 北京100095）

為滿足現(xiàn)代航空發(fā)動機輕量化結(jié)構(gòu)設(shè)計的要求，具有高比強、耐高溫等優(yōu)異特性的鈦合金得到大量應用[1-2]。然而，隨著先進航空發(fā)動機性能的進一步提升，高壓壓氣機內(nèi)部的運行工況更加復雜和苛刻，包括異物撞擊、刮磨、喘振、碰磨和熔滴等在內(nèi)的異?，F(xiàn)象突發(fā)的可能性大大增加，導致鈦合金葉片和機匣等部件局部溫度急劇升高而發(fā)生著火事故，即鈦火。鈦火的特點是，起火時間短、燃燒傳播速度快，具有突發(fā)性和不可預測性，難以采取有效的滅火措施，輕則造成葉片燒損、機匣燒穿，重則燒毀整個發(fā)動機。鑒于鈦火的巨大危害性，鈦火防控（含實驗測試、材料和表面等[3-5]的綜合技術(shù)）引起高度重視，并成功應用于F119和EJ200等國外著名航空發(fā)動機，遺憾的是國內(nèi)尚處于應用基礎(chǔ)研究階段，尤其對鈦火機理及邊界特性不清楚。

國外學者對不同加熱形式的金屬著火模型進行理論與實驗研究。Classman等[6]和Frank等[7]認為金屬著火屬于非均相氧化反應，并在熱自燃理論基礎(chǔ)上建立了金屬著火模型；Khaikin等[8]通過研究著火參數(shù)與氧化層變化的相關(guān)性，建立-金屬氧化層著火模型；Merzhanov等[9]在研究催化反應過程中指出等溫著火模型的局限性，建立了非等溫加熱（線性加熱）條件下的熱自燃理論，提出臨界加熱速率的概念，同時得出降低反應熱并提高散熱條件能有效降低臨界加熱速率；Barzykin[10]指出當著火延遲時間足夠長時非等溫加熱速率對著火參數(shù)的影響不能忽略；Wolf[11]針對鈦合金試件在不同加熱速率（28 K/s、58 K/s及100 K/s）下的著火過程進行實驗，建立相應的理論模型，認為模型計算結(jié)果與實驗結(jié)果的誤差是由反應過程中熱應力的產(chǎn)生所致；Sanchez等[12]改進了Frank-Kamenetskii理論并應用于線性加熱形式的自蔓延高溫合成過程，計算不同形狀燃燒室的著火參數(shù)，計算結(jié)果與實驗符合良好。

然而，以上研究采用的是簡化的零維模型，且未考慮加熱速率存在的衰減（例如在發(fā)動機實際摩擦過程中，由于摩擦因素隨著溫度升高而降低，加熱速率也會隨之衰減），不能很好的用于發(fā)動機鈦火特性的解釋。因此，本工作針對可能導致航空發(fā)動機鈦火發(fā)生的等溫加熱、非等溫線性加熱以及非等溫摩擦加熱等形式熱源，建立相應的三維有限元模型進行溫度場計算，進而得到著火參數(shù)的變化規(guī)律，同時通過與相關(guān)文獻中著火參數(shù)的實驗值對比，探討模型的不足及改進的方法，為航空發(fā)動機鈦火特性分析評估提供理論方法。

1 模型

1.1 物理模型

物理模型如圖1所示，試件上下兩個環(huán)狀表面分別代表化學反應面和加熱面，為了在子程序中記錄溫度場變化并進行相應的反應熱計算，定義上表面中的散熱面不發(fā)生反應，加熱面的熱量通過試件內(nèi)部導熱傳導至上表面，氣流溫度為環(huán)境溫度，氣流從試件上表面流過。模型針對臨界著火溫度及著火延遲時間這兩個重要參數(shù)進行計算，臨界著火溫度代表反應過程由緩慢氧化轉(zhuǎn)變?yōu)閯×胰紵臏囟?，著火延遲時間為由系統(tǒng)初始溫度至達到著火發(fā)生溫度的時間。在實際工況中，高壓壓氣機內(nèi)部的非正常升溫時間很短（譬如葉片與機匣的非正常摩擦過程至可持續(xù)燃燒時間小于10 s[13]），而散熱能力較強，提高著火延遲時間表示系統(tǒng)在未到達著火溫度之前就恢復至正常狀態(tài)，因而著火延遲時間的增加可以有效預防鈦火發(fā)生。模型假設(shè)如下：

圖1 著火模型示意圖Fig. 1 Schematic diagram of the ignition model

（1）假設(shè)在摩擦過程中持續(xù)生成微米級顆粒狀摩擦產(chǎn)物（即微凸體[3,14]），其中面容比最大的球狀微凸體首先發(fā)生著火，假設(shè)所有單位時間內(nèi)生成的球狀微凸體的集合為Npar，在微凸體與氧氣反應并產(chǎn)生氧化膜的過程中，其尺寸隨著反應的消耗逐漸減小，由于反應區(qū)域并非封閉系統(tǒng)，存在著強烈的氣流，當微凸體的半徑小于某一特定值時，微凸體在氣流的拖拽力或者試件本身的離心力作用下脫離反應表面，在此基礎(chǔ)上提出微凸體有效半徑re的概念，當微凸體半徑小于微凸體有效半徑后忽略其生成的熱量，同時由于優(yōu)先發(fā)生著火的微凸體其半徑為Npar中半徑最小者（rmin），因此假設(shè)re=rmin，其物理意義為單位時間內(nèi)微凸體放熱量為定值；

（2）如圖1所示，微凸體氧化層為TiO2，其厚度h隨著溫度的升高而增加，微凸體內(nèi)各點溫度一致；

（3）氣流流向與靜子試件平行，且流速不變，反應區(qū)的散熱通過與氣流之間對流換熱形式實現(xiàn)，忽略熱輻射的影響；

（4）忽略金屬顆粒在燃燒過程中的揮發(fā)以及熱解[15]；

（5）反應區(qū)表面滿足單分子吸附模型；

（6）放熱反應為理想狀態(tài)下鈦與氧氣充分燃燒產(chǎn)生TiO2的單步反應；

（7）在非等溫摩擦加熱過程中，摩擦因數(shù)為溫度的函數(shù)，隨著溫度的升高線性下降。

1.2 計算模型

反應面的加熱速率可通過阿累尼烏斯公式計算仿真過程中每個迭代步的熱值變化，從而得出相應的著火參數(shù)。首先根據(jù)氧化層生長理論[8,15-17]，計算相應的氧化層厚度，如式（1）所示：

式中：h為微凸體的氧化膜厚度，單位為μm，在初始條件下h=hini；K為指前因子,（kg/m2?s）；ci為氧濃度；E為微凸體的激活能，J?mol-1；Rg為普適氣體常數(shù)，J?mol-1?K-1；T為微凸體瞬時溫度，K。

初始氧化層hini可通過式（2）計算：

式中：R為微凸體半徑，μm；?為微凸體純度，為90%[15]；ρTi為Ti顆 粒 密 度，ρTiO2為TiO2密 度，kg?m3。

氧化熱QH采用式（3）計算：

式中：qr為單位質(zhì)量反應熱，MJ/kg；S為微凸體表面積，m2。

反應面微凸體能量守恒方程如下：

式中：Nu為努塞爾數(shù)；λ為導熱系數(shù)，W?m?1?K?1；m為微凸體質(zhì)量，kg；c為微凸體比熱，J?m?1?K?1；Qc為傳熱量。

等溫加熱過程中加熱面溫度方程如下：

非等溫線性加熱過程的加熱面溫度方程如下：

式中：a為加熱速率常數(shù)；t為當前時間；

非等溫摩擦加熱過程的加熱面溫度方程如下：

式中：Δt為單位時間步；為第m?1個時間步的加熱面溫度。

其生熱量QF方程[18-19]如下：

式中：B和F為摩擦因數(shù)；N為接觸應力，kPa；ω為轉(zhuǎn)子轉(zhuǎn)速，r?min?1；R1為旋轉(zhuǎn)半徑，μm。

從式（7），（8）可以看出，在摩擦加熱過程中，當前時間步的加熱溫度取決于上一個時間步的加熱面溫度，同時其加熱溫度會隨著摩擦因素的變化而逐漸衰減。

由于目前有限元軟件沒有相應的微凸體表面著火理論子程序，因此分別編制三個UDF子程序?qū)崿F(xiàn)反應面及加熱面的溫度變化，具體計算流程如圖2所示，將幾何模型、材料模型以及相應的邊界條件代入。

圖2 程序流程圖Fig. 2 Flow chart of program

同時，為了保證網(wǎng)格的精度，采用六面體模型進行網(wǎng)格劃分；接下來根據(jù)不同的加熱形式編制子程序1，并將其加載至加熱面網(wǎng)格；在流體域及固體域完成單個迭代步的溫度場計算后，通過子程序2將反應面鄰接網(wǎng)格的溫度場數(shù)據(jù)存入UDM（用戶自定義內(nèi)存）中；子程序3根據(jù)UDM中溫度場的數(shù)據(jù)進行反應熱的計算并與散熱量進行對比，若此時反應面的生熱量大于散熱量，表明已進入著火階段，則更新反應面的溫度數(shù)據(jù)，單次迭代計算結(jié)束。其中，子程序計算完成后的熱量以溫度的形式傳遞至相應的反應面與加熱面。

1.3 邊界條件

為了在保證計算流場在近壁區(qū)域具有良好收斂性的同時能夠保證足夠的精度，湍流方程采用K-E湍流模型[20]。有限元模型的節(jié)點數(shù)為268128，計算過程中瞬態(tài)時間步為0.1 s，單位時間步的迭代次數(shù)為20，環(huán)境溫度為823 K，由文獻[21-22]可知，鈦的著火溫度接近其沸點（＞3000 K），為了兼顧有限元模型的收斂性，選擇鈦合金的燃燒溫度為3000 K。邊界條件的設(shè)置包括試件著火反應區(qū)的幾何尺寸參數(shù)、環(huán)境條件以及材料的熱物性參數(shù)等。試件材料選用阻燃鈦合金，其高溫摩擦因數(shù)F為0.977，B為?6.56×10?4[23]。

表1 材料熱物性參數(shù)Table 1 Thermal property parameters of materials

表2 模型初始邊界條件Table 2 Initial boundary conditions of the model

2 結(jié)果分析

通過對有限元模型每一個時間步的計算結(jié)果進行分析，得到不同邊界條件下的溫度歷史，然后對每一個時間步的溫度場進行分析，即可得到不同時間試件的溫度分布，如圖3所示，靜子表面的溫度梯度最高處位于反應區(qū)，受到對流散熱的影響，溫度梯度隨著半徑的增加而下降。通過對溫度歷史的分析，可以得到不同邊界條件對著火參數(shù)的影響，進而得到局部最先著火區(qū)域以及熱量傳遞過程。

圖3 有限元模型溫度場分布Fig. 3 Temperature field distribution of finite element model

2.1 等溫加熱過程

作為應用最為廣泛的加熱方式，等溫加熱過程可以模擬航空發(fā)動機環(huán)境中異物撞擊或熔滴黏附于鈦合金構(gòu)件而發(fā)生的鈦火。同時，由于其數(shù)學描述最為簡單，等溫加熱也是熱自燃理論建立時所采用的方法，通過研究等溫加熱過程的溫度歷史可直觀的區(qū)分出著火和未著火區(qū)域并計算出臨界著火溫度。下面分別采用950 K、1000 K 以及1941 K（鈦熔滴溫度）作為加熱面的初始溫度進行數(shù)值模擬，分析不同條件下熱量的傳遞過程。

通過對圖3有限元模型溫度場的俯視圖及剖視圖進行分析，得到圖4和圖5所示的加熱過程的溫度場?？梢?，在加熱0.9 s時，反應面溫度為932 K，如圖4（a）所示，由于未達到臨界著火區(qū)域，其反應熱量可以忽略，因此熱量主要由加熱面提供，如圖5（a）所示；在加熱2 s時，反應面的平均溫度接近1000 K，如圖4（b）所示，由于此時達到了臨界著火區(qū)域，反應加熱量的作用表現(xiàn)明顯，而中心未參與反應的對流面溫度仍維持在960 K左右，如圖5（b）所示；當加熱時間延長至2.2 s時，加熱過程由化學反應熱主導，如圖4（c）所示，反應面的溫度超過了加熱面，達到1171 K，如圖5（c）所示；當加熱達到2.3 s時，反應面溫度達到3000 K，即達到持續(xù)燃燒溫度，如圖4（d）及圖5（d）以后進入穩(wěn)態(tài)燃燒階段，直至著火反應物或者氧化劑耗盡為止。

圖4 反應區(qū)剖視圖Fig. 4 Section view of reaction zone （a） 0.9 s; （b） 2 s; （c） 2.2 s; （d） 2.3 s

圖5 反應面俯視圖Fig. 5 Top view of reaction surface （a） 0.9 s; （b） 2 s; （c） 2.2 s; （d） 2.3 s

圖6為不同加熱溫度下加熱速率的變化情況（考慮收斂效率將950 K時的流速降低至10 m/s）?？梢?，當加熱溫度為950 K時，由于加熱溫度未達到臨界著火溫度，所以加熱時間超過3.1s后，反應面溫度維持在938 K附近不再變化，此時對流散熱量始終大于反應產(chǎn)熱量，整個加熱過程可以視為無內(nèi)熱源的穩(wěn)態(tài)傳熱過程。當加熱溫度為1000 K時，由于超過了臨界著火溫度，整個加熱過程反應面的溫度呈指數(shù)增長，反應面溫度在超過958 K之后溫度梯度明顯增加，說明臨界著火溫度約為958 K，在加熱2.3 s時達到持續(xù)燃燒溫度。當加熱溫度為1941 K時，加熱0.2 s時就已超過臨界著火溫度，在加熱0.4 s時達到持續(xù)燃燒溫度?？梢?，在高壓壓氣機工況條件下，若葉片局部著火產(chǎn)生的鈦熔滴在高速流體作用下快速擴散，在極短時間內(nèi)就會引發(fā)其他鈦合金構(gòu)件迅速達到可持續(xù)燃燒的邊界條件，即鈦熔滴擴散傳播的危險程度不容忽視。

圖6 等溫加熱溫度歷史Fig. 6 Isothermal heating temperature history

2.2 非等溫線性加熱過程

在實際工況條件下，當環(huán)境溫度為時間的函數(shù)時，加熱過程可以采用線性近似，通過采用不同加熱速率，分析不同加熱條件下試件溫度場的變化，并對加熱速率為58 K/s的實驗結(jié)果與數(shù)值計算結(jié)果進行對比。圖7是非等溫線性加熱溫度歷史的計算結(jié)果?？梢?，隨著加熱速率的增加，著火延遲時間逐漸減小，分別為1.5 s、1.1 s和0.9 s，而臨界著火溫度基本維持在950 K。該臨界著火溫度的仿真計算結(jié)果與Wolf[11]的實驗結(jié)果存在差異，主要原因在于發(fā)動機鈦火發(fā)生過程，加熱速率的不同導致構(gòu)件受到的熱應力也不同，結(jié)合微凸體理論可以推斷，微凸體尺寸會隨著加熱速率的變化而改變，因此下面將對不同微凸體有效半徑的加熱溫度歷史進行探討。

圖7 線性加熱溫度歷史Fig. 7 Linear heating temperature history

圖8是不同有效半徑線性加熱溫度歷史的計算結(jié)果。計算過程中，仿真計算模型采用靜態(tài)著火實驗條件，不考慮流速的影響，加熱速率為58 K/s，加熱溫度設(shè)置在反應面?？梢姡R界著火溫度及延遲時間隨著微凸體有效半徑的減小而下降，而在接近著火臨界條件之前，其溫度變化基本一致。當粒徑為55 μm時，臨界著火溫度約為900 K，著火延遲時間為1.3 s；當粒徑為33 μm時，臨界著火溫度約 為850 K，著 火 延 遲 時 間 為1.2 s；當 粒 徑 為16.5 μm時，臨界著火溫度約為765 K，著火延遲時間為1.1 s。文獻[11]中的臨界著火溫度為685 K，誤差值約為11.7 %，因此按照氧化層生長理論及微凸體理論，當加熱速率為58 K/s時，受到熱應力的影響其微凸體有效半徑在16.5 μm以下。由此可見，微凸體有效半徑的減小會顯著提高著火發(fā)生的概率，因此在鈦火阻燃設(shè)計中應盡量提高材料本體及表面摩擦磨損產(chǎn)生顆粒的尺寸，同時考慮從適航條件下對摩擦磨損顆粒的著火特性進行評估。

圖8 不同有效半徑線性加熱溫度歷史Fig. 8 Linear heating temperature history with different effective radius

2.3 非等溫摩擦加熱過程

在實際的高壓壓氣機工況中，加熱形式以葉片與機匣之間的非正常高速摩擦產(chǎn)熱為主，與線性加熱方式相比，由于摩擦因素隨著溫度升高而近似線性降低，因此其加熱速率也會相應衰減。下面分別對不同接觸應力、流速以及氧濃度的加熱過程進行計算分析。需要說明的是，當反應生熱量數(shù)值超過摩擦生熱量時，雖然會增加摩擦因素的衰減，然而此時反應面溫度已經(jīng)處于臨界著火溫度之上，達到持續(xù)燃燒溫度的時間僅為0.1～0.3 s，摩擦產(chǎn)熱的影響可以忽略，因此在計算分析過程中主要對到達著火溫度之前的加熱過程進行理論分析。

圖9是不同接觸應力摩擦加熱溫度歷史的計算結(jié)果。可見，臨界著火溫度維持在950 K左右，著火延遲時間隨著接觸應力的增加而下降，當接觸應力為2.65 kPa時，到達臨界著火溫度之前反應面的加熱速率維持在25 K/s，著火延遲時間為6.2 s；當接觸應力為13.25 kPa時，維持在70 K/s左右，著火延遲時間降至1.9 s；而當接觸應力提升10倍，達到26.5 kPa時提高至130 K/s左右，著火延遲時間降至1.4 s?？梢?，隨著接觸應力的提高發(fā)生鈦火的概率也會大幅增加，因而在航空發(fā)動機的設(shè)計過程中應考慮盡量減小轉(zhuǎn)子葉片與機匣之間的接觸應力來提高其阻燃性能。

圖9 不同接觸應力摩擦加熱溫度歷史Fig. 9 Friction heating temperature history with different contact stresses

流速的變化會改變努塞爾數(shù)，進而對散熱量產(chǎn)生影響。如圖10所示，臨界著火溫度及著火延遲時間均隨著流速的提高而呈上升趨勢，當流速為100 m/s時，臨界著火溫度約為995 K，而著火延遲時間為2.3 s；當流速為200 m/s時，臨界著火溫度約為1000 K，而著火延遲時間為2.4 s；當流速為300 m/s時，此時流速接近音速，臨界著火溫度提高至1040 K左右，而著火延遲時間也提升至2.8 s；由此可見相比于跨音速階段，鈦火更易于發(fā)生在低速飛行階段（例如爬升階段），因此在高壓壓氣機防鈦火結(jié)構(gòu)設(shè)計時應考慮低速環(huán)境下的阻燃性能。

圖10 不同流速摩擦加熱溫度歷史Fig. 10 Friction heating temperature history with different flow velocities

臨界著火溫度及延遲時間隨著氧濃度ci的提高而不斷降低，如圖11所示。當ci為30 %時，臨界著火溫度約為950 K，著火延遲時間為1.9 s；當ci為42 %時，臨界著火溫度約為930 K，著火延遲時間為1.7 s；當ci為50 %時，臨界著火溫度為約920 K，著火延遲時間為1.5s。計算結(jié)果中臨界著火溫度的變化規(guī)律與文獻[4,11,24]中相一致，均呈指數(shù)下降趨勢，例如在文獻[24]的著火模型中，當ci分別為30 %、42 %、50 %時，臨界著火溫度分別為945 K 、896 K 、851 K，與本模型的計算值誤差在10%以內(nèi)，其誤差主要是解析模型在求根公式中采用泰勒展開所致?？梢?，環(huán)境中氧的比例（氧分壓）對著火參數(shù)影響比較顯著，通過減少環(huán)境中氧的比例是鈦火阻燃設(shè)計的重要途徑。從氧分壓較低的環(huán)境發(fā)生鈦火概率小的角度，需要考慮從地面飛行至高空爬升過程的高壓壓氣機的阻燃設(shè)計。

圖11 不同氧濃度下摩擦加熱溫度歷史Fig. 11 Friction heating temperature history with different oxygen concentrations

3 結(jié)論

（1）等溫加熱過程，隨著加熱面溫度急劇升高，著火延遲時間迅速降低，在不改變其他邊界條件的情況下，當加熱面溫度為1000 K時，臨界著火溫度約為958 K，著火延遲時間為2.3 s，當加熱面溫度為1941 K時，著火延遲時間縮短為0.2 s。

（2）非等溫線性加熱過程，若不考慮加熱速率對微凸體的影響，加熱速率為28 K/s、58 K/s及100 K/s的著火延遲時間分別為1.5 s、1.1 s和0.9 s，而臨界著火溫度基本維持在950 K；若考慮加熱速率對微凸體的影響，臨界著火溫度及延遲時間隨著微凸體尺寸的減小而下降，當微凸體直徑減小至16.5 μm時，臨界著火溫度約為765 K，與實驗值誤差為11.7%，著火延遲時間為1.1 s。

（3）非等溫摩擦加熱過程，著火延遲時間隨著接觸應力的增加而下降，當接觸應力達到26.5 kPa時，加熱速率提高至約130 K/s，著火延遲時間為1.4 s；隨著流速的增大而呈上升趨勢，當流速增大至300 m/s時，臨界著火溫度約為1040 K，著火延遲時間為2.8 s；隨著氧濃度增大而不斷降低，當氧濃度增大至50%時，臨界著火溫度約為920 K，著火延遲時間為1.5 s。