金烜 沈赤兵

1) (國防科技大學空天科學學院,高超聲速沖壓發(fā)動機技術重點實驗室,長沙 410073)

2) (中國空氣動力研究與發(fā)展中心,綿陽 621000)

1 引言

變推力液體火箭發(fā)動機在空間運輸、交會對接、星際著陸等方面具有顯著的優(yōu)越性[1],而針栓式噴注器作為變推力發(fā)動機應用的典型,相對于傳統(tǒng)的同軸離心式或直流撞擊式噴注器具有結構簡單、工況和推進劑組合適應性良好、燃燒高效穩(wěn)定等優(yōu)點[2].上述優(yōu)點使針栓式噴注器獲得了一系列成功應用,包括先后將12 名宇航員送上月球的阿波羅登月艙下降發(fā)動機LMDE(推力變比10∶1)[3]、實現(xiàn)月球表面軟著陸的嫦娥三號探測器主發(fā)動機(推力變比5∶1)[4,5]、多次成功降落回收的獵鷹9 號火箭梅林1D 發(fā)動機(推力變比2∶1)[6].針栓式噴注器的噴霧特性關系著燃燒效率以及發(fā)動機工作可靠性[7].在液體火箭發(fā)動機的研制過程中,噴注器通常是所需周期最長的一個部件,其方案選擇和結構設計也往往是應首先解決的關鍵技術.

針栓式噴注器已有半個多世紀的研發(fā)歷史,但公開文獻中關于噴霧機理的介紹卻非常有限[8-11],尤其是徑向孔型氣液針栓式噴注器.目前該類型噴注器的射流破碎霧化過程的研究以數(shù)值模擬為主,試驗觀測僅作為輔助驗證.Radhakrishnan 等[12,13]借助商用軟件Fluent 提供的歐拉-拉格朗日體系描述徑向液體射流在軸向氣體來流作用下的二維破碎霧化過程,一次破碎和二次破碎過程分別用簡單的Single 噴注模型和Wave 破碎模型(Kelvin-Helmholtz 不穩(wěn)定主導)代替;Son 等[14]則利用Fluent開展歐拉多相流模型的二維軸對稱仿真,該方法對液滴的捕捉能力受限于網(wǎng)格密度.Radhakrishnan和Son 的仿真結果僅能在針栓噴嘴的噴霧形態(tài)和噴霧錐角方面與試驗結果吻合,無法模擬真實的氣液作用界面和表面波發(fā)展過程.張彬[15,16]將Fluent的新功能VOF-to-DPM 模型結合三維網(wǎng)格自適應技術,成功捕捉到了在氣膜作用下單股液體射流的兩種破碎過程:液體射流迎風面和背風面之間的壓差促使射流往下游彎曲,同時低密度氣體對高密度液體產(chǎn)生的加速作用(Rayleigh-Taylor(R-T)不穩(wěn)定性)導致射流迎風面出現(xiàn)表面波,其振幅隨時間發(fā)展最終在波谷處發(fā)生斷裂(柱狀破碎);氣流與液體射流的切向速度梯度則導致射流表面出現(xiàn)尺度遠小于表面波的層狀褶皺(Kelvin-Helmholtz(K-H)不穩(wěn)定),隨著褶皺拉伸變薄陸續(xù)有小液滴被氣流剝離(表面破碎).隨著局部動量比增大,柱狀破碎逐漸取代表面破碎開始主導射流的一次破碎過程,表面波長逐漸增大至最大值,同時連續(xù)射流的斷裂位置逐漸遠離噴孔.林森[17]采用兩相流大渦模擬程序[18,19]對單股矩形噴孔射流進行仿真得到了相似的破碎過程,同時研究發(fā)現(xiàn)當噴孔長寬比變大時,氣動力作用隨著氣液接觸面同步增大,射流破碎尺度減小,三維空間內(nèi)液滴分布更為均勻.

