射流在不可壓氣流中破碎過程高精度數值仿真

張波濤，張友平，張民慶

(1.西安航天動力研究所，陜西 西安 710100； 2.液體火箭發動機技術重點實驗室，陜西 西安 710100； 3.航天推進技術研究院，陜西 西安 710100)

0 引言

射流破碎后液滴粒徑小且均勻，就有利于氣態燃料和液態燃料的摻混，實現高效率燃燒。為了透徹研究液體射流在不可壓橫向氣流中破碎過程及破碎機理，近年來國內外學者對其進行了大量的研究。

橫向氣流中液體射流破碎過程研究工作主要包括射流柱破碎過程、穿透深度、射流軌跡和液滴分布等。朱英[1]采用高速攝像儀對射流破碎過程進行研究，重點測量了射流柱穿透深度、破碎長度和射流軌跡等參數。Pei-Kuan Wu[2-3]對射流霧化過程進行了實驗研究，得出了不同射流速度和射流角度與穿透深度之間的關系。K.A.Sallam[4]通過高速攝影測量射流柱表面波長、射流柱上出現的袋數和液滴粒徑等參數，重點對低韋伯數橫向氣流中射流破碎過程進行宏觀分析。王延勝[5]利用激光誘導熒光(PLIF)技術對航空煤油的穿透特性進行了實驗，獲得了射流穿透深度關于動量比、韋伯數和軸向距離的經驗關系式。王雄輝[6]采用高速攝像儀對低Weber數橫向氣流中射流破碎過程進行了實驗研究，觀察了表面波現象及射流破碎后形成液滴的尺寸及其速度。李龍飛[7]對真空環境下的射流霧化特性進行研究，得到真空環境下射流的閃蒸是射流破碎主要因素。劉靜[8]和林宇震[9]都對橫向氣流中液體射流霧化過程的研究進展進行了總結。

實驗采用的光學設備不能穿透射流柱，不易捕捉射流細節特征。因此研究者欲通過使用數值模擬對該問題進行研究。當前比較常用的是VOF方法和Level set方法。仝毅恒[10]采用Euler-Euler雙流體模型的VOF方法和CLSVOF方法對圓柱射流破碎過程和破碎位置等特性進行了研究，分析了射流破碎形式及其產生的原因。劉靜[11]采用一次霧化模型和二次霧化模型對超聲速橫向氣流中的液體噴射進行了數值模擬，研究了湍流度和附面層厚度對液霧穿透深度的影響。劉日超[12]等用LES結合VOF的方法，對射流破碎過程進行直接模擬，觀察到射流柱在進入橫向氣流中由于Rayleigh-Taylor(RT)不穩定性和Kelvin-Helmholtz(KH)不穩定性的共同作用迅速發生變形。由于射流破碎后液滴數目眾多且液滴尺寸范圍較大，傳統計算模型中固定不變的網格量難以滿足霧場動態發展過程中不同位置對空間分辨率的需求。基于Gerris的樹形自適應空間離散算法與分段線性的 VOF 方法可以很好的解決這一難點，四叉樹/八叉樹的樹形結構網格使得自適應加密算法可簡易靈活地實現，網格自適應函數可以用疏密程度不同的網格來解決不同空間分辨率的問題，在不損失計算精度的情況下顯著降低了計算量。李佳楠[13]在國內率先使用Gerris軟件對直流互擊噴注單元霧化過程進行數值模擬，分析了撞擊波的形成機理及孔徑比、動量比和射流速度等參數對霧化特性的影響。王凱[14]基于Gerris軟件建立錐形液膜霧化破碎過程數值仿真方法，對相鄰多個離心式噴嘴液膜撞擊霧化過程進行數值仿真，展示了噴霧場三維形態和細節特征，分析了液膜形成原因及液膜撞擊對霧化效果的影響?？梢娀贕erris的高精度數值算法可以很真實的反映霧化過程且捕捉到細節特征，但以往文獻中采用Gerris研究的均是噴注單元外流場為靜止氣體的霧化過程，目前國內還沒有采用Gerris研究射流在流動氣體中霧化的相關文獻。

為深入認識射流在不可壓氣流中的破碎過程，在國內率先采用基于Gerris的全三維VOF方法和樹形自適應的空間離散算法對其研究，可以解決射流柱細節特征和多尺度液滴難捕捉的問題，通過對文獻[1]中射流柱袋式破碎和文獻[2]中射流柱剪切破碎的經典實驗進行數值計算，分析射流在不可壓氣流中的破碎過程和捕捉霧場細節特征，全面驗證基于Gerris數值方法計算射流在不可壓氣流中霧化的高精度性和高可靠性，為后續更進一步認識其霧化機理和評估霧化效果奠定基礎。

1 數學物理模型

