高速飛行器燃油瞬態溫度預測的修正熱網絡法

程顯達, 鄭皓冉, 楊學森, 董 威

(上海交通大學 機械與動力工程學院,上海 200240)

由于燃油具有比熱容大、性質穩定以及代償損失小等優點[1],現已被作為熱端部件的熱沉廣泛應用于現代高速飛行器設計中[2].但以燃油作為熱沉必須關注燃油的溫度變化,如溫度過高,不僅容易破壞其熱穩定性,而且形成的非溶性沉積物和樹脂沉淀物也可能堵塞油箱[3-4].因此,在儲油系統結構設計過程中研究油箱內燃油溫度分布尤為重要.試驗和數值模擬是獲得高速飛行器燃油溫度變化的兩個重要手段,雖然前者可獲得準確的燃油溫度數據,但因其無法在儲油系統結構設計階段完成,故基于數值模擬的燃油溫度變化預測對儲油系統結構設計起重要作用[5].

國內外對燃油溫度場計算的研究方法主要包括集總參數法、計算流體動力學(CFD)法和熱網絡法等.集總參數法假設燃油溫度一致,重點關注燃油熱平衡分析和邊界條件研究.在燃油熱平衡分析方面,Teare等[6]提出一種快速的油箱瞬態熱分析模型,該模型可以獨立分析油箱的結構傳熱和油箱內部的換熱問題;Bodie等[7]提出一種應用于包括燃油管路在內的飛機燃油系統整體網絡熱分析模型,基于該模型探究了燃油溫度的變化情況;Mawid等[8]介紹了一種航空飛行器燃油系統模擬工具Athena,并通過試驗驗證了其有效性;張興娟等[9]利用熱量平衡方法構建了燃油溫度計算模型,分析飛機超聲速巡航狀態下燃油溫度的瞬態變化.在邊界條件研究方面,Ai等[10]分析了影響油箱溫度場的內部導熱、外部對流和輻射等邊界條件,并在此基礎上建立了燃油瞬態計算模型.集總參數法簡單高效,但忽略了燃油溫度在空間上的差異,無法得到局部區域上的燃油溫度,尤其是無法關注到靠近油箱壁面位置的高溫區域油溫[11],因此這一方法難以滿足高速飛行器油箱設計的需要.

使用CFD方法可對油箱整體結構建模和計算,從而獲得油箱內溫度場分布情況.陳悅[4]采用二維和三維的方法模擬了輸油箱和供油箱內燃油流動換熱過程,獲得了油箱內燃油的溫度分布特點及燃油出口溫度隨時間變化而變化的過程;趙璇[11]采用CFD方法分析了超聲速飛行器機翼油箱內外流場溫度分布,并將CFD方法與集總參數法結果進行對比.CFD方法可以精確捕獲油箱溫度場分布,但對飛行全過程進行模擬所需的計算開銷極大,難以滿足油箱設計階段快速迭代優化的需要.

為解決集總參數法得不到局部溫度和CFD方法計算效率過低的問題,熱網絡法被應用于溫度場計算,尤其是在固體溫度場計算中,如軸承[12]、電動機及發電機[13]、管路[14]和電路[15-16]等.康振燁等[17]將熱網絡法應用于燃油溫度場計算,將不同位置的燃油劃分到不同的基元模塊中,采用MATLAB/Simulink軟件搭建飛機油箱穩態熱分析仿真模型;呂亞國等[5]在此基礎上建立了飛機油箱非穩態熱分析模型,結果表明熱網絡法的計算結果與飛行測試數據吻合良好.上述研究假設油箱基元模塊節點之間的傳熱方式為固體導熱,這一假設僅適用于燃油溫度較低情況,對于具有高油溫的高速飛行器儲油系統,燃油內部會存在較強的冷熱摻混,假定燃油基元模塊間傳熱方式為固體導熱的方法會帶來較大誤差.不僅如此,目前的熱網格法未考慮燃油內部的宏觀流動,隨著燃油消耗,油箱中燃油存在整體向下的宏觀流動,對于每一個位置固定的基元模塊,都需要考慮燃油宏觀流動帶來的傳質影響.

