流動拓撲對液體射流失穩與霧化的影響

楊立軍 黃東騏 韓旺 李敬軒 富慶飛

(1. 北京航空航天大學 宇航學院, 北京 100083; 2. 北京航空航天大學寧波創新研究院, 寧波 315800)

推進劑霧化包括液體射流失穩、斷裂、破碎、液滴生成等一系列過程。 推進劑的霧化性能會顯著影響發動機性能;研究射流霧化機理有助于實現對霧化過程的精準調控,從而優化后續的蒸發燃燒過程,提升發動機的推力性能。 盡管前人對液體射流霧化開展了大量的研究,但對其霧化機理還存在認識不清的地方[1]。

射流霧化的研究方法主要分為2 類:實驗測量和數值模擬。 實驗測量方法的優勢在于容易開展參數可調的大規模研究,但也面臨著測量精度有限,測量物理量較少等挑戰。 數值模擬方法雖然計算量較大,但能夠提供更高的精度,同時也能提供更多實驗中不易獲得的數據,幫助開展更加細致的研究,因此受到廣泛的關注。

在霧化的直接數值模擬(DNS)研究中,通常采用水平集[2]、流體體積[3]等方法來捕捉氣液界面。 例如,Leboissetier 和Zaleski[4]最早對柴油燃料射流進行了DNS,該研究表明噴嘴內部的液體射流在層流條件下不會發生破裂霧化。 這一結果與實驗對比顯示出了較為良好的一致性。 對于噴嘴外部的射流,Klein[5]指出,當平面液體射流的韋伯數小于270 時,沒有發現射流破裂及霧化的行為。 Sander 和Weigand[6]發現,除了噴嘴出口處的湍流特性會對射流液面的穩定性造成影響,速度型也會顯著影響射流霧化的結果。 Ménard等[7]首次在兩相環境下研究了柴油燃料射流DNS的網格分辨率問題;Desjardins 和Pitsch[8]則使用DNS 研究了中等韋伯數下(We=500 ～2 000)湍流射流液面的破裂行為。 Shinjo 和Umemura[9]在層流環境及靜止空氣中,在大韋伯數下(We=14 000)使用高分辨率網格對圓柱柴油燃料射流進行了數值模擬。

以上數值模擬研究展現了較為清晰的射流霧化過程,發現韋伯數和流場環境是影響射流霧化的重要因素,但是有關射流霧化的詳細機制仍然不清。 基于此,本文嘗試通過流動拓撲理論分析液體射流霧化機理。 流動拓撲理論由Chong等[10-11]提出。 流動拓撲理論主要是基于速度梯度張量的3 個不變量來定義的不同流場拓撲結構,這種方法被廣泛運用于解決流體問題及燃燒問題。 Chong 等[12]運用該方法研究了壁面剪切流動,Elsinga 和Marusic[13]研究了均勻且各向同性的湍流流動。 Hasslberger 等[14]首次在氣泡空氣和水的兩相流流動中研究了流動拓撲的分布,并發現代表旋轉運動的拓撲結構主要存在于氣泡內部。 Watanabe 等[15]研究了在湍流與非湍流界面附近的拓撲結構。 Josef 等[16]研究了柴油燃料射流霧化過程中的流動拓撲結構影響;Han 等[17]研究了湍流流動拓撲結構對湍流非預混火焰切向擴散的影響,從流動拓撲結構角度給出了經典火焰面模型假設不成立的條件。

基于上述討論,本文采用水平集方法,對湍流液體平面射流直接進行三維數值模擬。 水平集方法有2 個優點:①相比于傳統的歐拉法,水平集方法可以直接在直角網格中對曲線和曲面進行數值計算,而不用將曲線和曲面參數化;②水平集方法易于追蹤物體結構的改變。 因此,使用水平集方法進行DNS 可以有效提高計算精度,并減少計算成本。 本文通過使用水平集方法得到高分辨率數值模擬數據,從流動拓撲角度揭示液體平面射流霧化的機理。

1 模型與算法

1.1 控制方程與算例設置

考慮在靜止氣體環境中的液體平面射流,其液體射流湍流初速度場在一個槽道中產生(見圖1,右側液體初始湍流場由左側充分發展的槽道流速度場提供)。 射流的質量方程和動量方程由經典的不可壓縮Navier-Stokes 方程給出:

圖1 湍流液體平面射流算例設置Fig.1 Turbulent liquid plane jet arithmetic setup

式中:u為速度場;ρ為密度;p為壓強;μ為黏度。

