魚雷燃燒室氣-液-固三相耦合換熱數值仿真

陳 肯,馮要飛,宗 瀟,伊進寶,肖炎彬,史小鋒

(中國船舶集團有限公司 第705 研究所,陜西 西安,710077)

0 引言

燃燒室是魚雷熱動力裝置中最重要的部件之一[1],推進劑在燃燒室內燃燒,將化學能轉變為熱能,并推動發動機做功[2]。為保證燃燒室殼體在高溫高壓的工作環境下具有足夠的熱強度,需通過冷卻水對燃燒室進行冷卻換熱[3]。在冷卻過程中,燃燒室因其結構復雜會出現局部區域過熱導致殼體強度失效的問題。為指導燃燒室的結構改進,需要對燃燒室的冷卻換熱過程進行研究。劉學等[4]利用數值仿真計算的方式,研究了燃燒室多孔介質幾何結構對其冷卻效果的影響;Chen 等[5]基于雙向熱流固耦合的方法,采用大渦模擬(large eddy simulation,LES)對非渦流環形燃燒室中的自激旋轉模式方位不穩定性進行了研究,并在此基礎上構建了頻率預測匹配模型;Laurent 等[6]建立了液體火箭發動機的湍流噴霧模型,并根據建立的模型對發動機進行了數值分析;Zhukov 等[7]采用了擴展的渦耗散模型,用CFX 軟件對多孔噴嘴頭的火箭燃燒室進行了仿真。上述研究聚焦于陸用及航空航天用燃燒室的冷卻換熱過程,對水下渦輪發動機燃燒室的冷卻換熱過程不適用。

趙衛兵等[8-9]利用理論解析的方式計算了魚雷燃燒室在低頻不穩定燃燒狀態下的溫度特性,并針對三組元推進劑的燃燒室流場進行了基于Fluent 的數值研究。Qin 等[10]采用雙向熱流固耦合的方式,比較了部分進氣式渦輪機燃燒室的質量流量和效率的試驗與仿真結果,分析了5 kW 功率下渦輪機的輸出特性,并建立了渦輪機輸出參數與輸入邊界條件間的一維設計計算模型。上述工作對水下燃燒室的工作狀況及冷卻特性進行了數值計算研究,但上述研究并未同時涉及氣-液-固三相的耦合換熱過程。

文中以燃燒室2 種不同冷卻水流道結構作為仿真對象,利用Fluent 軟件首次開展燃燒室內高氯酸羥胺(HAP)三組元推進劑燃燒,燃氣、冷卻水流動換熱,以及殼體參與導熱的氣-液-固三相傳熱換熱過程的熱流固耦合數值仿真,分析相同工況及冷卻水流量下燃燒室不同結構的冷卻換熱特性差異及造成局部區域過熱的原因,為進一步研究和改進燃燒室的結構提供理論參考。

1 冷卻換熱理論分析

燃燒室的冷卻換熱過程涉及到的物理模型包括: 燃料的霧化、燃燒模型,燃氣及冷卻水的流動模型,燃氣、冷卻水的對流換熱模型,固體域的導熱模型,以及燃氣、固體域的輻射傳熱模型。

1.1 流體控制方程

采用N-S 方程描述流體域的內流場。

連續性方程

式中: ρ為流體密度;t為時間;u為流體速度;Sm為流體控制體在單位時間質量因源項產生的變化量。

動量方程

式中:P為流體控制體表面的壓力;μ為流體的動力粘度;g為重力加速度;F為作用于流體控制體的質量力合力。

能量方程

式中:e為流體控制體的能量;κ為導熱系數;Q為控制體內單位時間熱源項提供的熱量。

式中:x代表一維換熱過程中的位置變化;i、j指代節點坐標。

文中,重力對仿真計算的影響較小,故不考慮其影響;冷卻水域無源項,故Sm、F、Q均等于0。在燃燒域中,Sm為注入的燃料在單位時間蒸發的流量;F為燃氣連續相受到的除重力外的體積力,F=0;Q為液態燃料蒸發吸收的熱量以及燃燒產生的化學能之和。