針對高速飛行器燃油瞬態溫度預測問題,本文基于文獻[5]和文獻[17]提出一種滿足設計精度的高效修正熱網絡方法,通過修正熱傳導方程的方法考慮油箱內部冷熱摻混的影響,通過修正熱平衡方程的方法考慮燃油宏觀流動的影響.基于該熱網絡法構建了高速飛行器儲油系統瞬態熱分析模型,在不同的熱邊界條件和儲油系統幾何結構下校核其普適性.

1 高速飛行器的油箱模型

典型高速飛行器油箱結構的三維模型如圖1所示.油箱內側面裝有隔熱材料以防止燃油溫度過高;油箱左上和左下側貼近飛機外壁面,飛行過程中會受到高速氣流的氣動加熱作用;油箱右下側貼近發動機進氣道,會受到發動機的熱傳導加熱;其余3個面以及兩端端面與其他部件相連,可近似為絕熱條件.在飛行過程中,通過高壓氮氣擠壓氣囊為發動機提供帶有壓力的燃油[18].

圖1 油箱示意圖Fig.1 Diagram of fuel tank

計算中高速飛行器油箱初始條件和模型簡化如下:① 油箱在垂直于紙面方向上結構相似,且端面為絕熱條件,因此將油箱模型簡化為二維模型處理,油箱最上端和最下端表面為燃油入口和出口,如圖2所示;② 初始時刻油箱內充滿燃油,燃油和隔熱材料溫度均為15 ℃,1 200 s內燃油均勻消耗完;③ 外部的氣動加熱效果等效為對流加熱邊界條件,其對流換熱系數和來流溫度隨時間的變化如圖3所示,t為時間,T為溫度,α為對流傳熱系數;④ 將內部發動機進氣道的熱傳導作用等效為溫度邊界條件,溫度隨時間的變化與外部相同;⑤ 忽略氣囊、蒙皮等結構對油箱溫度場的影響.

圖2 油箱計算區域示意圖Fig.2 Diagram of fuel tank calculation domain

圖3 隨時間變化的邊界溫度和對流傳熱系數Fig.3 Variation of boundary temperature and heat transfer coefficient with time

2 修正熱網絡方法

2.1 熱網絡基元模塊分類

對基元模塊進行分類是因為不同的基元模塊存在不同的熱平衡方程.油箱內存在燃油和氮氣兩相,因氮氣并非熱設計的重點關注對象,在一定范圍內可以接受氮氣升到比較高的溫度,且其密度較小,熱擴散系數比較大,溫度更易趨于一致,故采用集總參數法的方式處理,即將油箱中填充氮氣的區域作為一個完整的基元模塊,定義為氮氣基元模塊.

對于燃油則需要對其內部區域進行劃分,以其形狀規則的油箱結構為例闡述其區域劃分,如圖4所示.首先,將油箱腔體等距劃分為4×4個區域,將每個區域內的燃油作為一個單獨的基元模塊,認為每一個基元模塊內部溫度相同.初始時刻,油箱內充滿燃油,隨著燃油消耗,油氣界面不斷下降,處于最上層的燃油基元模塊體積不斷縮小,其他燃油基元模塊體積恒定.根據燃油消耗過程中某一時刻基元模塊體積是否發生改變,將燃油基元模塊定義為兩種類型:體積發生改變的模塊定義為半油基元模塊;體積恒定不變的模塊定義為滿油基元模塊.隨著燃油持續消耗,當油氣界面到達基元模塊交界面位置時,最上層的燃油基元模塊體積變為0后消失,原本處于第二層的滿油基元模塊轉變為半油基元模塊,體積不斷縮小,處于下層的基元模塊仍為滿油基元模塊,以此類推直至燃油全部耗盡.

圖4 半油基元模塊和滿油基元模塊隨燃油消耗變化示意圖Fig.4 Variation of partial fuel module and full fuel module with fuel consumption

除了對油箱內部氮氣和燃油區域劃分基元模塊外,對隔熱材料區域也可進行劃分,將其定義為隔熱材料基元模塊.

2.2 基元模塊熱平衡方程

2.2.1半油基元模塊 根據基元模塊內能的變化量等于與周圍基元模塊的傳熱量和傳質量,可將半油基元模塊和滿油基元模塊的熱平衡方程表述為

(1)