根據不可壓橫向氣流中液體射流破碎的物理過程，基于八叉樹結構的正方體單元對空間進行離散，使用網格自適應加密方法處理多尺度問題[15]，采用有限體積法求解不可壓N-S方程和分段線性VOF方法求解氣液界面重構。通過將表面張力轉化為某一區域連續的體積力并結合高度函數曲率估計實現表面張力的精確求解，結合隱式大渦模擬(Implicit Large Eddy Simulation)近似模擬小于最小網格尺度渦的耗散過程。

1.1 控制方程

在計算中假定氣液流動過程是等溫的且不考慮蒸發過程，因此無需求解能量方程。求解的不可壓、兩相流、含有表面張力的N-S方程組為:

1.2 VOF方法

VOF方法最早由Hirt和Nichols[16]提出，對于兩相流動，通過定義第一相體積分數α(x，t)來描述氣液界面。將密度連續方程用體積分數α的連續方程替換為：

密度和粘度也可由體積分數表示為:

式中：ρ1，ρ2和μ1，μ2分別為第一相和第二相的密度和粘度。

1.3 計算域與參數設置

Gerris計算三維算例的網格均為正方體單元，如圖1所示每個Box由前、后、左、右、上、下共六個面組成，計算域可以根據計算模型由多個Box塊連接而成。計算區域可以不用三維建模軟件建立，直接在參數腳本文件中給定Box數目、Box連接方向和計算模型在計算域的位置。為了提高計算精度，文中采用八個邊長為11.25 mm的Box，計算域如圖2所示，上下兩層，每層4個Box，計算域長、寬、高分別為45 mm，11.25 mm和22.5 mm，橫向氣流以一定的速度從兩層最左側Box的左側面流入，液體射流以一定角度從下表面上的圓孔射入。計算時最高網格等級設為10級，最小網格為10.99 μm，計算域和邊界條件的設置均通過編寫大量代碼實現。

圖1 基本塊Box離散方向的定義Fig.1 Definition for discrete direction of basic block Box

圖2 計算域Fig.2 Calculation domain

射流破碎過程仿真邊界條件如表1所示，算例1和算例2屬袋式破碎，在算例1與文獻[1]中同工況實驗結果對比驗證計算方法準確的基礎之上，算例2捕捉射流柱表面細節特征。算例3～算例5屬剪切破碎，在與文獻[2]中同工況實驗對比驗證后，分析射流角度對破碎效果的影響。

橫向氣流中液體射流破碎過程的數值模擬屬于多相流模擬的范疇，在計算過程中需要給定的物性參數有液相的密度、粘度系數、表面張力系數、氣相的密度及粘性系數等參數。數值仿真中液相為水，氣相為空氣。物性參數如表2所示。

表1 邊界條件

表2 物性參數表

2 結果分析

2.1 計算結果與實驗結果對比

圖3給出了算例1工況下數值仿真結果和文獻[1]中實驗結果對比圖，計算所得到的霧場宏觀結構與實驗結果非常相近。

圖3 算例1的計算結果與試驗結果對比Fig.3 Comparison of numerical calculation and experiment results of Case 1

圖4為算例3～算例5工況下的數值仿真結果與文獻[2]中實驗結果的對比圖，圖4中顏色代表速度大小，從圖4中可以看出，隨著射流角度的增大，射流柱縱向穿透深度增大，橫向破碎長度減小，由于氣流速度與射流柱橫向分速度差變大，射流柱上的波動和慣性力更加明顯，霧化效果較好。當角度較小時，射流柱上的氣動力相對作用將減小，射流柱變直，接近于自由射流。

圖4 算例3～算例5的計算結果與試驗結果對比Fig.4 Comparison of numerical calculation and experiment results of Case 3， Case 4 and Case 5

為了進一步驗證數值計算的準確性，本文取射流柱外圍為射流軌跡，將數值計算得到的射流軌跡結果與實驗結果對比，圖5(a)和圖5(b)分別為袋式破碎和剪切破碎射流軌跡圖，數值計算結果與實驗結果相對誤差在5%左右。

2.2 射流柱袋式破碎過程

圖6為射流柱在橫向氣流作用下變形、彎曲到進一步破碎成大液滴的破碎全過程，由于采用VOF相界面捕捉方法，顯示結果是體積分數為α=0.5的等值面。從圖6中可以看出，射流柱在初始發展階段沒有明顯的液滴剝離，只在射流柱表面有小褶皺，隨著射流柱的發展，射流柱表面的小褶皺在氣動力作用下剝離形成小液滴，接著迎風面向內凹陷形成波谷，隨著射流柱的變形，迎風面和背風面間的波谷厚度不斷變薄，有袋式破碎的特征，最終在表面張力和氣動力共同作用下斷開，形成小液塊和液滴。

