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

張波濤，張友平，張民慶

(1.西安航天動(dòng)力研究所，陜西 西安 710100； 2.液體火箭發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710100； 3.航天推進(jìn)技術(shù)研究院，陜西 西安 710100)

0 引言

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

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

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

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

1 數(shù)學(xué)物理模型

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

1.1 控制方程

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

1.2 VOF方法

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

密度和粘度也可由體積分?jǐn)?shù)表示為:

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

1.3 計(jì)算域與參數(shù)設(shè)置

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

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

圖2 計(jì)算域Fig.2 Calculation domain

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

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

表1 邊界條件

表2 物性參數(shù)表

2 結(jié)果分析

2.1 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比

圖3給出了算例1工況下數(shù)值仿真結(jié)果和文獻(xiàn)[1]中實(shí)驗(yàn)結(jié)果對(duì)比圖，計(jì)算所得到的霧場宏觀結(jié)構(gòu)與實(shí)驗(yàn)結(jié)果非常相近。

圖3 算例1的計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.3 Comparison of numerical calculation and experiment results of Case 1

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

圖4 算例3～算例5的計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.4 Comparison of numerical calculation and experiment results of Case 3， Case 4 and Case 5

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

2.2 射流柱袋式破碎過程

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

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

圖6 袋式破碎過程數(shù)值計(jì)算結(jié)果Fig.6 Numerical calculation results of bag breakup process of jet column

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

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

圖7 體積分?jǐn)?shù)云圖Fig.7 Contours of volume fraction α

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

2.3 射流柱剪切破碎過程

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

圖10～圖12分別表征的是射流中間截面的體積分?jǐn)?shù)云圖、速度云圖和渦量云圖，從圖中可以看出射流在出口處即有小液滴剝離，射流柱迎風(fēng)面在氣動(dòng)力的作用下出現(xiàn)小尺度波，同時(shí)可以看到在射流柱迎風(fēng)面很薄的區(qū)域內(nèi)速度很小，這是因?yàn)闄M向氣流遇到射流柱發(fā)生滯止，速度減小，壓強(qiáng)增大，背風(fēng)面比迎風(fēng)面速度大，壓強(qiáng)小，迎風(fēng)面的壓強(qiáng)大于背風(fēng)面的壓強(qiáng)使射流柱彎曲。隨著射流柱的發(fā)展，非穩(wěn)態(tài)波動(dòng)變大，并形成大量渦結(jié)構(gòu)，加劇霧化過程，縮短射流柱破碎長度。

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

圖9 射流柱剪切過程數(shù)值計(jì)算結(jié)果Fig.9 Numerical calculation results of shear breakup process of jet column

圖10 射流柱中心截面體積分?jǐn)?shù)云圖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 結(jié)論

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

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

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

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

本文數(shù)值計(jì)算因重點(diǎn)放在驗(yàn)證Gerris數(shù)值算法計(jì)算射流在流動(dòng)氣體中破碎的高可靠性和高精度性，沒有定量分析射流柱表面細(xì)節(jié)特征，因此對(duì)射流柱表面細(xì)節(jié)特征定量分析作為下一步的工作，本文的工作也進(jìn)一步為應(yīng)用和改進(jìn)Gerris代碼奠定了基礎(chǔ)。

參考文獻(xiàn)：

[1] 朱英， 黃勇， 王方， 等. 橫向氣流中的液體圓形射流破碎實(shí)驗(yàn)[J]. 航空動(dòng)力學(xué)報(bào)， 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技術(shù)的航空煤油橫向射流穿透深度研究[J]. 推進(jìn)技術(shù)， 2015， 36(9): 1395-1403.

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

[7] 李龍飛， 池保華， 楊偉東， 等. 真空環(huán)境下液體射流霧化特性的實(shí)驗(yàn)研究[J]. 火箭推進(jìn)， 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ìn)展[J]. 力學(xué)進(jìn)展， 2009， 39(3): 273-283.

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

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

[11] 劉靜， 徐旭. 超聲速橫向氣流中噴霧的數(shù)值模擬[J]. 火箭推進(jìn)， 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]. 推進(jìn)技術(shù)， 2016， 37 (11): 2135-2141.

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

[14] 王凱， 李鵬飛， 楊國華， 等. 相鄰離心式噴嘴液膜撞擊霧化過程仿真[J]. 推進(jìn)技術(shù)， 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.