式中:Ti為基元模塊i的瞬時溫度;ρcVi為基元模塊i的密度、比熱容和體積的乘積;Qi為基元模塊i與周圍基元模塊之間的傳熱量;Ei為由于燃油整體向下流動導致的基元模塊i與周圍基元模塊之間的傳質量.

假定燃油密度和比熱容恒定,半油基元模塊體積和溫度均隨時間變化而變化,式(1)左邊展開為

(2)

式中:wi為基元模塊的相界面面積與總的相界面面積之比;Vfuel為燃油總體積.

假設燃油均勻平穩流出,燃油流線方向與縱向的熱網絡劃分線平行.對于燃油基元模塊i,僅在縱向上存在傳質作用,對于半油基元模塊僅存在流出項,傳質項可以寫為

(3)

式中:dmi/dt為流出基元模塊i的質量流量,與基元模塊i質量的變化量ρdVi/dt大小相等,符號相反.

燃油體積的變化率采用下式計算:

(4)

(5)

式中:Vfuel,0為初始時刻燃油體積.對于基元模塊i,其體積變化速率dVi/dt為總體積變化速率dVfuel/dt的wi倍,其體積大小也可通過式(5)的方式積分得到.

將式(2)～(5)帶入式(1)可得半油基元模塊的熱平衡方程為

(6)

2.2.2滿油基元模塊 滿油基元模塊熱平衡方程的基本形式同式(1),不同的是其體積不隨時間發生改變,其傳質項除了考慮流出基元模塊的質量流量,還應考慮從上方流入的基元模塊質量流量:

(7)

式中:dmin/dt為流入基元模塊i的質量流量,因滿油基元模塊質量守恒,其值與流出的質量流量dmi/dt相等;Tin為流入基元模塊i質量流量的溫度.

最后整理可得滿油基元模塊的熱平衡方程為

(8)

2.2.3氮氣基元模塊 氮氣基元模塊內的質量可通過理想氣體狀態方程表述為

(9)

式中:p為壓力,保持1.5個標準大氣壓恒定;Rg為氮氣的氣體常數;氮氣體積VN2為油箱總體積與燃油體積之差.對式(9)兩邊取微分得到

(10)

式中:TN2為氮氣溫度;氮氣體積變化率dVN2/dt與燃油體積變化率dVfuel/dt大小相等,符號相反.

氮氣的熱平衡方程為

(11)

式中:cV為定容比熱容;cP為定壓比熱容;TN2,in為進入油箱的氮氣溫度.進一步整理可得

(12)

2.2.4隔熱材料基元模塊 對于油箱隔熱材料等固體區域,也可以劃分熱網絡,其控制方程為

(13)

式中:下標s表示固體區域.

2.3 油箱基元模塊劃分

基于燃油的流動方向可將油箱結構離散化,劃分為若干個基元模塊,熱網絡劃分方式如圖5所示.其中,因隔熱材料在平行于邊界面方向上溫差較小而在邊界面法線方向上存在較大溫差,故將每個隔熱材料在沿邊界法向方向上劃分為3個等寬度的基元模塊,編號為①～⑨;然后,對燃油角落區域的油箱結構進行熱網絡劃分,角落區域直接受到熱載荷,其溫度與油箱中央區域相差較大,需單獨處理,編號為⑩～,同時角落區域控制體流動情況較為復雜,故在計算過程中進行適當簡化,在控制方程(1)中不考慮傳質項的影響;最后,對燃油中央區域進行基元模塊劃分.從CFD計算結果來看,中央區域溫差不大,基元模塊不必劃分過多.中央區域基元模塊橫向劃分線為水平線,縱向劃分線為穿過橫向線的1/4和3/4等分線,即在中央區域任畫一條水平線,縱向劃分線均穿過其1/4和3/4等分點,編號為～.忽略由溫度梯度引起的燃油自然對流,只考慮燃油消耗引起的燃油宏觀流動,則燃油流線與縱向劃分線趨于平行,因此基元模塊之間只需考慮縱向傳質.定義氮氣區域為編號的基元模塊.

3.3 主題活動“人類的起源與發展” 在開展本次活動前，學生需要具有的前概念是： 人類在自然界中的位置、人類起源過程中新舊特征的更替、現代進化理論和人類進化的歷程等。此外，教師在活動前出示南方古猿、能人、直立人和智人頭骨圖及相關簡介。