圖5 液體射流軌跡Fig.5 Trajectories of liquid jets

圖6 袋式破碎過程數值計算結果Fig.6 Numerical calculation results of bag breakup process of jet column

為了更清晰的觀察射流柱表面特征和截面特征，對算例2中射流柱彎曲前截面進行分析。將三維霧場沿著射流柱所在的xy平面和yz平面剖切，得到了前視圖和左視圖的射流柱二維體積分數分布云圖，如圖7所示。圖7(a)顯示了射流柱在橫向氣流的作用下沿射流方向越來越窄，圖7(b)顯示了射流柱沿射流方向逐漸變寬，迎風面稍有變大。這是因為橫向氣流遇到射流柱發生滯止，在迎風面上形成一個滯留區，滯留區內速度為零，使得迎風面壓力增大，同時射流柱背風面壓力較低。在壓差作用下，射流柱變得寬而薄，隨著射流柱變形，射流柱受到氣流作用力越來越大，導致射流柱在流動方向表現為向下游彎曲。

為了進一步認識射流柱破碎機理，圖8是不同視角觀察到的射流柱彎曲前迎風面放大圖，可以清晰看到射流柱表面形成的表面波，表面波是導致射流柱不穩定的根本原因。由于射流柱左側迎風面受到氣流垂直作用，氣液密度差和交界面上的垂直加速度導致RT波的產生，RT波將射流柱兩側液體向前推移，以致射流柱變得薄且向下游彎曲。橫向氣、液相對速度由1逐漸變小，氣液橫向速度差導致KH波的產生。橫向氣、液相對速度差也是影響射流柱破碎的重要參數。箭頭指出了射流柱發展過程中在射流柱表面產生的表面波，1為射流柱在射流發展方向與氣流存在速度差通過相互剪切而生成的KH表面波，2為射流柱在橫向與氣流存在速度差產生的KH表面波，隨著KH波的發展射流柱側面有液滴剝離，從圖中也可以清晰的看出表面波處顏色較深，即表面波處速度較大。

圖7 體積分數云圖Fig.7 Contours of volume fraction α

圖8 射流柱表面速度分布圖Fig.8 Contours of velocity distribution on jet column surfaces

2.3 射流柱剪切破碎過程

以30°為例分析射流柱剪切破碎過程，圖9是射流角度為30°工況下的射流破碎過程圖，在射流柱剛噴入橫向氣流中的初始階段，即有小液滴從射流柱表面剝離，這一現象始終貫穿于射流柱發展的全過程。射流柱迎風面和背風面均出現小尺度波動，但迎風面表面波尺度大于背風面表面波尺度，這是由于迎風面受到的氣動力更大。小尺度波動在氣動力作用下發展為明顯的大尺度波動，最終由于氣動力大于射流柱表面張力和粘性力而導致射流破碎。這類似于射流柱袋式破碎在氣動力作用下迎風面向內凹陷形成波谷并被拉長后斷裂的特征。同時，隨著射流柱的發展，射流柱逐漸破碎成許多大液塊，發生一次霧化，而不像射流柱袋式破碎先斷裂為單一大液塊后在氣動力的作用下發生二次霧化破碎為眾多小液滴。在整個射流柱破碎過程中，射流柱剪切破碎和袋式破碎一樣存在R-T和K-H兩種不穩定波，但剪切破碎的不穩定波更加明顯。

圖10～圖12分別表征的是射流中間截面的體積分數云圖、速度云圖和渦量云圖，從圖中可以看出射流在出口處即有小液滴剝離，射流柱迎風面在氣動力的作用下出現小尺度波，同時可以看到在射流柱迎風面很薄的區域內速度很小，這是因為橫向氣流遇到射流柱發生滯止，速度減小，壓強增大，背風面比迎風面速度大，壓強小，迎風面的壓強大于背風面的壓強使射流柱彎曲。隨著射流柱的發展，非穩態波動變大，并形成大量渦結構，加劇霧化過程，縮短射流柱破碎長度。

Gerris可以對全場所有的液滴數目、體積進行統計后并輸出。算例3～算例5工況下的全場液滴SMD分別為60.23 μm，59.96 μm和59.12 μm，射流角度對粒徑影響很小。圖13是射流角度為30°時計算的全場液滴粒徑分布圖，從圖13中可以看出，液滴粒徑分布范圍廣，但呈現兩端少、中間多的分布。

圖9 射流柱剪切過程數值計算結果Fig.9 Numerical calculation results of shear breakup process of jet column

圖10 射流柱中心截面體積分數云圖Fig.10 Volume fraction contour on center section of jet column

圖11 射流柱中心截面速度云圖Fig.11 Velocity contour on center section of jet column