1.2 燃料霧化及燃燒控制方程

1.2.1 燃料霧化

液滴直徑取中間質量的定義方式,根據液滴霧化的累計分布R-R 函數,可得液滴直徑在R-R 函數分布情況下,按中間質量方式定義的直徑表達式為[2]

式中:dm為中間質量直徑;d*為液滴破碎后的標稱直徑;n為破碎系數。

1.2.2 液滴運動及蒸發

液滴經破碎霧化進入燃燒室內的運動控制方程為[11]

式中:up為顆粒速度;ρp為顆粒密度;FX為其他作用力;dp為顆粒直徑;Re為顆粒雷諾數;CD為阻力系數。

未達到蒸發溫度的溫度控制方程為[12]

達到蒸發溫度后的液滴質量控制方程為

式中:mp,cp,Tp,εp,Ap分別為顆粒的質量、比熱、溫度、發射率和表面積;h為對流傳熱系數;T∞為氣相溫度;σ為波爾茲曼常數;θR為輻射溫度;hfg為潛熱。

1.2.3 氣體燃燒

組分輸運模型可以模擬燃料蒸發產生氣體的燃燒過程。組分輸運體積反應模型渦耗散(eddydissipation)模型的反應速率Ri,k的控制方程為[6]

式中:mq、mR分別為所有產物組分和某一特定反應物的質量比重;Mw,i、Mw,R分別為第i、R種物質的分子量;分別是反應物和生成物的化學當量系數;A和B是經驗常數,分別取4 和0.5。

HAP 三組元推進劑完全燃燒的化學反應方程式為

2 模型前處理及網格劃分

2.1 計算域簡化及提取

進行燃燒室熱流固耦合仿真計算,對燃燒室的幾何結構進行適當簡化,可以降低網格劃分的難度,降低網格數量,提高網格質量,幫助仿真模型更快更好地收斂。

利用提取功能處理簡化完成的燃燒室固體域,得到研究對象的全體計算域,如圖1 所示。

圖1 計算域提取圖Fig.1 Extraction graph of computational domain

燃燒室在工作過程中,固體之間的接觸部位會產生縫隙,燃氣通過縫隙泄漏進入內襯與殼體之間的保溫層,擾動保溫層內的氣體流動,保溫層區域由死腔變為流動區域,主要換熱形式由自然對流變為強制對流。將固體域變形產生的縫隙修正到提取的計算域模型中,保溫層域與燃燒域通過0.25 mm 的縫隙連通。

文中燃燒室采用2 種不同冷卻水流道結構,差異主要集中在冷卻水域的碗狀底部以及最后一段的冷卻水出口流道上,如圖2 所示。不同結構參數對比見表1。

表1 不同結構參數對比Table 1 Parameters comparison for different structures

圖2 不同結構冷卻水域流道差異Fig.2 Channel differences of the cooling water area for different structures

2.2 計算域網格劃分

采用ICEM 軟件完成冷卻水域、燃燒域與保溫層域的網格劃分。對于冷卻水域的螺旋流道部分,使用結構化網格生成其中1 條流道的網格,通過周期陣列得到8 條螺旋流道的結構化網格,對冷卻水域的碗狀底部、出口流道及燃燒域全域,利用八叉樹算法生成四面體非結構化網格。

采用Mesh 模塊完成內襯與殼體域的網格劃分。采用Mechanical 算法生成內襯域及殼體域的四面體非結構化網格。

設置接觸面網格尺寸,進行網格的局部加密,保證不同計算域在交界面處的網格單元節點個數與單元尺寸大小大致相等。將劃分后的網格導入Fluent 中,利用Make Polyhedral 功能將四面體網格轉為多面體網格,顯著降低網格單元與節點數量,提高仿真計算的收斂精度與速度。網格劃分結果如圖3 所示(其中,圖(b)為圖(a)圈內放大部分)。