圖5 油箱基元模塊Fig.5 Cell module of fuel tank

2.4 基元模塊間傳熱量計算

基元模塊i與周圍基元模塊傳熱量計算公式為

(14)

式中:Tj為相鄰基元模塊j的溫度;Rij為基元模塊i和基元模塊j之間的熱阻.傳熱量計算的關鍵在于兩個基元模塊之間熱阻的計算.隔熱材料內部熱阻計算公式為

R=L/(λA)

(15)

式中:L為兩基元模塊節點之間的距離;λ為材料的導熱系數;A為導熱面積.

因燃油受熱之后內部存在冷熱摻混,燃油內部之間的熱阻若直接按照固體導熱熱阻計算,內部換熱將遠高于被看作固體導熱時的換熱強度,但燃油內部的冷熱摻混現象較為復雜,難以直接求解,因此提出摻混系數的概念.將燃油節點之間的導熱系數λ調整為kλ(k>1),摻混系數k表征由于燃油內部冷熱摻混導致的換熱量增加,燃油基元模塊之間的熱阻計算公式表示為

R=L/(kλA)

(16)

燃油與壁面、氮氣與壁面熱阻計算公式為

R=1/(αA)

(17)

式中:α可以采用經驗公式的方法求解.對于形狀復雜的壁面,也可采用CFD計算方法提取或修正經驗公式.

為了將對流傳熱系數的經驗公式用于油箱壁面,將油箱壁面等效為3種類型,如圖6所示.

圖6 等效的油箱壁面示意圖Fig.6 Diagram of equivalent fuel tank wall

將圖6中類型1等效為熱面朝下的水平壁面,其對流傳熱系數計算公式[19]為

(18)

式中:D為特征長度;Gr為格拉斯霍夫;Pr為普朗特數.

將圖6中類型2等效為熱面朝上的水平壁面,其對流傳熱系數計算公式[20]為

(19)

將圖6中類型3等效為與豎直方向夾角為θ的傾斜平面,其對流傳熱系數計算公式[20]為

(20)

對于燃油與氮氣控制體之間的熱阻計算,因氮氣熱阻遠高于燃油熱阻,故計算過程中僅考慮氮氣對流熱阻而忽略燃油熱阻,對流傳熱系數由式(19)計算.

2.5 方程求解方法

油箱溫度場計算流程如圖7所示.

圖7 油箱溫度場計算流程圖Fig.7 Flow chart of fuel tank temperature field calculation

在初始溫度場已知情況下,先對各個基元模塊的類型加以判定,包括半油基元模塊、滿油基元模塊、氮氣基元模塊和隔熱材料基元模塊,然后求解下一時刻基元模塊體積和體積變化率,液面高度H可通過對下式迭代求解得到:

(21)

式中:A(h)為油箱內高度為h的位置的橫截面面積,在幾何模型確定的情況下為已知量.如液面高度經過基元模塊交界面,則意味著一些半油模塊的消失和一些滿油基元模塊轉變為半油基元模塊,此時需要更新燃油區域的基元模塊類型;隨后,計算各個基元模塊之間的傳熱量;最終,可以得到每個控制體對應的熱平衡方程.各個控制體的熱平衡方程組成一個常微分方程組:

(22)

式中:i=1, 2, …, 29.方程組的初始條件為Ti=288.15K.采用龍格-庫塔法求解式(22),最終得到各個基元模塊隨時間的溫度變化.

3 結果與討論

3.1 燃油溫度的CFD計算結果

CFD的結果可以為網格法提供合理的修正系數及方法驗證.采用 ANSYS Fluent軟件進行高速飛行器燃油溫度變化的計算,流場控制方程采用雷諾平均Navier-Stokes方程,湍流模型采用標準k-epsilon湍流模型計算.CFD作為網格法驗證的基礎,本文首先利用規則形狀油箱的實驗對采用的CFD計算方法進行驗證.

實驗中油箱外形尺寸為500 mm×300 mm×500 mm,由側板、蓋板、剛性隔熱墊和隔熱氈等部件組成,其外形如圖8所示,實驗加熱系統如圖9所示.實驗中加熱表面溫度保持400 ℃,對比分析上表面加熱滿油、側表面加熱滿油、側表面加熱半油3個實驗狀態下的溫度測量結果和CFD計算的結果,如圖10～圖12所示(dbo為距底部距離),得出CFD計算結果與實驗結果的趨勢一致,溫度計算值與實驗測量值吻合較好.