界面Γ將液相和氣相分隔開。 值得注意的是,在氣液多相流中,物理量在兩相界面處存在突變。 令液相密度ρ=ρl,氣相密度ρ=ρg,液相動力黏度μ=μl,氣相動力黏度μ=μg。 因此,在界面Γ處滿足:[ρ]Γ=ρl-ρg,[μ]Γ=μl-μg。 而由于速度連續,在界面處[u]Γ=0。 影響氣液界面傳輸計算精度的主要是如下的壓力跳躍:

[P]Γ=σκ+ 2[μ]ΓnT·(Δu)·n(3)

式中:σ為表面張力系數;κ為界面曲率;n為界面法向量。

本文使用低馬赫多相湍流燃燒求解器NGA[18]對式(1) ～式(3)方程進行數值求解。 在進行數值模擬時,為解決界面處的壓力跳躍和密度跳躍,本文使用了精確水平集守恒方法ACLS,在處理表面張力時應用了虛擬流動方法GFM[19]。

圖1 展示了流動計算域。 射流初始高度為H,入口處平均速度為U0,整個計算域在流向和展向上(即x,z方向)是周期性邊界條件。 在流動高度方向(即y方向)是出口邊界條件。 整個計算域在x-y-z方向上尺寸是4H×6H×4H,對應每個方向的網格數目是256 ×384 ×256,網格大小為H/64[8]。 而提供湍流初速度流場的槽道流計算域如圖1 所示,其計算尺寸是4H×H×4H,對應每個方向的網格數目是256 ×64 ×256,槽道流的高度為H,與后續射流初始高度對應。 射流初始高度處于-H/2 ＜y＜H/2 的范圍內,液體射流的雷諾數Re=5 000,韋伯數We=1 000,其中雷諾數Re=ρU0H/μ,韋伯數We=ρU20H/σ。 氣液密度比被設置為ρl/ρg=40,為使氣液運動黏度相等,設置μl/μg=40。

圖2 展示了不同時刻下氣液界面的演化。 本文氣液界面通過液體體積百分數為0.5 的等值面表征,圖2(a) ～(d)分別對應的時間是0. 01,0.03,0.06,0.10 ms 時刻,對應界面失穩,表面波增長,界面出現破碎,射流初次霧化的時間點。

圖2 射流氣液界面的演化Fig.2 Evolution of jet gas-liquid interface

1.2 流動拓撲理論

本節介紹流動拓撲分析方法。 根據Perry 和Chong 等[10-12]的研究,速度梯度張量Aij由式(4)給出:

同時,可將速度梯度張量分解為Sij和Wij兩部分之和,其中Sij=0.5(Aij+Aji),Wij=0.5(Aij-Aji),分別稱為對稱張量和反對稱張量。 速度梯度張量Aij的3 個不變量由式(5)給出:

顯然,第1 不變量P代表流體體積的變化,由質量守恒知P=0,因此本文將主要在Q-R平面內分析流動拓撲結構。 式(5)是一元三次方程,其判別式為

判別式將Q-R平面分為2 個區域,當Δ＞0,方程有1 個實數根和2 個復數根,當Δ＜0,方程有3 個實數根,當Δ=0 時,對應2 條直線,這2 條直線分離不同的拓撲結構區域(見圖3):

當2 個復數根為純虛數時,對應r2=0。 這3條直線r1a、r1b、r2將整個Q-R平面分為4 個區域,分別命名為S1 ～S4,如圖3 所示。 利用速度梯度張量的第2 不變量Q可以解釋每個部分所表示的物理意義,Q可以被分成2 個部分:

圖3 不同拓撲區域分布示意圖Fig.3 Different topology area distribution diagram

式中:QS和QW分別為Sij和Wij的第2 張量 不 變量。QS代表每單位質量的動能耗散為熱能,代表變形;QW與渦量相關,代表旋轉,因此Q＞0,代表系統以旋轉為主,Q＜0,則代表系統以變形為主。同時R＞0,代表壓縮行為,R＜0,代表拉伸行為。通過流動拓撲分析,可以看出在不同的拓撲區域內,流動受到不同的機制控制。 總而言之,S1 拓撲區域代表不穩定焦點結構(UFC),主要受旋轉壓縮機制影響。 S2 拓撲區域代表不穩定節點/鞍點結構(UN/S/S),主要受到壓縮變形機制影響。S3 拓撲區域代表穩定節點/鞍點結構(SN/S/S),主要受到拉伸變形機制影響。 S4 拓撲區域代表穩定焦點結構(SFS),主要受到旋轉拉伸機制影響。

1.3 氣液界面曲率

為了研究流動拓撲與氣液界面的耦合關系,還需研究氣液界面曲率。 下面對界面曲率給出定義。 本文中所考慮的界面是液體體積分數等值面,以α表示。α=1 時表示全為液相,α=0 時表示全為氣相。 該等值面的法向量n定義為