圖3 計算域網格劃分結果Fig.3 Grid division result of computational domain

2.3 網格無關性驗證

通過改變網格的大小得到不同的網格劃分方案,進行網格無關性驗證,最終用于網格無關性驗證的4 種網格劃分方案見表2,網格無關性驗證結果見圖4。

表2 網格無關性驗證網格劃分方案Table 2 Grid division scheme of grid independence verification

圖4 網格無關性驗證結果Fig.4 Verification results of grid independence

圖4 以計算域冷卻水域出口處溫度分布作為衡量網格數量對計算結果精度影響的指標。從圖中可以看出,冷卻水域的可接受網格尺寸為0.4 mm,燃氣域的可接受網格尺寸為3 mm,在可接受的網格尺寸基礎上,繼續加密網格的數量不會對仿真計算的結構造成影響,方案2 的網格劃分可滿足對計算精度的要求。

3 仿真模型選擇及驗證

3.1 湍流模型選擇

流體機械中雷諾數大于105認為處于高雷諾數湍流狀態,針對湍流流動狀態,仿真模型適合選用兩方程模型中的k-ε模型進行流體運動的仿真。k-ε模型包括標準k-ε、RNG(renormailization group)k-ε和可實現k-ε等3 種模型。

以結構1 為仿真的計算域,采用網格劃分方案2 的方式劃分網格。表3 為進行湍流模型驗證時設置的計算域邊界條件。

表3 湍流模型選擇邊界條件設置Table 3 The boundary condition setting of turbulence model selection

監測3 種不同湍流模型的燃氣域出口溫度和冷卻水域出口溫度,計算冷卻水溫升,并分別與同工況下的試驗結果對比,求出燃氣出口溫度與冷卻水溫升的相對誤差,如表4 所示。

表4 不同湍流模型仿真與結果Table 4 Simulation and experiment results of different turbulence models

從表中可知,對于燃燒域的出口溫度、冷卻水域的出口溫升,可實現k-ε模型的計算精度最高,因此湍流模型選擇可實現k-ε模型。

3.2 其他仿真模型選擇

燃燒室實際工作過程中,燃料液滴在很短的時間內燃燒完全。文中選用組分輸運模型作為燃燒模型,該模型的計算原理為氣相反應物接觸后瞬間完成體積反應,能很好地適用燃燒室的工作狀況。

燃燒室的燃燒過程中注入的液滴顆粒會參與到輻射換熱過程。由于數據/領域對象(data/domain object,DO)模型適用于所有光學深度的輻射過程及液滴顆粒燃燒過程的輻射模擬,故被選用進行計算。

3.3 仿真模型驗證

以燃燒室結構1 為驗證仿真模型準確性的計算域,以其在3 種不同工況下的工作環境作為仿真邊界條件,進行熱流固耦合仿真計算。因冷卻水的溫度變化數值能直觀反映整個系統的冷卻換熱功率,故選用冷卻水溫升作為衡量模型準確性的指標。分別對比3 種工況下的仿真結果與試驗在結果的誤差,以驗證模型的準確性。

仿真結果、試驗結果以及計算的模型誤差大小見表5。

表5 不同工況下仿真與試驗結果對比Table 5 Results comparation between the simulation and experiment on the different working conditions

從表中可以看出3 種不同工況下,仿真模型的誤差大小均控制在10%以內,建立的仿真模型準確。

4 燃燒室不同冷卻水流道結構冷卻換熱仿真分析

燃燒室2 種冷卻水流道結構冷卻工況對比時,除結構差異外其余邊界條件參數設置見表6。

表6 除結構差異外其余邊界條件設置Table 6 The other boundary condition setting except the structural difference

表7 給出了表6 邊界條件下針對不同結構的仿真結果。可知,結構2 相較于結構1 冷卻水的出口溫度更高,燃氣溫度更低,冷卻水溫升更高,殼體的內表面均溫及最高溫度均更低。結構2 的換熱功率為101.3 kW,結構1 的換熱功率為93.3 kW。結構2 相較于結構1 換熱功率提高了8.5%。