圖8 實驗油箱外形示意圖Fig.8 Schematic diagram of experimental fuel tank

圖9 實驗加熱系統Fig.9 Experimental heating system

圖10 上表面加熱滿油狀態油溫分布對比Fig.10 Comparison of fuel temperature distribution when heated on upper surface with full fuel

圖11 側表面加熱滿油狀態油溫分布對比Fig.11 Comparison of fuel temperature distribution when heated on side surface with full fuel

圖12 側表面加熱半油狀態油溫分布對比Fig.12 Comparison of fuel temperature distribution when heated on side surface with partial fuel

在CFD方法驗證基礎上,進行實際油箱CFD數值計算,計算網格如圖13所示.對邊界層區域網格進行加密處理,網格總數為 17 048 個,滿足網格無關性要求.

圖13 油箱計算網格Fig.13 Grid of fuel tank

選取RP-3燃油為研究對象,密度采用Boussinesq假設處理,參考密度為700 kg/m3,熱膨脹系數為 0.003 128 2 K-1,氮氣密度按照理想氣體方程計算.隔熱材料熱導率為0.05 W/(m·K),密度為128 kg/m3,比熱容為929.7 J/(kg·K),入口壓力為p,溫度為15 ℃,出口采用速度邊界條件以表示燃油勻速率流出.計算中氮氣和燃油界面采用流體體積(VOF)多相流模型.

1 200 s 內油箱內氮氣和燃油速度(v)分布云圖和流線圖如圖14所示.從速度云圖中可見,由于燃油和氮氣受熱密度發生變化,二者均產生了明顯的自然對流.氮氣對流速度約為0.1 m/s,對流較為強烈;燃油內部對流較弱,部分區域特別是底部流動較小,上端對流速度約在0.01 m/s量級.

圖14 不同時刻油箱內速度和相界面分布云圖Fig.14 Velocity and phase contours at different times

1 200 s 內油箱內溫度分布云圖如圖15所示,3個隔熱材料在邊界面法線方向上從邊界上的834 K快速過渡到內部的400～500 K,而在平行于邊界面方向上溫度接近.因氮氣具有較大的熱擴散系數,邊界的熱量更容易傳輸到內部,而燃油熱擴散系數較小,熱量不容易傳入,故總體上氮氣溫度高于燃油.

圖15 不同時刻油箱內溫度分布云圖Fig.15 Temperature contours at different times

3.2 與CFD的結果對比及摻混系數影響分析

取中心點位于圖5中點1的典型絕熱材料基元模塊③和中心點位于圖5中點2的燃油基元模塊進行結果分析:點1靠近燃油邊界區域,可反映油氣交界面區域的溫度情況;點2位于整個油箱的中心位置,可反映燃油中央區域的溫度情況.

典型絕熱材料控制體和典型燃油控制體熱網絡法以及CFD計算結果如圖16和圖17所示.當摻混系數k=1時,將燃油控制體之間的傳熱方式按照固體傳熱計算,熱網絡法得到的隔熱材料溫度與CFD計算結果相近,但燃油溫度顯著偏低,這是由于摻混系數k=1時,忽略了燃油對流換熱導致的內部換熱量增加,從而導致誤差較大.

圖16 絕熱材料基元模塊CFD和熱網絡法計算結果Fig.16 Calculation results of thermal network method and CFD method for insulator module

圖17 燃油基元模塊CFD和熱網絡法計算結果Fig.17 Calculation results of thermal network method and CFD method for fuel module

基于CFD計算結果對摻混系數k進行修正,定義:

(24)

式中:tm=100, 200, …, 1 200 s,Ti,CFD為t時刻基元模塊i中心區域溫度的CFD計算結果;Ti,k為t時刻通過熱網絡方法得到的摻混系數為k時基元模塊i的溫度計算結果.最終的摻混系數k取式(24)中的最小值點,本文中k=5 000,表明燃油內部存在較為強烈的對流換熱,此時CFD計算結果與熱網絡法計算結果相對誤差小于5%,說明引入摻混系數k的熱網絡法可以有效對燃油溫度場進行計算.

