流動拓撲對液體射流失穩(wěn)與霧化的影響

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

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

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

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

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

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

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

1 模型與算法

1.1 控制方程與算例設(shè)置

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

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

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

界面Γ將液相和氣相分隔開。 值得注意的是,在氣液多相流中,物理量在兩相界面處存在突變。 令液相密度ρ=ρl,氣相密度ρ=ρg,液相動力黏度μ=μl,氣相動力黏度μ=μg。 因此,在界面Γ處滿足:[ρ]Γ=ρl-ρg,[μ]Γ=μl-μg。 而由于速度連續(xù),在界面處[u]Γ=0。 影響氣液界面?zhèn)鬏斢嬎憔鹊闹饕侨缦碌膲毫μS:

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

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

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

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

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

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

1.2 流動拓撲理論

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

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

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

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

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

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

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

1.3 氣液界面曲率

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

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

2 結(jié)果與討論

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

圖5 第2、第3 張量不變量的聯(lián)合概率密度分布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

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

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

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

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

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

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

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

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

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

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

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

圖11 不同拓撲區(qū)域內(nèi)的應(yīng)變率概率密度分布Fig.11 Probability density distribution of strain rate in different topological regions

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

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

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

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

3 結(jié) 論

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

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

2) 隨著霧化進行,氣相區(qū)域主要產(chǎn)生凸界面,液相區(qū)域主要產(chǎn)生凹界面,氣液界面處界面主要為片狀結(jié)構(gòu)。 除UFC 外,其余拓撲結(jié)構(gòu)主要產(chǎn)生壓縮應(yīng)變率,在SFS 區(qū)域產(chǎn)生的壓縮應(yīng)變率較小,而在UN/S/S 區(qū)域和SN/S/S 區(qū)域,產(chǎn)生的壓縮應(yīng)變率較大。

3) 隨著霧化進行,產(chǎn)生壓縮應(yīng)變率的概率會略微大于拉伸應(yīng)變率概率。

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