表7 不同結構仿真結果對比Table 7 The comparison of simulation results for different structures

探究結構2 冷卻效果提升的原因,具體分析2 種結構溫度場與流場的差異。

4.1 燃燒域冷卻特性

圖5 是燃燒室2 種冷卻水流道結構的燃燒域溫度分布云圖,可知2 種結構的燃燒域溫度分布幾乎一致。2 種結構燃燒域出口段的溫度隨位置的變化曲線如圖6 所示,燃氣溫度隨位置的變化規律大致相同,溫度在數值上的最大差距為10 K,相較于燃燒產生的燃氣1 400 K 左右的溫度,10 K 溫差對于傳熱的影響可以忽略。

圖5 不同結構燃燒域溫度分布云圖Fig.5 Temperature contours of combustion area for different structures

圖6 不同結構燃燒域出口端溫度分布曲線Fig.6 Temperature curves of combustion area in the outlet for different structures

圖7 是燃燒室不同結構燃燒域的速度流線圖及在中截面處的速度矢量圖,可知2 種結構的流場狀態幾乎相同,燃氣最大流速均約61 m/s,且均在燃燒域下腔與出口流道交匯處達到最大。燃燒域的流場狀況相似,可以佐證溫度場的分布狀況相似。

圖7 不同結構燃燒域速度流線圖及矢量圖Fig.7 Velocity streamlines and vectors of combustion area for different structures

4.2 冷卻水域冷卻特性

圖8 是燃燒室2 種結構冷卻水域溫度分布云圖,可知冷卻水在螺旋流道處的溫度分布狀況基本一致。圖9 為燃燒室2 種結構在螺旋流道段的溫度分布圖,螺旋流道終端溫度均約295 K,且溫度隨位置的變化規律大致相同。

圖8 不同結構冷卻水域溫度分布云圖Fig.8 Temperature contours of the cooling water area for different structures

圖9 不同結構冷卻水域螺旋流道溫度分布曲線Fig.9 Cooling water temperature curves in the spiral channel for different structures

由圖8 可知,2 種結構的碗狀底部及出口流道段,冷卻水的溫度分布存在較大差別。結構1 碗狀底部的冷卻水平均溫度小于結構2 的溫度。結構1 碗狀底部的6 條流道中,5 條流道的左側存在一處大面積的高溫區域,右側則存在一處較為細長的條狀低溫區域,在碗狀流道的最低端交匯后迅速消失。結構2 碗狀底部的18 條流道中,同樣存在部分流道內形成局部的高溫區域。高溫區域的大小相較于舊結構更小,溫度更高,分布更為隨機。但整體上分布于流道上半部分較寬的區域。與出口流道呈90°夾角的2 條流道的整體溫度相較于其他流道更高。2 種結構的冷卻水域出口流道的溫度分布整體上呈現越靠近冷卻水域出口處溫度越高的趨勢,這是冷卻水流動換熱的原因,符合傳熱學規律。結構1 相較于結構2 在出口流道段會形成許多局部高溫或者低溫區域。

圖10 是監測的2 種結構冷卻水域出口流道一段長度上的溫度分布隨位置變化規律曲線,可以看出,結構2 冷卻水溫度越靠近出口溫度越高,而結構1 的冷卻水溫度則在靠近出口處時出現了2 個局部高溫。

圖10 不同結構冷卻水域出口段溫度分布曲線Fig.10 Temperature curves of the cooling water in the outlet area for different structures

圖11 是燃燒室2 種結構在冷卻水域碗狀底部差異處的溫度云圖。表8 為監測所得碗狀底部流道不同結構差異處截面溫度均值。

表8 不同結構差異處溫度值Table 8 The temperatures at the differences for different structures

圖11 不同結構差異處冷卻水域溫度分布云圖Fig.11 Temperature contours of the cooling water area at the differences for different structures