徑向孔型氣液針栓式噴注器的射流破碎霧化過程與超/亞聲速氣流中的液體橫向射流較為相似,差異主要在于液體射流穿透深度與氣流來流厚度的相對大小,而研究人員針對后者已開展了大量工作,包括一次破碎模式和二次破碎模式的劃分[20-22]和射流與液滴表面波產(chǎn)生機理與發(fā)展規(guī)律[19,23-25],結果表明液體射流的破碎霧化主要涉及R-T 不穩(wěn)定和K-H 不穩(wěn)定這兩類機制,為前者的研究提供了參照.試驗觀測手段包括背景光成像法[8,24]、激光陰影法[26]和全息成像法[22],其中背景光成像的應用最為普及;捕捉氣液界面的數(shù)值仿真方法包括體積分數(shù)法VOF (volume of fluid)[14]、水平集法LS (level set)[27]、CLSVOF(coupled LS and VOF)法[17-19,28]和VOF-to-DPM 法[16,29],其中CLSVOF和VOF-to-DPM 的效果更好.

針栓式發(fā)動機工作時噴霧往往集中于近噴孔區(qū)域內(nèi),液體射流從圓孔噴出后受氣膜作用發(fā)生彎曲變形,伴隨著氣液界面不穩(wěn)定波的產(chǎn)生和發(fā)展,連續(xù)射流逐漸破碎成液絲和液滴并迅速蒸發(fā).因此近噴孔區(qū)域存在表面波主導破碎過程(一次破碎)和快速霧化過程(二次破碎),屬于流體力學和兩相流領域的研究難點,也是真正實現(xiàn)氣液針栓式火箭發(fā)動機噴霧與燃燒“耦合計算”的主要問題.由于該區(qū)域復雜的兩相流動耦合將產(chǎn)生大量空間尺度快速變化的特征結構,需要新的思路來克服徑向孔型針栓式噴注器不同液體射流之間的相互干擾及其導致的背景光成像法有效性變差等問題.本文設計了一種采用空氣和水為介質(zhì)的針栓噴注單元作為研究對象,采用兩相流大渦模擬和背景光成像相結合的方法對近噴孔區(qū)域的噴霧場的建立過程及其動態(tài)特性進行仿真與試驗研究.

2 試驗與仿真方案

2.1 氣液針栓噴注單元工況設置

如圖1 左側所示的徑向孔型氣液針栓式發(fā)動機工作時,液體推進劑(Fluid-1)通過針栓中心流道后自針栓頭端的一系列孔呈放射狀徑向噴出,氣體推進劑(Fluid-2)從噴注器集氣腔流經(jīng)針栓外側的環(huán)縫軸向噴出形成氣膜,徑向射流與軸向氣膜呈90°交叉撞擊后推進劑發(fā)生破碎霧化與混合燃燒.為了更好地觀測徑向液體射流與軸向氣膜之間的碰撞過程,設計用于冷試的針栓噴注單元,采用過濾水和干燥空氣來作為模擬介質(zhì).軸向噴注連續(xù)氣膜,噴注面為長方形,徑向孔噴注單股液體射流,將三維噴霧降維處理,增強可視化,以方便分析.

圖1 針栓噴注單元結構示意圖(單位:mm)Fig.1.Conceptual illustration of the gas-liquid pintle injector element (Unit:mm).

4 種不同尺寸的圓形噴孔見圖1 右側,其余主要結構參數(shù)如下:狹縫寬度w取7 mm,狹縫高度h取2.23 mm,跳過距離Ls取12 mm,針栓長度L取22 mm.針栓噴注單元裝配體的具體結構如圖2 所示,其中水流通道的軸向密封采用O 形圈實現(xiàn),集氣腔(蓋板和底座之間)的密封采用紫銅墊片實現(xiàn).

圖2 針栓噴注單元裝配圖Fig.2.Final geometry of the pintle injection element assembly.

噴霧試驗在靜止空氣環(huán)境中開展,環(huán)境溫度取300 K,環(huán)境壓力取101.3 kPa.在節(jié)流條件下空氣和水的質(zhì)量流量分別按如下公式計算:

式中q表示質(zhì)量流量,下標G 和L 分別表示氣體和液體(在本文即為空氣和水);Cd表示流量系數(shù);A表示噴注通道的節(jié)流面積;Pe表示噴出壓力,這里取環(huán)境壓力;Pi表示噴前壓力;T表示流體溫度;R表示氣體常數(shù);k表示比熱比.