圖4 使用平均曲率和高斯曲率對等值面結構進行分類[11]Fig.4 Classification of structure of equivalence surface using mean curvature and Gaussian curvature[11]

2 結果與討論

本節討論液體射流霧化過程中流動拓撲結構演化以及流動拓撲結構對氣液界面曲率的影響。首先,圖5 展示了第2 不變量Q和第3 不變量R的聯合概率密度分布圖(JPDF)。 為方便對比,Q使用QW進行無量綱化(Q/ ＜QW＞),同理R使用QW的3/2 次方進行無量綱化(R/ ＜QW＞3/2)。選取t=0.01 ms 和t=0.10 ms 的結果進行展示。如圖5所示,2 個時刻最大的概率密度分布均出現在原點附近,因此將圖5 中的結果原點處在圖6中進行放大展示。 可以看出,在液體射流初始時刻,S1 拓撲結構出現的概率略微不同,S1 概率稍大,而隨著時間推移,到了初次霧化的時刻,S1 拓撲結構出現的概率較大,表明流動主要受UFC 拓撲結構的影響,即在霧化過程中,流體微團主要受旋轉壓縮作用機制的影響;而且隨著霧化過程的進行,這種機制的重要性越發明顯。 值得注意的是,不同時刻的圖像與之前研究中觀察到的水滴形狀比較相似,且在霧化后,區域面積變化較小,呈現出一定的自相似性質。

圖5 第2、第3 張量不變量的聯合概率密度分布Fig.5 Joint probability density function of the second and the third tensor invariance

圖6 圖5 在原點附近的放大圖Fig.6 Fig.5 enlarged view near origin point

為了進一步研究流動拓撲對不同液體體積分數等值面包括氣液界面的影響,本文引入條件應變率來對結果進行分析,定義如下:

式中:a為等值面切向應變率,該值的正與負分別對應于拉伸應變率和壓縮應變率。

圖7 展示了射流霧化過程中應變率在Q-R空間的條件平均分布。 如圖7 所示,在UFC(S1)拓撲分區,主要引起拉伸應變率(應變大于0)。而在SFS(S4)拓撲分區,主要引起壓縮應變率(應變小于0)。 而在UN/S/S(S2)、SN/S/S(S3)區域,則由應變主導霧化。 另外,圖7 表明,初始時刻,應變能主要是拉伸應變率主導。 而隨著霧化進行,壓縮應變率與拉伸應變率效應共存。 其中,UFC(S1)流動拓撲區域產生拉伸應變率的概率最高。

圖7 不同時刻下應變率在Q-R 平面內的條件平均分布Fig.7 Conditional average distribution of strain rate in Q-R plane in different moments

圖8 展示了同一時刻t=0.06 ms 下,不同液體體積分數(α=0.2,0.5,0.8)等值面上應變率的差異。 如圖8 所示,對于UFS 拓撲結構,當氣相占比較大時(α= 0. 2),具有較高的拉伸應變率,而液相占比較大(α=0.8),拉伸應變率較小。對于其他流動拓撲結構,均為產生壓縮應變率的概率較高。 其中SFS 結構,壓縮應變率較小,而在變形主導的區域(S2,S3),壓縮應變率較大。

圖8 不同等值面應變率在Q-R 平面內的條件平均分布Fig.8 Conditional average distribution of strain rate in Q-R plane on different iso-surfaces

選取同一時刻t=0.06 ms 下不同的液體體積分數等值面(α=0.2,0.5,0.8) 曲率進行分析,通過κg-κm聯合概率密度分布圖像探討不同界面的幾何結構。 如圖9 所示,等值面呈現出不同的結構。 在氣液界面處(α=0.5 處),概率密度分布較為對稱,最高概率密度出現在原點附近,說明主要呈現出片狀的結構。 而在α=0. 2 時,氣相占比較大,平均曲率大于0 概率較大,這種情況下主要呈現出凸形表面結構。 在α=0. 8 時,液相占比較大,平均曲率小于0 概率較大,主要呈現出凹形表面結構,這些表明,局部界面幾何結構會受到局部拓撲結構影響很大。

圖9 不同等值面上的κg-κm 聯合概率密度分布Fig.9 Joint probability density distribution of κg-κm on different equivalence sufales