圖12 射流柱中心截面渦量云圖Fig.12 Vorticity contour on center section of jet column

圖13 全場液滴粒徑概率密度分布Fig.13 Probability density distribution of droplet diameter in whole flow field

3 結論

通過數值求解三維不可壓N-S 方程組對射流在不可壓橫向氣流中的破碎過程開展數值模擬工作，得出的主要結論如下:

1)通過采用分段線性的VOF方法和樹形自適應加密算法對朱英和P K Wu的實驗進行數值模擬，可以詳細展現整個破碎過程及霧場細節特征，計算得到的射流軌跡和液滴空間分布都與實驗結果吻合，全面驗證了基于Gerris數值方法計算射流在不可壓橫向氣流中霧化的高可靠性和高精度性。

2)射流柱剪切破碎全場液滴平均直徑約為60μm左右，隨著角度從30°增大到75°粒徑沒有明顯變化，射流角度對粒徑影響較小，由于氣流速度射流橫向分速度差變大，只會加快射流霧化。

3)射流柱在氣動力和表面張力共同作用下產生RT波和KH波，表面波是導致射流柱不穩定的根本原因，表面波處的速度較大，射流柱剪切破碎比袋式破碎形成的表面波更多。

本文數值計算因重點放在驗證Gerris數值算法計算射流在流動氣體中破碎的高可靠性和高精度性，沒有定量分析射流柱表面細節特征，因此對射流柱表面細節特征定量分析作為下一步的工作，本文的工作也進一步為應用和改進Gerris代碼奠定了基礎。

參考文獻：

[1] 朱英， 黃勇， 王方， 等. 橫向氣流中的液體圓形射流破碎實驗[J]. 航空動力學報， 2010， 25(10): 2216-2266.

[2] WU Peikuan， KEVIN A. Breakup processes of liquid jets in subsonic crossflows: AIAA 1996-3024 [R]. USA: AIAA， 1996.

[3] FULLER R P， WU Peikuan. Effects of injection angle on atomization of liquid jets in tranverse airflow [R]. AIAA journal， 2000， 38(1): 2216-2266.

[4] SALLAM K A. Bag breakup of nonturbulent liquid jets in crossflow [J]. International journal of multiphase flow， 2008， 34(2008): 64-72.

[5] 王延勝， 林宇震， 李林， 等. 基于PLIF技術的航空煤油橫向射流穿透深度研究[J]. 推進技術， 2015， 36(9): 1395-1403.

[6] 王雄輝， 黃勇， 王方， 等. 橫向氣流中液體圓柱射流的破碎特性和表面波現象[J]. 航空動力學報，2012，27 (9):1979-1987.

[7] 李龍飛， 池保華， 楊偉東， 等. 真空環境下液體射流霧化特性的實驗研究[J]. 火箭推進， 2010， 36(5): 27-30.

LI Longfei， CHI Baohua， YANG Weidong， et al. Experimental study on atomization characteristic of liquid jet in vacuum environment [J]. Journal of rocket propulsion， 2010， 36(5): 27-30.

[8] 劉靜， 徐旭. 高速氣流中橫向液體射流霧化研究進展[J]. 力學進展， 2009， 39(3): 273-283.

[9] 林宇震， 李林， 張弛， 等. 液體射流噴入橫向氣流混合特性研究進展[J]. 航空學報， 2014， 35(1): 46-57.

[10] 仝毅恒. 橫向氣流中液體射流噴注特性和破碎過程研究[D]. 長沙：國防科學技術大學， 2012.

[11] 劉靜， 徐旭. 超聲速橫向氣流中噴霧的數值模擬[J]. 火箭推進， 2006， 32(5): 32-36.

LIU Jing， XU Xu. Numerical simulation of liquid jet in supersonic cross flows[J]. Journal of rocket propulsion， 2006， 32(5): 32-36.

[12] 劉日超， 樂嘉陵， 楊順華， 等. 亞聲速橫向氣流中液體射流破碎過程的直接模擬[J]. 推進技術， 2016， 37 (11): 2135-2141.

[13] 李佳楠， 費俊， 楊偉東， 等. 直流互擊式噴注單元霧化特性準直接數值模擬[J]. 推進技術， 2016， 37 (4): 713-724.

[14] 王凱， 李鵬飛， 楊國華， 等. 相鄰離心式噴嘴液膜撞擊霧化過程仿真[J]. 推進技術， 2017， 38 (2) : 408-415.

[15] GERRIS P S. A tree-based adaptive solver for the incompressible euler equations in complex geometries [J]. Journal of computational physics， 2003， 190(2): 572-600.

[16] HIRT C W，NICHOLS B D. Volume of fluid(VOF)method for the dynamics of free boundaries [J]. Journal of computational physics， 1981， 39(1): 201-225.