冷試工況安排見表1,針對圖1 右側所示的4 種不同直徑的射流噴孔,各自設定的水流質(zhì)量流量與噴孔面積成正比(見表1 第2 列),并與6 種不同的空氣噴注壓降(即Pi與Pe之差,從大到小依次為2.13,0.99,0.62,0.42,0.31,0.24 MPa,與表1第2 行中從左至右的6 種不同空氣質(zhì)量流量相對應)進行組合共形成24 個冷試工況;v表示噴注速度,根據(jù)(1)式和(2)式計算可知vG在不同工況下略有變化,而vL保持33.48 m/s 不變.

表1 冷試工況設置Table 1.Operating conditions of cold tests.

2.2 噴霧試驗與測量系統(tǒng)

冷態(tài)噴霧試驗在常溫常壓下開展,霧化試驗系統(tǒng)如圖3 所示,包括模擬介質(zhì)供應裝置、噴霧收集器、試驗件和壓力流量測控裝置等基本組成部分.試驗過程中空氣直接從高壓氣瓶中噴出,貯箱內(nèi)的水則經(jīng)高壓氮氣增壓后噴出.空氣和水進入針栓噴注單元時的入口壓力和流量由減壓閥和氣動閥控制.噴霧收集器負責將液霧及時抽吸排出以免環(huán)境中液霧彌漫影響觀測結果.

圖3 試驗系統(tǒng)示意圖Fig.3.Schematic diagram of the experimental platform.

噴霧試驗中需測量的常規(guī)參數(shù)包括各工質(zhì)在試驗系統(tǒng)中的貯箱壓力、管路壓力和噴前壓力,以及各路工質(zhì)在管路中的溫度和流量.空氣和水的體積流量測量均采用渦輪流量計,壓力測量采用壓阻式壓力傳感器,溫度測量采用熱電阻.需要注意的是空氣質(zhì)量流量的計算首先要基于流量計附近的溫度和壓力,根據(jù)理想氣體方程得到空氣的密度.另外噴霧場本身不發(fā)光且不易于被照亮,因此非常適合背景光成像進行拍攝.考慮到衡量氣液撞擊過程中連續(xù)液體射流破碎時間的特征時間尺度t*定義見(3)式,表1 所列工況的特征時間尺度t*對應的頻率f均為10 kHz 量級.為了通過高速攝影結果準確識別該特征頻率f,要求采用頻率fs＞ 2f,因此本文中黑白高速相機的拍攝幀頻和曝光時間分別取200 kHz 和1.25 μs.

2.3 兩相流大渦模擬

在超/亞聲速氣流中射流一次破碎和液滴二次破碎的研究中,肖鋒等[18,19,25]提出了一種兩相流大渦模擬算法,采用不可壓流求解器和可壓流求解器來分別計算液體和氣體:采用有限體積法求解不可壓流控制方程,時間離散和空間離散分別通過一階正向投影方法和二階中心差分方法實現(xiàn);采用有限差分方法求解可壓流控制方程,其中時間離散通過二階TVD(total variation diminishing)和Runge-Kutta 方法實現(xiàn),時間步長由取值0.4 的CFL (courant-friedrichs-lewy)數(shù)決定,空間離散采用五階加權WENO(weighted essentially non-oscillatory)方法和二階中心差分方法分別離散對流項和黏性項.構建笛卡爾網(wǎng)格進行求解時,不可壓求解器中物理量交錯分布(速度分量位于網(wǎng)格面,其他變量位于網(wǎng)格中心),而可壓求解器中物理量均位于網(wǎng)格拐角;氣相為不可壓流求解器提供壓力邊界條件,液相為可壓流求解器提供速度邊界條件.另外該算法采用CLSVOF 法追蹤氣/液相界面,綜合了LS 和VOF 的優(yōu)點,既保證了相界面的連續(xù)性,又保證了質(zhì)量守恒,已通過經(jīng)典的液體射流一次破碎試驗結果[22,30]驗證了算法的準確性.