3.3 熱網絡法對邊界條件普適性驗證

為驗證熱網絡法對不同邊界條件的適用性,保持摻混系數k=5 000 不變,取其他可能達到的最高和最低熱邊界(穩定狀態溫度增減200 K)條件重新計算,并與CFD結果對比,結果如圖18和圖19所示.

圖18 不同邊界條件下隔熱材料基元模塊計算結果Fig.18 Calculation results for insulator module under different boundary conditions

圖19 不同邊界條件下燃油基元模塊計算結果Fig.19 Calculation results for fuel module under different boundary conditions

由圖可見,當邊界條件溫度改變后,熱網絡法仍可準確獲得溫度計算結果,與CFD結果相比相對誤差在5%以內,驗證了熱網絡法在一定邊界條件范圍內的普適性.如需計算溫度偏離較大的情況,則可以改變邊界溫度重新進行CFD計算獲得對應溫度下的摻混系數,根據兩點或多點的摻混系數值,構造摻混系數關于溫度的函數,提高溫度適用范圍.

3.4 熱網絡法對幾何條件普適性驗證

油箱設計過程中,其外形結構受油箱外其他設備限制不會發生很大改變,油箱結構優化設計通常是對隔熱材料厚度的微小調整.

為了驗證熱網絡法對不同隔熱材料厚度的適用性,保持摻混系數k=5 000 不變,將原本偏厚的隔熱材料厚度減少一半重新計算,并與CFD結果對比,計算結果如圖20和圖21所示.

圖20 不同隔熱材料厚度下隔熱材料基元模塊計算結果Fig.20 Calculation results for insulator module with different thicknesses of insulation material

由圖可見,隔熱材料厚度減少一半后,熱網絡法仍可準確獲得溫度計算結果,與CFD結果相比相對誤差在4%以內,同樣驗證了熱網絡法對隔熱材料厚度的普適性.如需計算幾何結構調整較大的情況,可以改變幾何結構重新進行CFD計算獲得對應幾何結構下的摻混系數,根據兩點或多點的摻混系數值,構造摻混系數關于幾何結構的函數,提高溫度適用范圍.

此外,CFD計算結果表明,當隔熱材料厚度減少后,中央區域燃油溫度有所下降,這一結果也被熱網絡法有效捕捉,這是由于在當前物性和邊界條件下,隔熱材料厚度減少,燃油總量提升,熱沉作用更加明顯;另一方面,厚度減少后,隔熱材料與燃油的對流換熱面積也相應減少,導致燃油溫度降低.但隔熱材料厚度過小會導致隔熱材料和燃油交界面處的燃油溫度提升,在油箱熱設計過程中需要對上述因素進行權衡考慮.

最后,以CPU主頻3.70 GHz、10核運行的計算機比較CFD方法和熱網絡方法的計算時間.CFD程序單次計算時間約為30～50 h,而使用相同計算機運行單次熱網絡法程序的計算時間在1 min以內,在計算速度上熱網絡方法遠超CFD方法.從油箱熱設計的全周期考慮,在本方法提出前,需要數十次CFD計算獲得優化后的幾何結構,本方法則僅需少量CFD計算獲得修正系數,將優化迭代過程依托于修正熱網絡方法,從而極大縮短研發周期.

4 結語

針對高速飛行器飛行全階段燃料溫度瞬態預測問題,基于修正的熱網絡法構建的燃油瞬態熱分析模型具有與CFD接近的計算精度,可準確計算出油箱內燃油溫度場的分布和變化情況,對不同邊界條件、幾何條件均有較好普適性,計算效率遠高于CFD方法.在實際油箱結構優化過程中,可以通過修正熱網絡法對油箱結構進行快速迭代設計,尋找滿足條件的最優方案.

修正的熱網絡法中摻混系數k表征了燃油受熱之后內部冷熱摻混的劇烈程度,其主要影響因素是燃油物理性質、邊界溫度和幾何條件等.通常情況下,航空燃油種類固定,燃油物性變化不大,摻混系數k通過一次CFD校正后,在油箱設計的迭代優化過程中無需再次調整.通過將摻混系數調整為與之相關性較強的變量函數,可以將本方法進一步拓展至其他復雜邊界條件和幾何條件的適用范圍,應用于其他形式的高速飛行器油箱溫度場計算.