從圖11 可知,2 種結構的冷卻水在流入碗狀底部前,在碗狀底部入口處的溫度基本相等;在流出碗狀底部進入出口流道入口時,結構2 的冷卻水溫度高于結構1。從表8 中可知,結構2 在碗狀流道段與出口流道段的換熱效率均大于結構1;且在碗狀底部流道段的換熱功率大于出口流道段,碗狀底部的結構差異對于換熱效率的影響貢獻也大于出口流道段結構差異對于換熱效率的影響貢獻。

圖12 是燃燒室2 種結構冷卻水域的速度流線圖,圖13 是2 種結構冷卻水域螺旋流道處水的流速隨位置的變化規律。從圖12 可看出在螺旋流道處,2 種結構的冷卻水流動均近似于層流流動;從圖13 可看出,在螺旋流道段2 種結構冷卻水的流速分布規律相同且流速均為9 m/s 左右。

圖12 不同結構冷卻水速度流線圖Fig.12 Velocity streamlines of the cooling water area for different structures

圖13 不同結構螺旋流道水流速度分布曲線Fig.13 Cooling water velocity curves in the spiral channel for different structures

從圖12 可知,結構1 碗狀底部5 條流道的左側區域均形成了較大的漩渦區域,右側則存在狹長的近似于層流流動的區域。較大的漩渦區域導致冷卻水滯止無法及時排出冷卻水管道,造成較大面積的局部高溫區。狹長的層流區域則通過對流換熱的形式將熱量帶出,形成溫度相對較低的區域。

結構1 的5 條流道的層流段在碗狀底部最低處交匯,高速水流在底部對沖造成截流現象,速度突變并在第6 條流道以及冷卻水域出口流道處形成大范圍的湍流區與漩渦區,造成冷卻水域出口流道段的冷卻水溫度分布不均,產生許多局部高溫、低溫區。

結構2 碗狀底部各條流道的上半部分存在許多小的漩渦區域,且漩渦區的大小相較于舊結構更小。小的漩渦區導致許多局部的高溫區域產生,其中與出口流道呈90°角的2 條流道的漩渦區域最大,幾乎與流道寬度一致,致使通過此2 條流道流出的冷卻水量相較于其他流道更少,水溫整體較高。

結構2 冷卻水域出口12 條流道的流道寬度較為狹窄,水的流動方式近似于層流,水溫通過對流傳熱過程升高,傳熱距離越長水溫越高。

5 結論

文中通過對比分析燃燒室2 種冷卻水流道結構的冷卻特性,分析其燃燒域及冷卻水域的溫度場及流場,得到以下結論:

1) 在燃燒域內部結構不變的前提下,殼體內冷卻水域的結構變化對于燃燒室燃氣的流動及溫度影響不大。燃氣出口溫差僅為10 K,燃氣流速大小均為61 m/s,分布基本相同。冷卻水域結構的差異對于冷卻水溫度場及流場影響較大。結構2 的整體冷卻效率相較于結構1 提升了8.5%。

2) 隨著冷卻水域結構底部流道數量的增加,冷卻水域底部的湍流及漩渦狀態減弱,冷卻水域的湍流度降低,水的流速增大,存在的局部高溫區域變小。增加底部流道數量,燃燒室換熱功率提高,且影響效果明顯。

3) 隨著出口流道數量以及冷卻水由底部流入橫向出口流道的通道數量增加,出口處的湍流、漩渦及局部高溫區域消失,碗狀底部存在的截流現象消失。增加出口流道段的流道數,燃燒室的換熱功率提高,但對換熱功率的影響效果不如增加碗狀底部流道處的流道數的影響效果明顯。

文中研究證明了結構2 的冷卻換熱效果優于結構1,未來可針對結構2 在不同工況下的殼體熱應力開展數值計算研究,論證結構2 的冷卻水流道結構是否滿足燃燒室的殼體強度使用需求,以確定燃燒室的冷卻換熱結構是否需要進一步優化改進。