針栓噴注單元的水射流在空氣膜撞擊下的一次破碎過程與超/亞聲速氣流中的橫向射流相類似,本文將采用肖鋒等[18,19,25]提出的兩相流大渦模擬算法對近噴孔區(qū)域圓孔液體射流的表面波主導破碎過程進行仿真研究.本文以圖1 和圖2 為參照構建了如圖4 所示的長方體計算域,其左側平面為包含射流噴孔出口面的基底面(設為壁面邊界條件),上側平面的左端為高2.23 mm 的氣膜狹縫出口面(設為空氣入口面,對狹縫長度不作限制以利于計算域分區(qū)),上側平面的其余部分為腔蓋面(設為壁面邊界條件),計算域的其余外平面均設為出口邊界條件.紅色的液體射流從圓孔噴出后在氣膜作用破碎霧化,當液滴尺寸小于所在位置的網(wǎng)格即消失不見,計算過程中以水射流的初始方向為x軸,空氣射流的初始方向為y軸,展向為z軸.整個計算域尺寸(x×y×z)為50 mm × 100 mm ×40 mm,被劃分為網(wǎng)格數(shù)量相等的256 個塊后可在256 個核上實現(xiàn)并行計算;為降低計算成本,針對射流發(fā)生一次破碎的近噴孔區(qū)域進行加密處理(即通過劃分小尺寸塊以提高網(wǎng)格密度).

圖4 計算域劃分示意圖Fig.4.Schematic diagram of the computational domain division.

3 結果分析與討論

3.1 表面波主導破碎過程

近噴孔區(qū)域液體橫向射流的破碎霧化是氣動力、液體黏性力和表面張力的相互作用結果,其中氣動力促使液體表面擾動的生成與發(fā)展,液體黏性力對液體表面擾動的發(fā)展有阻尼作用,而表面張力則有利于維持液體的聚集狀態(tài).因此破碎霧化類型通常基于韋伯數(shù)(見(4)式)和Ohnesorge 數(shù)(見(5)式)進行劃分[22],前者表示氣動力與表面張力之比,后者表示液體黏性力與表面張力之比.本節(jié)針對表1 中的工況CT14#研究圓孔射流在氣膜作用下的表面波主導破碎過程,并分析破碎霧化典型特征和流場渦結構,該工況的韋伯數(shù)為4350,Ohnesorge數(shù)為0.00324,因此一次破碎過程形態(tài)上均屬于剪切破碎.

式中,d表示圓孔直徑,σ為表面張力,μL為液體黏度.

仿真過程中在噴孔C(d=0.97 mm)出口面指定具有33.48 m/s 法向初始速度的水射流(qL=24.54 g/s),在狹縫出口面指定具有265.22 m/s 法向初始速度的空氣射流(qG=33.41 g/s).參照文獻[17,19]以加密區(qū)網(wǎng)格尺度為25 μm(約0.026d)的計算網(wǎng)格為基準(總網(wǎng)格量在4800 萬左右),同時額外設置兩套對比網(wǎng)格A 和B,加密區(qū)網(wǎng)格尺度分別為35 μm 和50 μm.

文獻[24]提出“噴霧分數(shù)”概念以表示任一空間點(以像素點等價代替)被液體射流/噴霧占據(jù)的概率值,定義為

其中gi代表第i個樣本圖像中待計算空間點處的灰度值,n代表樣本圖像總數(shù)(取值越大越接近真實結果).通過對工況CT14#的1500 幅連續(xù)噴霧圖像進行統(tǒng)計,提取出如圖5(a)所示的近噴孔區(qū)域的時均射流/液霧邊界帶,邊界帶右上方和左下方分別表示恒氣相區(qū)(γ=0)和恒液霧區(qū)(γ=1),噴霧分數(shù)介于0 到1 之間的邊界帶對應射流/噴霧邊界的振蕩,其中包含若干條噴霧分數(shù)等值線,任意一條噴霧分數(shù)等值線均可被視為是液體射流/噴霧的時均邊界.本文將網(wǎng)格加密區(qū)內(nèi)若干時刻的射流/液霧分布結果(通過紅色的氣液界面表示)疊加得到其射流/液霧穿透深度,對比圖5(b)—(d)可以發(fā)現(xiàn)計算網(wǎng)格和對比網(wǎng)格A 的結果與噴霧分數(shù)γ=0.1 等值線表示的噴霧場外邊界基本吻合(誤差將在下文進一步解釋),而對比網(wǎng)格B 的結果與試驗結果的差異較為顯著,驗證了計算網(wǎng)格、對比網(wǎng)格A 和該計算方法對氣液針栓式噴注器一次破碎過程仿真的適用性.這里采用γ=0.1 等值線而非γ=0 等值線,是因為后者受小尺度液滴影響在邊界帶中不太穩(wěn)定,不適合代表真正的噴霧分布范圍.如無特別說明,本節(jié)均采用計算網(wǎng)格以提高對液滴的捕捉能力.