取不同時刻研究所有液體體積百分數等值面上的κg-κm聯合概率密度分布,以此探討霧化過程的不同時刻整個射流流場的曲率分布,從而獲得整個流場界面的大致形狀。 除了選取初次霧化之前的t=0.01ms,t=0.03 ms 和t=0.06 ms,還選取了霧化發生后的時間t=0.10 ms,t=0.15 ms和t=0.20 ms。 如圖10(a) ～(f)所示,霧化初始時刻主要存在的是凸面、凹面、片狀結構,其中片狀結構概率最高。 隨著霧化進行,出現了管狀的等值面結構,同時平均曲率及高斯曲率的絕對值不斷增加,這意味著隨著霧化進程,也會出現凹雙曲面、凹拋物面、凸雙曲面、凸拋物面結構,但在每個時刻,仍然是出現片狀結構的概率最高。 而在初次霧化后,逐漸呈現出圖像不關于κm=0 軸對稱的現象,產生凹雙曲面的概率略微大于產生凸雙曲面的概率。 由圖9 可知,液相占比較大時產生凹面的概率較大,可能是由于霧化后期液滴彌散至整個空間所致。

圖10 不同時刻下整個流場的κg-κm 聯合概率密度分布Fig.10 Joint probability density distribution of κg-κm for whole flow field at different moments

圖11 展示了應變率在不同拓撲區域內的概率分布圖(PDF)。 可以看出,所有流動拓撲結構都有助于在不同區域內形成壓縮應變率和拉伸應變率。 但在每個區域內,產生零應變的概率均為最高。 而隨著霧化進行,每個拓撲區域產生的壓縮應變率都會大于其產生的拉伸應變率,即偏向于負的應變率。 由圖11 可以看出,S1 和S3 的概率密度顯著大于S2 和S4 的概率密度。 隨著霧化進行,S1 拓撲會產生更大的壓縮應變率。 另外,S1 和S3 對拉伸與壓縮應變率的貢獻最高,由于S1 拓撲比S3 拓撲概率更大,S1 概率增大,應變率變大,曲率變小。

圖11 不同拓撲區域內的應變率概率密度分布Fig.11 Probability density distribution of strain rate in different topological regions

圖12 展示了平均曲率在不同拓撲區域內的概率分布。 可以看出,所有流動拓撲結構都會產生正平均曲率和負平均曲率,每個區域內,產生零曲率的概率最大,表明產生片狀結構的概率最高,這與之前的結論一致。 同時,曲率的概率分布大致關于κm=0 對稱,表明在霧化過程中,產生凸界面結構與凹界面結構的概率大致相等。

圖12 不同拓撲區域內的平均曲率概率密度分布Fig.12 Mean curvature probability density distribution in different topological regions

圖13 展示了t=0.10 ms 時刻應變率和平均曲率的聯合概率密度函數分布。 圖13 表明,隨著應變增加,等值面界面趨向于片狀,管狀結構,在壓縮應變率區域(應變小于0),曲率下降的速度比拉伸區域(應變大于0)快,在壓縮應變率區域,應變和曲率的負相關關系更為明顯。 這里負相關關系指一個變量隨著另一個變量的減少而增加的關系?；谏鲜鼋Y論,可以總結出如下因果鏈,其中,上下箭頭代表該量絕對值的增加或者減少,水平箭頭代表因果鏈:PDF(UFC)↑→應變率↑→曲率↓→霧化性能↓。

圖13 整個流場中應變率與平均曲率的聯合概率密度函數分布Fig.13 Joint probability density function distribution of strain and mean curvature over the entire flow field

3 結 論

本文對湍流液體平面射流進行了三維直接數值模擬。 通過流動拓撲理論,分析了射流霧化過程中流動拓撲的影響機制。 通過分析氣液界面的曲率,揭示了流動拓撲與界面結構之間的關系。從流動拓撲的角度闡述了射流霧化氣液界面不斷演化的機制。 研究結果表明:

1) 所有流動拓撲結構都有助于產生壓縮應變率和拉伸應變率。 射流霧化過程主要受到UFC 拓撲結構的影響,從而產生拉伸應變率,且UFC 的增多會導致更多的等值面界面壓縮應變率,進而導致氣液界面曲率下降,降低霧化性能。

2) 隨著霧化進行,氣相區域主要產生凸界面,液相區域主要產生凹界面,氣液界面處界面主要為片狀結構。 除UFC 外,其余拓撲結構主要產生壓縮應變率,在SFS 區域產生的壓縮應變率較小,而在UN/S/S 區域和SN/S/S 區域,產生的壓縮應變率較大。

3) 隨著霧化進行,產生壓縮應變率的概率會略微大于拉伸應變率概率。

根據上述結論, 總結出如下的因果鏈:PDF(UFC)↑→應變率↑→曲率↓→霧化性能↓。該因果鏈說明,隨著UFC 拓撲結構的增加,會使得拉伸應變率增大,界面受到拉伸,曲率則會下降,不利于霧化過程的進行,因此霧化性能下降率。 為了驗證以上因果鏈的普適性,本文作者團隊正在進行多參數變化的直接數值模擬研究,旨在全面地揭示不同雷諾數與韋伯數下流動拓撲對液體射流霧化的影響機制。