圖5 關于工況CT14#中射流/液霧分布范圍的仿真結果(紅色氣液界面)與試驗數(shù)據(jù)(藍色的噴霧分數(shù)γ=0.1 等值線)對比(a) 時均射流/液霧邊界帶;(b) 計算網(wǎng)格;(c) 對比網(wǎng)格A;(d) 對比網(wǎng)格BFig.5.Comparison of LES-predicted liquid jet/spray distribution (gas-liquid interface colored by red) for case CT14# with experimental data (i.e.the blue spray fraction iso-line of γ=0.1):(a) time-averaged boundary band;(b) computing grid;(c) contrast grid A;(d) contrast grid B.

從圖6 給出的瞬時壓力和速度云圖可以看出,射流上方的空氣壓力和流速相比于氣膜狹縫出口分別出現(xiàn)了顯著的下降和上升,說明氣膜到達噴孔位置前經(jīng)歷了膨脹加速過程,其在x方向的作用范圍也大大超過了狹縫高度(2.23 mm).由于橫向噴出的水射流對高速空氣來流產(chǎn)生阻擋作用,射流前方出現(xiàn)一道較強的脫體弓形激波①.在弓形激波最靠近壁面處,其與空氣來流的夾角接近90°,這是由于水射流離開噴孔后的初始階段基本未發(fā)生變形,其對高速空氣來流的阻擋作用與橫向放置的圓柱相似.高速空氣來臨通過弓形激波后壓力與速度分別劇烈升高和下降,進而在到達氣液交界面時發(fā)生滯止并繞過水射流迎風面往下游流動,氣液動量交換將促使射流發(fā)生彎曲變形;同時繞流運動產(chǎn)生的流動分離現(xiàn)象導致射流背風面出現(xiàn)低壓區(qū),故射流前后壓差也將產(chǎn)生一個垂直射流表面并促進彎曲變形的加速度.

圖6 工況CT14#中瞬時液體射流形態(tài)對應的 (a)壓力云圖與(b)速度云圖Fig.6.(a) Instantaneous pressure and (b) velocity contours of case CT14#,together with the liquid jet structure.

需要注意的是弓形激波并非從壁面生成,弓形激波后側的高壓區(qū)產(chǎn)生逆壓梯度并經(jīng)由邊界層的亞聲速區(qū)往上游傳遞,導致水射流上游的邊界層發(fā)生轉捩和變厚,形成圖6(b)中較為明顯的分離區(qū)②,因此弓形激波是在其上方產(chǎn)生(激波結構只存在于超聲速氣流中).

隨著水射流遠離噴孔,其逐漸往空氣來流方向彎曲且氣液速度差有所減小,而弓形激波由初始的正激波轉變?yōu)樾奔げㄇ覐姸戎饾u減弱,當弓形激波到達右側的高速空氣來流作用邊界時消失.從圖6(a)中還可以發(fā)現(xiàn),在離開噴孔一定距離后,存在個別小激波③從水射流迎風面突起處往外延伸至弓形激波上,小激波生成的位置及強度有一定隨機性,分析認為高速空氣來流經(jīng)弓形激波減速后相比水射流依然保持較大速度差,射流表面的流動受到突起結構的阻礙是形成小激波的原因.

如圖6 所示,連續(xù)射流發(fā)生彎曲變形過程中始終存在一個與其流向相垂直的氣動力,即低密度氣體對高密度液體產(chǎn)生法向加速度,因此水射流將受到R-T 不穩(wěn)定性的影響,表現(xiàn)為迎風面表面波沿水射流流向的發(fā)展[22].圖7 給出了試驗和不同網(wǎng)格仿真得到的典型R-T 不穩(wěn)定表面波沿射流迎風面的發(fā)展規(guī)律,其中試驗結果通過顯化處理得到,即在灰度值線性拉伸的基礎上再通過冪變換進行非線性拉伸,從而凸顯表面波和連續(xù)射流斷裂位置等特征結構,本文采用的冪變換函數(shù)為

圖7 工況CT14#中典型R-T 不穩(wěn)定表面波形態(tài)的試驗[(a)—(b)]與仿真[(c)—(h)]對比:(a) t1 (試驗結果);(b) t1+15 μs (試驗結果);(c) t2 (計算網(wǎng)格仿真結果);(d) t2+15 μs(計算網(wǎng)格仿真結果);(e) t3 (對比網(wǎng)格A 仿真結果);(f) t3+15 μs (對比網(wǎng)格A 仿真結果);(g) t4 (對比網(wǎng)格B 仿真結果);(h) t4+15 μs (對比網(wǎng)格B 仿真結果)Fig.7.Comparison of experimental [(a)—(b)] and simulation [(c)—(h)] results of case CT14# on the distribution of typical R-T unsteady surface waves:(a) t1 (experimental result);(b) t1+15 μs (experimental result);(c) t2 (simulation result of computing grid);(d) t2+15 μs (simulation result of computing grid);(e) t3 (simulation result of contrast grid A);(f) t3+15 μs (simulation result of contrast grid A);(g) t4 (simulation result of contrast grid B);(h) t4+15 μs(simulation result of contrast grid B).

其中gR和gT分別代表冪變換前后圖中的相對灰度值,c為非線性系數(shù)(c＜ 1 時可顯化小灰度值區(qū)域,c＞ 1 時可顯化大灰度值區(qū)域,此處c取10).

從圖7 可以看出,隨著表面波波長與振幅的不斷增長,水射流在展向和流向同時受到拉伸而不斷變薄,結合圖6(b)所示的速度云圖和圖8 所示的渦結構可以認為連續(xù)射流發(fā)生斷裂前波谷處易形成低速漩渦結構(柱狀破碎渦,column breakup vortices),加劇氣液兩相之間的相互作用;當表面張力無法抵抗氣動力時連續(xù)射流將波谷處發(fā)生斷裂,同時在壓差作用下高速氣流自斷裂位置進入,將大量剝離出的液滴往壁面附近輸運從而形成拉絲現(xiàn)象,圍繞液滴和液絲生成的渦將對二次破碎和混合特性產(chǎn)生重要影響.另外噴孔附近還捕捉到由未發(fā)生變形液束與壁面邊界層之間的相互作用激發(fā)形成的馬蹄渦(horseshoe vortices).總的來說,氣液針栓式噴注器的表面波主導破碎過程與超聲速氣流中液體橫向射流的破碎過程有一定相似性[24].

圖8 工況CT14#中瞬時射流結構(紅色)與渦結構(綠色)的疊加圖,其中渦結構由速度梯度第二不變量[31]的等值面表示Fig.8.Instantaneous liquid jet structure (colored by red) of case CT14#,superimposed with the vertical structures identified by isosurface of the second invariance of the velocity gradient tensor (colored by green).

圖7 中通過對比表面波波長還可發(fā)現(xiàn):計算網(wǎng)格和對比網(wǎng)格A 仿真得到的表面波形態(tài)與試驗結果基本相同,同時時間間隔較短的情況下前后兩幅圖像具有明顯的時間相關性(即具有相同的特征結構),有利于對表面波主導破碎過程中特征結構的時間演化特性和速度變化規(guī)律開展研究;而對比網(wǎng)格B 仿真得到的表面波形態(tài)與試驗結果差異較大且時間相關性較差,因此該網(wǎng)格無法滿足計算精度.圖7 中計算網(wǎng)格和對比網(wǎng)格A 的仿真結果相對于試驗結果的不足主要在于仿真無法得到沿連續(xù)射流分布的濃密液霧,原因為仿真計算難以捕捉小于網(wǎng)格尺寸的液滴顆粒,特別是氣液切向速度梯度誘發(fā)K-H 不穩(wěn)定而剝離的大量小尺寸液滴;噴孔附近較為濃密的液霧是由于亞聲速流動的分離區(qū)內(nèi)存在旋轉方向相反的分離渦與再附渦,能將小尺寸液滴沿分離區(qū)往上游傳遞,但在仿真結果中同樣不是很明顯.

3.2 破碎霧化過程的動態(tài)特性分析

如前所述,表面波主導破碎過程對氣液針栓噴注單元的噴霧特性有重要影響,本節(jié)選取典型的工況CT17#進一步分析近噴孔區(qū)域的破碎霧化過程的動態(tài)特性.對工況CT17#的連續(xù)噴霧圖像采用POD 方法(原理介紹參見文獻[7])分析其擬序結構,圖9(a)和圖9(b)所示為其連續(xù)噴霧的某時刻瞬態(tài)圖像和時均結果,圖9(c)—(h)為表征噴霧流場的特征結構的若干POD 模態(tài)(對應時間系數(shù)的功率譜見圖10),其中模態(tài)云圖幅值越大(即顏色越深)說明該處的動態(tài)特征越強,同時比其余幅值較低處對于整體振蕩的貢獻越大.

圖9 工況CT17#的(a)瞬時噴霧圖像,(b)時均結果及(c)—(h)若干POD 模態(tài)Fig.9.(a) Spray snapshot,(b) time average and (c)—(h) several POD modes of case CT17#.

圖10 工況CT17#若干POD 模態(tài)的時間系數(shù)功率譜疊加圖Fig.10.Superposition of the power spectrum densities from several POD modes of case CT17#.

從圖9 中的POD 模態(tài)可以看出工況CT17#的噴霧振蕩主要包含兩類特征:模態(tài)1 和模態(tài)2 中沿噴霧外邊界法向出現(xiàn)幅值及其符號的劇烈變化,表明噴霧流場存在整體擴張/收縮過程,雖然相應功率譜的幅值較大,但主頻在1 kHz 以下,主要受上游管路振蕩引起的氣液介質(zhì)流量變化影響;模態(tài)3 和模態(tài)4 沿噴霧迎風面呈現(xiàn)出比較規(guī)則的空間周期性結構,主要源于連續(xù)射流斷裂后形成的液塊或液霧團在氣流作用下往下游運動,這兩個模態(tài)的功率譜基本一致,無明顯主頻,但在3—6 kHz范圍內(nèi)有較大能量.后續(xù)的模態(tài)10,基本可認為是前4 個模態(tài)的組合或高階諧振(包括中低頻和高頻等多個主頻),未提供值得關注的新特征;而模態(tài)100 在空間結構和功率譜均趨于均勻化,類似于高速攝影采樣頻率下的噪聲.另外前4 個模態(tài)噴霧振蕩的能量貢獻(energy contribution,EC)較大并快速下降,后續(xù)模態(tài)的能量貢獻則處于緩慢減小的趨勢.

模態(tài)3 和模態(tài)4 的空間結構與撞擊式噴嘴研究中的“撞擊波”概念[32]類似,即當韋伯數(shù)超過臨界值后,將從撞擊點輻射出撞擊波.從圖10 已知模態(tài)3 和模態(tài)4 的時間系數(shù)功率譜基本一致,通過交叉功率譜CPSD(cross-power spectrum density,相關介紹參見文獻[33])進一步分析兩者時間系數(shù)的相關聯(lián),圖11 的結果顯示,在交叉功率譜的能量較大的頻段(3—6 kHz)兩者的相位差為90°左右,且圖9(e)和圖9(f)顯示這兩個模態(tài)空間結構相差四分之一波長,故可認為兩者相互耦合形成了正弦或余弦狀的行波結構.通常認為單個行波由兩個駐波疊加而成,模態(tài)3 和模態(tài)4(以Φ3和Φ4表示,相應時間系數(shù)為a3和a4)疊加后形成的行波可近似為

圖11 工況CT17#的POD 模態(tài)3 和4 時間系數(shù)的交叉功率譜幅值和相位角Fig.11.Amplitude and phase angle from CPSD (Mode-3 and Mode-4) of case CT17#.

圖12 展示了工況CT05#,CT11#和CT23#中耦合產(chǎn)生行波結構的POD 模態(tài),其中工況CT23#中相關模態(tài)的位置有所變化,可認為該行波結構在橫向射流中廣泛存在.從POD 模態(tài)中提取噴霧迎風面行波結構的波長可避免從單個噴霧圖像計算波長產(chǎn)生的誤差.如圖13 所示,表1 中不同工況下無量綱行波波長(λ/d,波長與孔徑的比值)和韋伯數(shù)之間呈冪次律關系,經(jīng)擬合得到吻合度良好(R2=96.98%)的如下關系式:

圖12 工況CT05#,CT11#和CT23#中耦合產(chǎn)生行波結構的POD 模態(tài)Fig.12.POD modes that generate the traveling wave structure in case CT05#,CT11# and CT23#.

圖13 無量綱行波波長與韋伯數(shù)變化的關系Fig.13.Variation of dimensionless surface wavelength versus Weber number.

在R-T 不穩(wěn)定表面波主導的一次破碎過程中,對應最大增長率的表面波波長λR-T的計算公式為

式中,σ和ρL分別為液體表面張力和密度,a為液體射流的加速度,

其中CD為有效阻力系數(shù),根據(jù)文獻[19]對亞聲速氣流中液體橫向射流一次破碎過程大渦模擬研究可知CD∝We-0.1.將加速度a的表達式代入(10)式,同時考慮到uG遠大于uL,(4)式中uL的影響可忽略,無量綱表面波波長(λR-T/d)的表達式簡化后與We的關系呈—0.45 冪次的規(guī)律:

因此可認為(9)式所代表的一次破碎后液塊或液霧團在近場區(qū)域迎風面的波動是連續(xù)液體射流斷裂前的R-T 不穩(wěn)定表面波的延續(xù).

最后通過瞬時噴霧圖像重構來說明POD 方法的有效性.假設兩張圖像(xp和xq)的像素分辨率均為m×n,兩者之間的乘積(xp·xq)定義如下:

由前M個模態(tài)重構K個時刻的噴霧圖像,重構誤差εM如下所示:

式中,Ut為t時刻的瞬時噴霧圖像,Φi為第i個模態(tài),ai,t為模態(tài)i的在t時刻的時間系數(shù)(需要注意的是Φ0為K個時刻噴霧圖像的時均結果,a0,t始終為1).

針對圖9(a)所示的瞬時噴霧圖像,采用前M個模態(tài)進行圖像重構,重構結果及其誤差見圖14.在圖9(b)所示的時均結果基礎上,前20 個模態(tài)的疊加結果主要重構出了最顯著的結構特征,重構誤差已降至0.5 以下;隨著參與重構的模態(tài)數(shù)量增至200,連續(xù)射流表面波和一次破碎后液塊或液霧團的波動結構逐漸顯現(xiàn),重構誤差為0.15 左右,其下降速率逐漸減弱;而當M取1000 時已重構得到較為清晰的液滴分布,證明POD 方法應用于噴霧圖像模態(tài)分析的有效性,此時重構誤差僅0.03 左右,但其下降速率進一步減小,也說明單個模態(tài)的影響已微乎其微.圖15 則展示了工況CT05#,CT11#,CT17#和CT23#中瞬時噴霧圖像重構誤差的統(tǒng)計結果,隨著參與圖像重構的模態(tài)數(shù)量的增加,每個工況中重構誤差的減小趨勢相同,同時相同空氣流量下通過改變孔徑增大液體流量將導致重構誤差增大.

圖14 工況CT17#的瞬時噴霧圖像重構Fig.14.Image reconstruction of spray snapshot from case CT17#.

圖15 部分工況中重構誤差的統(tǒng)計結果Fig.15.Statistical results of the reconstruction deviation in several cases.

4 結論

本文以空氣和水為模擬介質(zhì),通過兩相流大渦模擬和噴霧背景光成像,分析了橫向氣膜作用下液體射流在近噴孔區(qū)域的破碎霧化過程及其動態(tài)特性,具體得到如下結論.

1)近噴孔區(qū)域的氣液流場建立過程為:亞聲速氣流離開狹縫后膨脹加速為超聲速氣流,液體射流垂直進入氣體來流,射流上游產(chǎn)生脫體弓形激波.液體射流在氣體來流的作用下逐漸向下游彎曲,相應的弓形激波由初始的正激波轉變?yōu)樾奔げㄇ覐姸戎饾u減弱.而射流迎風面在壓差作用下出現(xiàn)R-T不穩(wěn)定表面波,表面波的發(fā)展導致連續(xù)射流在波谷處發(fā)生斷裂,氣流自斷裂位置進入后將大量剝離出的液滴往壁面附近輸運從而形成拉絲現(xiàn)象.

2) POD 方法可有效重構冷態(tài)試驗中獲得的瞬時噴霧圖像,POD 模態(tài)表明近噴孔區(qū)域的噴霧振蕩主要由噴霧場的整體擴張/收縮過程(低頻)和液塊或者液霧團在迎風面的運動構成,其中后者是具有高頻振蕩的行波結構.無量綱行波波長和R-T不穩(wěn)定表面波波長與韋伯數(shù)之間的冪次律關系十分相近,可認為該行波結構受連續(xù)液體射流斷裂前的R-T 不穩(wěn)定影響而產(chǎn)生.