李佳楠，雷凡培，楊岸龍，周立新

1. 西安航天動力研究所 液體火箭發(fā)動機(jī)技術(shù)重點實驗室，西安 710100 2. 中國船舶集團(tuán)有限公司，北京 100097

燃燒不穩(wěn)定是液體火箭發(fā)動機(jī)研制過程中經(jīng)常遇到的重大技術(shù)問題[1]，其重要特征表現(xiàn)為燃燒室中與聲學(xué)固有振型一致的大幅值背壓振蕩，振幅達(dá)到穩(wěn)態(tài)室壓的10%～100%，振蕩頻率從數(shù)百Hz到15 000 Hz以上[2]。燃燒不穩(wěn)定會導(dǎo)致發(fā)動機(jī)極高加速度的機(jī)械振動，對發(fā)動機(jī)的機(jī)械系統(tǒng)造成破壞，還會破壞推力室內(nèi)的冷卻系統(tǒng)，導(dǎo)致超高熱流流向噴注器以及壁面，引起內(nèi)壁面的燒蝕，對發(fā)動機(jī)造成毀滅性打擊。由于燃燒不穩(wěn)定自身復(fù)雜的耦合特性，以及推力室高溫、高壓、高熱流難以獲取數(shù)據(jù)的惡劣環(huán)境，人們對燃燒不穩(wěn)定的激發(fā)、維持機(jī)理仍沒有清晰的認(rèn)識，噴霧燃燒的各個物理化學(xué)子過程都有可能成為燃燒不穩(wěn)定激發(fā)與維持的關(guān)鍵子過程。

研究燃燒不穩(wěn)定的方法有自發(fā)激勵與強(qiáng)迫響應(yīng)2種方式[3]，自發(fā)激勵必須考慮燃燒室中所有的物理化學(xué)過程，強(qiáng)迫響應(yīng)則是通過在邊界施加擾動，研究強(qiáng)迫擾動下的噴霧燃燒過程。自發(fā)激勵的優(yōu)點是完全揭示了燃燒火焰與背壓振蕩之間的激發(fā)與反饋機(jī)理，強(qiáng)迫響應(yīng)的優(yōu)點則在于可以研究各種不同的子過程對于不同振幅、不同頻率的聲學(xué)壓力振蕩的響應(yīng)。相關(guān)研究表明霧化不僅影響穩(wěn)態(tài)的燃燒過程，而且與燃燒不穩(wěn)定的激發(fā)、維持有著千絲萬縷的聯(lián)系。撞擊式噴嘴的Hewitt準(zhǔn)則[4-6]表明，撞擊式噴嘴燃燒穩(wěn)定性由2個重要參數(shù)控制，一是射流直徑，二是噴射速度。通過改變噴射速度可以改變發(fā)動機(jī)的穩(wěn)定性閾值，實現(xiàn)發(fā)動機(jī)由穩(wěn)定向不穩(wěn)定的轉(zhuǎn)變，Hewitt準(zhǔn)則預(yù)測的穩(wěn)定性邊界與多種全尺寸發(fā)動機(jī)的實際工作特性十分吻合[6]。Qin等[7]通過數(shù)值計算發(fā)現(xiàn)，液滴粒徑分布與空間分布對燃燒不穩(wěn)定產(chǎn)生重要影響，局部聚集的細(xì)小液滴會形成局部準(zhǔn)定容燃燒過程，并在推力室空間內(nèi)形成壓力尖峰，空間隨機(jī)產(chǎn)生的壓力尖峰最終會誘導(dǎo)出聲學(xué)振型，產(chǎn)生燃燒不穩(wěn)定，在這個過程中，噴霧是最為關(guān)鍵的環(huán)節(jié)。Bai等[8]發(fā)現(xiàn)，噴霧自激振蕩必然引起火焰振蕩以及燃料噴前壓力的振蕩，當(dāng)自激振蕩的頻率與燃燒室的固有頻率一致時，會誘發(fā)燃燒室的壓力振蕩，表明霧化的自激振蕩是發(fā)生不穩(wěn)定燃燒的非常關(guān)鍵的誘因。本文將采用強(qiáng)迫響應(yīng)的方法研究強(qiáng)迫擾動條件下的射流撞擊霧化特性。

為了揭示燃燒不穩(wěn)定復(fù)雜的耦合機(jī)理，國內(nèi)外的研究機(jī)構(gòu)采用試驗或者數(shù)值模擬的方法開展了聲學(xué)振蕩下的噴霧或者燃燒特性研究。Baillot等[9]試驗研究了橫向聲場作用下的氣液同軸剪切噴嘴的霧化特性，當(dāng)射流位于速度波腹位置時，在非線性聲輻射壓力的作用下，射流由圓柱液體變?yōu)楸馄揭耗?，液核長度減小，并且破碎更加劇烈。當(dāng)推力室內(nèi)多個噴嘴的霧化特性同時改變時，非穩(wěn)態(tài)燃燒過程會發(fā)生顯著變化，就有可能激發(fā)出燃燒不穩(wěn)定。Carpentier等[10]在Kundt管中研究了橫向聲學(xué)壓力振蕩對于單束射流的影響，當(dāng)射流位于壓力波腹位置時，射流的形態(tài)幾乎沒有發(fā)生變化。而當(dāng)射流位于速度波腹位置時，與Baillot等[9]觀察到的現(xiàn)象類似，射流變?yōu)楸馄揭耗?。?dāng)射流位于壓力波腹與速度波腹之間時，射流在非對稱聲輻射壓力的作用下發(fā)生偏轉(zhuǎn)。Ficuciello等[11]開展了橫向聲學(xué)振蕩條件下氣液同軸射流的霧化試驗，進(jìn)一步驗證了聲學(xué)振蕩對于霧場的直接作用，發(fā)現(xiàn)了聲學(xué)耦合作用下的液滴聚集現(xiàn)象，并基于非線性聲學(xué)理論對試驗觀察到的現(xiàn)象進(jìn)行了闡釋。Hardi等[12]開展了低溫液氧射流試驗，當(dāng)存在橫向聲學(xué)振蕩時，液氧的液核長度減小了近70%。為進(jìn)一步揭示聲波對于霧化過程的影響，基于氣液兩相流大渦模擬算法，Rutard等[13]開展了亞臨界條件下氣液同軸剪切噴嘴在橫向振蕩壓力場下的霧化特性數(shù)值模擬，同樣發(fā)現(xiàn)了射流扁平以及液核縮短的現(xiàn)象，并分析了橫向聲學(xué)振蕩對于噴霧動力學(xué)以及液滴粒徑的影響。Dighe和Gadgil[14-16]研究了橫向聲波對于射流撞擊形成液膜破碎的影響，結(jié)果表明垂直于液膜的聲學(xué)壓力振蕩促進(jìn)了液膜的破碎，生成的液滴數(shù)目增多，液滴平均粒徑減小。Hakim等[17]采用大渦模擬研究了橫向聲學(xué)擾動條件下低溫推進(jìn)劑的火焰特性，當(dāng)火焰位于速度波腹位置時，火焰長度也會減小，并且變成扁平，同時在速度波腹位置出現(xiàn)周期性振蕩的現(xiàn)象，這與橫向壓力振蕩條件下觀察到的冷態(tài)噴霧場有些類似。

筆者團(tuán)隊[18]總結(jié)出了背壓振蕩影響撞擊霧化的兩種作用機(jī)制，背壓振蕩除了對噴霧燃燒產(chǎn)生直接作用外，還會通過改變噴注壓降影響噴射過程，進(jìn)而對后續(xù)的霧化、燃燒等子過程產(chǎn)生影響。當(dāng)推力室中出現(xiàn)縱向的不穩(wěn)定或者噴嘴位于橫向不穩(wěn)定的壓力波腹位置時，這種機(jī)制將會成為背壓振蕩影響霧化的主要作用機(jī)制，這也是本文重點關(guān)注的問題。當(dāng)噴注壓降周期性變化時，噴射速度也會發(fā)生周期性改變，此時，運動速度較快的流體就會追擊運動速度較慢的流體，從而使流體在空間呈現(xiàn)出堆積的現(xiàn)象，這就是Klystron效應(yīng)[19-20]。由系統(tǒng)動力學(xué)分析[21-22]可以證明，燃燒室的壓力振蕩與管路中的壓力振蕩都會引起噴射速度的周期性變化，二者在幅值與相位上存在定量關(guān)系，在一定程度上可以等效處理，由于后端擾動的研究開展較為困難，可以開展前端擾動下的霧化試驗。

Crane等[23-24]設(shè)計了機(jī)械擾動裝置，研究了前端壓力擾動對射流破碎的影響，建立了擾動頻率、幅值與射流破碎長度的關(guān)系，對低速射流破碎的Rayleigh分析作了改進(jìn)，發(fā)現(xiàn)在強(qiáng)迫擾動的作用下，低速射流破碎機(jī)理由表面張力主導(dǎo)向Klystron效應(yīng)主導(dǎo)轉(zhuǎn)變。Chigier[25]應(yīng)用壓電式驅(qū)動裝置對射流前端施加擾動，同樣也發(fā)現(xiàn)了Klystron效應(yīng)，射流速度的周期性變化導(dǎo)致射流形態(tài)出現(xiàn)圓盤狀的結(jié)構(gòu)，在射流本征頻率附近的擾動能增加液滴尺寸的均勻性?？抵覞萚26]采用水力擾動裝置研究了有/無振蕩2種情況下氣液同軸離心式噴嘴自激振蕩的噴霧形態(tài)。在理論分析方面，Heister等[27]采用邊界元方法(BEM)研究了氣體縱向振蕩條件下的射流過程，縱向擾動會使射流速度發(fā)生周期性變化，周期性擾動的速度會導(dǎo)致射流形成Klystron效應(yīng)，通過不同頻率的擾動速度分析發(fā)現(xiàn)，射流對其本征頻率處的擾動響應(yīng)最強(qiáng)。在數(shù)值模擬方面，Srinivasan等[28-29]基于OpenFOAM，應(yīng)用VOF(Volume of Fluid)界面捕捉方法實現(xiàn)了射流破碎的數(shù)值模擬，通過在射流入口給定正弦的擾動速度，研究了射流平均速度以及速度擾動的頻率、幅值對射流破碎的影響。數(shù)值模擬研究表明，在相同的Strouhal數(shù)下速度擾動振幅對射流破碎產(chǎn)生顯著影響，增加振蕩振幅和頻率都能夠增加射流的噴霧錐角。Yang和Turan[30]基于開源程序Gerris[31-32]，研究了低速與高速條件下，擾動頻率、幅值對于射流霧化的影響。對于低頻與中頻擾動而言，隨著擾動頻率的增大，未受擾動的射流長度將會減小，液滴平均粒徑將會增大。對于高頻擾動而言，霧化特性參數(shù)則變化較小，表明射流對高頻擾動的響應(yīng)不敏感。擾動幅值對于液滴數(shù)目以及液滴粒徑具有重要影響，擾動幅值增大，液滴數(shù)目將會增大，而液滴平均粒徑將會減小。參照Bazarov等[33]設(shè)計的壓力擾動裝置，Yang等[34-35]設(shè)計了前端壓力擾動裝置，通過對噴前壓力施加擾動，試驗研究了噴射流量變化的Klystron效應(yīng)對撞擊霧化的影響，本文將在前人的基礎(chǔ)上繼續(xù)對強(qiáng)迫擾動下的撞擊霧化過程開展系統(tǒng)深入的研究。

對霧化過程的研究，傳統(tǒng)上主要采用試驗的方法開展，而與試驗相比，數(shù)值模擬具有獨特的優(yōu)勢，可以對液膜、液絲的破碎過程進(jìn)行細(xì)致的觀察，捕獲有效的數(shù)據(jù)信息，還可以對試驗現(xiàn)象進(jìn)行深入分析，揭示試驗現(xiàn)象背后所蘊含的物理本質(zhì)，本文將采用數(shù)值模擬與試驗結(jié)合的方法開展研究。霧化過程屬于多相、多尺度的流動過程，需要處理相界面追蹤、自由界面大變形等問題，對網(wǎng)格尺度、時間尺度要求很高，計算量很大。本文采用的數(shù)值模擬方法在處理霧化這種多相、多尺度的流動過程中具有諸多優(yōu)勢：采用樹形自適應(yīng)加密算法對空間進(jìn)行離散，同時結(jié)合分段線性的VOF方法，在實現(xiàn)相界面精確捕捉的同時，與固定網(wǎng)格相比顯著降低了計算量，時間步長也隨網(wǎng)格尺度相應(yīng)調(diào)整，實現(xiàn)了時間尺度自適應(yīng)；由高度函數(shù)結(jié)合連續(xù)表面力模型實現(xiàn)表面張力的精確求解，這在霧化過程的數(shù)值計算中尤為重要。

關(guān)于撞擊式噴嘴的霧化機(jī)理及霧化特性，國內(nèi)外已經(jīng)開展了大量的研究工作，基本掌握了其在穩(wěn)態(tài)條件下的工作特性，但對于其在非穩(wěn)態(tài)條件下的霧化特性目前還知之甚少，導(dǎo)致發(fā)動機(jī)設(shè)計中無法全面準(zhǔn)確把握該種型式噴嘴的噴霧特性。當(dāng)發(fā)動機(jī)在工況調(diào)節(jié)狀態(tài)下工作或出現(xiàn)燃燒不穩(wěn)定性時，因無法掌握其噴霧特性而給性能評估以及故障分析和定位造成嚴(yán)重阻礙。由于燃燒不穩(wěn)定是極其復(fù)雜的物理化學(xué)過程，本文將重點關(guān)注背壓振蕩影響撞擊霧化第一種機(jī)制，采用試驗結(jié)合數(shù)值模擬的方法研究噴射速度的周期性變化對于撞擊霧化特性的影響，這對于進(jìn)一步認(rèn)識霧化在燃燒不穩(wěn)定中所起的作用以及揭示燃燒不穩(wěn)定機(jī)理都將起到積極的推動作用。

1 試驗系統(tǒng)

首先介紹一下試驗系統(tǒng)，試驗系統(tǒng)主要由高壓水路供應(yīng)系統(tǒng)、水力擾動裝置、撞擊式噴嘴、脈動壓力數(shù)采系統(tǒng)、高速攝影系統(tǒng)以及同步觸發(fā)裝置組成，其結(jié)構(gòu)示意圖如圖1所示。液體貯箱內(nèi)的工質(zhì)水經(jīng)過高壓氮氣擠壓之后分為兩路，一路流經(jīng)撞擊式噴嘴完成霧化過程，另一路流經(jīng)水力擾動裝置排出。

水力擾動裝置是產(chǎn)生噴前壓力擾動的核心部件，其結(jié)構(gòu)示意圖如圖2所示，主要由定子與轉(zhuǎn)盤兩部分組成，定子與轉(zhuǎn)盤相互貼合。轉(zhuǎn)盤邊緣沿周向均勻布置一圈小孔，由電機(jī)帶動轉(zhuǎn)動，當(dāng)轉(zhuǎn)盤上的孔與定子的孔重疊時就有流體流經(jīng)擾動裝置，當(dāng)定子上的孔被完全遮擋時，就沒有流體流經(jīng)

圖1 試驗系統(tǒng)結(jié)構(gòu)簡圖Fig.1 Sketch of test system

擾動裝置，由此來產(chǎn)生噴前壓力的擾動，通過調(diào)節(jié)電機(jī)轉(zhuǎn)速來改變壓力擾動的頻率，試驗達(dá)到的最高頻率約為3 563 Hz。

由Kistler脈動壓力傳感器采集噴前的脈動壓力，由XP5壓力傳感器記錄周期性變化的噴注壓降，采樣頻率均設(shè)定為20 480 Hz。由高速攝影對霧場進(jìn)行背光拍攝，設(shè)定高速相機(jī)的拍攝頻率為每秒拍攝20 480張圖片，圖像分辨率為512 pixel×512 pixel。通過同步觸發(fā)裝置實現(xiàn)脈動壓力測量與霧場拍攝的同步，同步測量的原理為脈動壓力傳感器對脈動壓力不間斷采集，由信號發(fā)生器產(chǎn)生一個階躍信號觸發(fā)高速攝影工作。高速攝影觸發(fā)的時間精度為2 μs，DG535信號發(fā)生器的時間精度為ns量級，而數(shù)采系統(tǒng)采樣的時間精度為10 μs，根據(jù)木桶原理整個同步觸發(fā)裝置的時間精度應(yīng)為10 μs，目前試驗所做的最高擾動頻率約為3 563 Hz，10 μs的時間精度等級可以滿足試驗需求。

試驗采用的撞擊式噴嘴結(jié)構(gòu)示意圖如圖3所示，兩束圓柱射流軸線的交點定義為撞擊點，液體工質(zhì)加壓噴射之后在撞擊點處相互撞擊完成霧化

圖3 撞擊式噴嘴結(jié)構(gòu)示意圖Fig.3 Sketch of impinging jet injector

過程。圖3中標(biāo)注了撞擊式噴嘴的部分結(jié)構(gòu)參數(shù)，選用了射流直徑d為1 mm與0.8 mm兩種孔徑的噴嘴，撞擊角2θ=60°，噴嘴出口到撞擊點的噴射距離L約為6.9 mm。

2 數(shù)值模擬方案及算例驗證

2.1 數(shù)值模擬方案

基于開源程序Gerris[31-32]實現(xiàn)兩股射流撞擊霧化過程的數(shù)值模擬，Gerris已在多種不同結(jié)構(gòu)形式的噴嘴霧化[36-38]中得到應(yīng)用，將采用的數(shù)值模擬方案描述如下。由八叉樹結(jié)構(gòu)形式的笛卡兒網(wǎng)格對空間進(jìn)行離散，結(jié)合自適應(yīng)加密算法實現(xiàn)特定區(qū)域的局部加密。由有限體積法直接數(shù)值求解不可壓Navier-Stokes方程組，由分段線性的VOF方法對氣液界面進(jìn)行重構(gòu)。表面張力采用Brackbill等[39]提出的CSF(Continuum Surface Force)方法，將相界面上一定寬度區(qū)域的表面張力等效為連續(xù)的體積力，在動量方程中加入這一項。湍流的處理采用隱式大渦模擬方法[40-41]，由數(shù)值黏性充當(dāng)亞格子黏性模擬亞網(wǎng)格尺度的渦耗散?；贚inux系統(tǒng)，以MPI(Message Passing Interface)為并行平臺實現(xiàn)三維兩相流的并行計算。由于采用了動態(tài)的自適應(yīng)網(wǎng)格，為了提高計算效率，采用了動態(tài)負(fù)載平衡技術(shù)。

圖4 數(shù)值模擬計算域Fig.4 Computational domain of numerical simulations

求解的三維不可壓Navier-Stokes方程組為

(1)

(2)

(3)

采用經(jīng)典的VOF方法[42]對氣液界面進(jìn)行重構(gòu)，在流場中引入體積分?jǐn)?shù)α的概念，表征某一網(wǎng)格內(nèi)第一相流體體積與網(wǎng)格體積之比：α=1表示該網(wǎng)格內(nèi)充滿第一相流體；α=0表示網(wǎng)格內(nèi)充滿第二相流體；0<α<1表示網(wǎng)格內(nèi)存在兩相流體，故存在自由界面。由于體積分?jǐn)?shù)的引入，密度的對流方程式(1)可以等效為體積分?jǐn)?shù)α的對流方程

(4)

流體的密度及黏度系數(shù)由網(wǎng)格內(nèi)2種流體的體積分?jǐn)?shù)加權(quán)計算

ρ=ρlα+ρg(1-α)

(5)

μ=μlα+μg(1-α)

(6)

式中：ρl與ρg分別為液相與氣相的密度；μl與μg分別為液相與氣相的黏度系數(shù)。

VOF方法的一個優(yōu)點是不需要對界面的破裂、融合等現(xiàn)象作特殊的處理，因為這些拓?fù)浣Y(jié)構(gòu)的改變都隱式地包含在VOF方法中[43]。因此，VOF方法非常適合應(yīng)用于存在諸多破碎、聚合現(xiàn)象的霧化過程的計算。

采用經(jīng)典的時間分裂投影方法對控制方程進(jìn)行離散，離散后的方程為

(7)

(8)

(9)

(10)

同時需要求解泊松方程

(11)

2.2 算例驗證

關(guān)于Gerris計算穩(wěn)態(tài)撞擊霧化的算例驗證工作在相關(guān)文獻(xiàn)[36,44]中已經(jīng)開展過，本文將進(jìn)一步驗證建立的數(shù)值模擬方案計算非穩(wěn)態(tài)撞擊霧化過程的有效性。對試驗工況進(jìn)行數(shù)值模擬的關(guān)鍵是獲得噴嘴出口的周期性變化速度，以此作為數(shù)值計算的邊界條件?，F(xiàn)有的流量計都是穩(wěn)態(tài)流量計，無法獲得高頻脈動的噴射速度，將通過流量壓降的關(guān)系計算高頻脈動的噴射速度。驗證算例采用直徑為0.8 mm的噴嘴，由XP5表壓傳感器記錄得到隨時間周期性變化的噴注壓降，由于測量點與噴嘴出口距離較近，忽略流體慣性產(chǎn)生的時滯，并假定噴嘴的流量系數(shù)不隨噴注壓降的變化而變化。由壓降與流量的計算公式(12)得到噴嘴的周期性變化流量。由式(13)換算成噴嘴出口的周期性變化速度u。要得到周期性變化的噴射速度還需要標(biāo)定噴嘴的流量系數(shù)Cd，由式(12)可以看出噴嘴流量的平方與噴注壓降呈線性關(guān)系，可以記錄一系列一一對應(yīng)的流量與壓降，通過曲線擬合得到流量系數(shù)。

(12)

(13)

XP5表壓傳感器記錄得到的周期性變化的噴注壓降如圖5所示，經(jīng)過曲線擬合可近似認(rèn)為噴注壓降為正弦形式的振蕩。噴嘴流量系數(shù)的標(biāo)定曲線如圖6所示，與上述分析一致，流量的平方與噴注壓降近似呈線性關(guān)系，經(jīng)過擬合可以得到噴嘴的流量系數(shù)Cd≈0.555。最終得到的噴嘴出口噴射速度隨時間的變化如圖7所示，噴嘴出口的周期性變化速度可以寫成式(14)的形式，其中平均噴射速度u0=16.1 m/s，速度脈動幅值u′=0.186，擾動頻率f=1 257 Hz。

u=u0(1+u′sin(2πft))

(14)

圖5 壓降信號與擬合曲線Fig.5 Pressure drop signals and fitted line

圖6 流量與壓降的對應(yīng)關(guān)系Fig.6 Relationship between mass flow rate and pressure drop

圖7 周期性變化的噴射速度Fig.7 Periodically varying injection velocity

由上述得到的周期性變化速度作為數(shù)值模擬的邊界條件，基于2.1節(jié)建立的數(shù)值模擬方案開展數(shù)值模擬。首先進(jìn)行網(wǎng)格無關(guān)性驗證工作，分別設(shè)定網(wǎng)格最高加密等級為7級、8級與9級，對應(yīng)的最小網(wǎng)格尺度分別為234 μm、117 μm與58.6 μm，計算得到的撞擊式霧場分別如圖8所示。當(dāng)網(wǎng)格分辨率較低時，數(shù)值模擬捕捉到的液滴尺寸較大，無法準(zhǔn)確捕捉液膜的波動破碎過程，同樣無法辨識霧化的Klystron效應(yīng)。隨著網(wǎng)格分辨率的提高，計算得到的霧場結(jié)構(gòu)越來越精細(xì)。當(dāng)網(wǎng)格最高加密等級為9級的時候，數(shù)值模擬能夠捕捉到霧場中液膜、液絲、液滴等精細(xì)結(jié)構(gòu)，同時也捕捉到了非穩(wěn)態(tài)霧化中周期性液滴群聚集現(xiàn)象。因此最高加密等級為9級可以滿足霧場空間發(fā)展對網(wǎng)格分辨率的要求，這與文獻(xiàn)[45]中網(wǎng)格無關(guān)性驗證的結(jié)論一致，在后續(xù)的數(shù)值模擬中設(shè)定網(wǎng)格的最高加密等級為9級。

計算得到了兩束射流從噴射到撞擊形成液膜，液膜波動破碎形成弓形液絲，弓形液絲進(jìn)一步收縮形成液滴的完整霧化過程，整個過程與Anderson等[46]建立的三步霧化模型完全一致。數(shù)值模擬結(jié)果與試驗結(jié)果的對比如圖9所示，圖9(a)為霧場正面的對比，圖9(b)為霧場側(cè)面的對比。從宏觀上看，計算得到的霧場結(jié)構(gòu)與試驗拍攝的霧場結(jié)構(gòu)是非常接近的。霧場正面形成以撞擊點為頂點的具有一定噴霧角的霧化區(qū)域，霧場側(cè)面液膜在撞擊波的作用下波動破碎，并且在周期性變化的噴射速度作用下調(diào)制出Klystron效應(yīng)，霧場下游出現(xiàn)周期性脫落的大尺度液絲，在3.1節(jié)將重點分析強(qiáng)迫擾動下的撞擊霧化特性。圖10展示了計算得到的非穩(wěn)態(tài)噴霧場結(jié)構(gòu)，根據(jù)試驗以及數(shù)值模擬結(jié)果，其噴霧場可以大致劃分為2個區(qū)域，即一次霧化區(qū)域與二次霧化區(qū)域。一次霧化區(qū)域主要包含液膜破碎的區(qū)域，液膜初步破碎形成了大尺寸的液絲以及大顆粒的液滴，是相對粗糙的霧化區(qū)域。二次霧化區(qū)域則更靠近霧場的下游，液絲與大顆粒液滴在氣動力以及表面張力的作用下進(jìn)一步收縮破碎形成更小尺度的球形液滴，霧化最終形成的液滴尺寸不僅取決于一次霧化形成的大顆粒液滴尺寸，還取決于二次霧化形成的小顆粒液滴尺寸。

圖9 數(shù)值模擬噴霧圖像與試驗的對比Fig.9 Comparison of spray images between simulations and tests

圖10 數(shù)值模擬得到的非穩(wěn)態(tài)噴霧場Fig.10 Unsteady atomization field from numerical simulations

圖11 YZ平面的速度分布Fig.11 Distribution of velocity on plane YZ

霧場沿XY截面與YZ截面的體積分?jǐn)?shù)分布灰度圖分別如圖12(a)和圖12(b)所示，黑色區(qū)域的體積分?jǐn)?shù)為0，是氣相區(qū)域，白色區(qū)域表示體積分?jǐn)?shù)為1，為液相區(qū)域，氣液界面處的體積分?jǐn)?shù)0<α<1。 從兩幅圖中都可以明顯觀察到液膜的波動特性，液膜向下游運動過程中，在Kelvin-Helmholtz不穩(wěn)定的作用下發(fā)生破碎。在圖12(a)中定義了液膜的破碎長度與撞擊波的波長λ，其中液膜破碎長度是一個動態(tài)的概念，液膜完全破碎之后，射流撞擊形成液膜的破碎長度定義為撞擊點到液膜完全破碎位置的距離。從圖12(a)可以看出隨著與撞擊點距離的增加，液膜的波動幅值在逐漸增大，達(dá)到破碎的臨界條件之后就開始破碎形成液絲。由于噴射速度是周期性變化的，速度較快的流體會追擊速度較慢的流體，數(shù)值模擬也捕捉到了射流中流體堆積形成的“鼓包”結(jié)構(gòu)。由圖12(b)可以看出，液膜下游黑色與白色區(qū)域交替出現(xiàn)，這正是由于液膜的波動造成的，并且在越靠近下游的區(qū)域，弓形液絲的厚度也在減小，表明液膜在向下游運動過程中厚度是在不斷減小的，這也是造成液膜不穩(wěn)定性的另外一個因素[47]。在液膜破碎之后的下游區(qū)域，由于霧化的Klystron效應(yīng)出現(xiàn)了大尺度的弓形液體結(jié)構(gòu)。

圖12 不同截面的α灰度圖Fig.12 Grayscale images of α of different slices

圖13 液膜破碎區(qū)域速度分量的分布云圖Fig.13 Contour of velocity component of liquid sheet breakup region

圖14 液膜破碎區(qū)域渦量分量的分布云圖Fig.14 Contour of vortex component of liquid sheet breakup region

數(shù)值模擬統(tǒng)計液滴粒徑的方法是由液滴的體積進(jìn)行換算，由包含在氣液界面內(nèi)的網(wǎng)格體積與網(wǎng)格內(nèi)液相體積分?jǐn)?shù)相乘求和得到液滴的體積，假定液滴形狀為球形，由球形的體積計算公式計算得到液滴的直徑。統(tǒng)計不同網(wǎng)格分辨率條件下霧場液滴尺寸分布的概率密度函數(shù)(Probability Density Function, PDF)，其定義為某個粒徑區(qū)間內(nèi)的液滴概率，不同網(wǎng)格分辨率的液滴粒徑統(tǒng)計結(jié)果如圖15所示。當(dāng)網(wǎng)格加密等級較低時，數(shù)值模擬會捕捉到一些大尺寸液滴。隨著網(wǎng)格分辨率的提高，數(shù)值模擬捕捉到的液滴尺寸開始向小粒徑區(qū)域偏移。需要指出的是，當(dāng)網(wǎng)格分辨率不夠時，就會發(fā)生數(shù)值破碎而產(chǎn)生偽液滴，這種情況會隨著網(wǎng)格分辨率的提高而得到改善，但是并不能夠完全消除。如果非線性不穩(wěn)定也能產(chǎn)生飛濺的小液滴的話，由于網(wǎng)格分辨率的缺陷也會導(dǎo)致對這種小液滴捕捉不夠準(zhǔn)確，因此需要對這種小液滴進(jìn)行剔除。數(shù)值模擬的有效性應(yīng)體現(xiàn)在可以較為完整地展現(xiàn)整個霧化過程，并且數(shù)值破碎產(chǎn)生的小液滴所占的比例足夠小。假定小于4個網(wǎng)格捕捉到的液滴屬于偽液滴，使用Gerris自帶的RemoveDroplets函數(shù)將小于4個網(wǎng)格捕捉到的液滴全部剔除，因此統(tǒng)計曲線的起始點并不是零點。由圖15可以看出，計算得到的液滴粒徑分布范圍大致在0～350 μm，液滴粒徑概率密度分布的峰值在50～200 μm之間，概率密度數(shù)值基本在0.02以上。

圖15 不同網(wǎng)格分辨率的液滴尺寸分布Fig.15 Droplet diameter distributions of different mesh resolutions

3 結(jié)果分析

3.1 強(qiáng)迫擾動下的射流撞擊霧化特性

在前面已經(jīng)提到，背壓振蕩與噴前壓力擾動都會引起噴射速度的周期性變化，二者在幅值與相位上存在一定關(guān)系，可以通過開展噴前壓力擾動試驗，并結(jié)合系統(tǒng)動力學(xué)分析將噴前壓力擾動轉(zhuǎn)化為背壓振蕩對噴射的影響。采用水力擾動裝置產(chǎn)生噴前壓力的擾動，開展了前端擾動條件下的射流撞擊霧化試驗，試驗采用孔徑為1 mm的噴嘴。首先對比強(qiáng)迫擾動霧化與自然霧化之間的差別，然后著重探討強(qiáng)迫擾動條件下的撞擊霧化特性。圖16為脈動壓力傳感器記錄得到的典型噴前脈動壓力p′，水力擾動裝置迫使噴前壓力出現(xiàn)周期性振蕩，從而使噴注壓降發(fā)生周期性改變，最終導(dǎo)致噴射速度周期性變化，并調(diào)制出霧化的Klystron效應(yīng)。對該脈動壓力信號作快速傅里葉變換(FFT)，可以得到該脈動壓力的頻譜圖，如圖17所示，脈動壓力的頻率為1 782 Hz。由高速攝影對存在噴前壓力擾動與不存在噴前擾動的噴霧場分別進(jìn)行拍攝，圖18(a)為不存在擾動時的自然霧化圖像，圖18(b)為存在噴前壓力擾動的霧化圖像。從宏觀上看，二液滴空間分布比較均勻，而對于存在擾動時的噴霧場，在液膜下游則出現(xiàn)了大量液滴群聚集的弓形結(jié)構(gòu)，觀察不同時刻的高速攝影圖像可以發(fā)現(xiàn)弓形液滴群的出現(xiàn)表現(xiàn)出一定的周期性。Rayleigh準(zhǔn)則[48]指出了發(fā)生燃燒不穩(wěn)定的2個必要條件：① 存在周期性的釋熱波動；② 壓力振蕩與釋熱波動同相位。當(dāng)同時滿足這2個條件時，壓力振蕩就會在正反饋機(jī)制下持續(xù)放大。霧化的周期性必然導(dǎo)致周期性的燃燒釋熱，如果周期性的燃燒釋熱與周期性的壓力振蕩耦合在一起形成正反饋，則有可能調(diào)制出燃燒不穩(wěn)定，由此可以推斷，背壓振蕩調(diào)制出的周期性霧化有可能成為燃燒不穩(wěn)定的驅(qū)動機(jī)制之一。

圖16 噴前脈動壓力Fig.16 Pulsating pressures before injection

圖17 脈動壓力的頻譜分布Fig.17 Spectral distribution of pressure fluctuations

圖18 撞擊式噴嘴的噴霧場Fig.18 Atomization field of impinging liquid jets

(15)

為獲得更多的數(shù)據(jù)信息，本文基于開源程序Gerris開展了數(shù)值模擬，重點關(guān)注的是噴射速度的周期性變化對于霧化特性的影響，因此，模擬強(qiáng)迫擾動霧化的方法是按照式(14)在邊界給定隨時間周期性變化的噴射速度。設(shè)定射流的平均速度u0=31.5 m/s，速度的擾動幅值u′=0.05，擾動頻率f=1 782 Hz，其他結(jié)構(gòu)參數(shù)與試驗所用的噴嘴參數(shù)一致。數(shù)值計算得到的自然噴霧場與強(qiáng)迫擾動下的噴霧場對比如圖22所示，二者之間的差別與試驗觀察到的霧場區(qū)別一致，自然噴霧場的液滴空間分布比較均勻，而強(qiáng)迫擾動下的液膜下游出現(xiàn)了大尺度的弓形結(jié)構(gòu)，并且在霧場出現(xiàn)了局部液滴聚集的周期性現(xiàn)象。在一個周期T內(nèi)，霧場隨時間的變化如圖23所示，由一個周期內(nèi)的霧化過程可以對撞擊霧化的Klystron效應(yīng)形成機(jī)制進(jìn)行闡述。當(dāng)噴射速度出現(xiàn)周期性變化時，運動速度較快的流體追擊運動速度較慢的流體，在射流中出現(xiàn)“鼓包”的現(xiàn)象，如圖23(a)，之后“鼓包”向下游運動并在撞擊點位置相互撞擊，如圖23(d)。該現(xiàn)象類似于射流速度更高、直徑更大的兩束射流撞擊，撞擊之后的液膜厚度與寬度增加，如圖23(e)。 液膜變得更加不穩(wěn)定，大尺度的弓形液團(tuán)形成并在液膜下游脫落，如圖23(f)～圖23 (h)。液團(tuán)在向下游運動過程中快速破碎并形成大量在空間聚集的液滴群，整個物理過程在時間上表現(xiàn)出顯著的周期性特征，類似于射流的Klystron效應(yīng)，可以認(rèn)為是撞擊霧化表現(xiàn)出的Klystron效應(yīng)。

圖19 3個觀察窗的大小與位置Fig.19 Positions and sizes of three interrogation windows

圖20 自然霧化隨時間變化的灰度值與相應(yīng)的FFT結(jié)果Fig.20 Time varying grayscale values and corresponding FFT results for natural atomization

圖21 強(qiáng)迫擾動霧化隨時間變化的灰度值與相應(yīng)的FFT結(jié)果Fig.21 Time varying grayscale values and corresponding FFT results for atomization with forced perturbations

圖22 數(shù)值模擬得到的不同噴霧場的對比Fig.22 Comparison of different atomization fields from numerical simulations

圖24 不同位置處速度分量隨時間的變化Fig.24 Time varying velocity component at different positions

與試驗統(tǒng)計霧場信息的方法類似，在撞擊點下游某個位置處設(shè)定觀察窗，統(tǒng)計該區(qū)域的液滴信息。統(tǒng)計撞擊點下游10d與20d位置處，厚度為2d，寬度為30d區(qū)域內(nèi)的液滴數(shù)目以及Sauter平均直徑(SMD)隨時間的變化。自然霧化的液滴信息統(tǒng)計結(jié)果如圖25所示，可以看出液滴數(shù)目以及Sauter平均直徑的變化與試驗觀測的灰度信息類似，呈現(xiàn)出隨時間無規(guī)律變化的特征。而當(dāng)噴射速度周期性變化時，觀察窗內(nèi)的液滴數(shù)目以及液滴的Sauter平均直徑則表現(xiàn)出周期性變化的特征，如圖26所示，并且周期性波動的頻率與施加的強(qiáng)迫擾動的頻率一致。對某個時刻兩個觀察窗內(nèi)的液滴粒徑分布進(jìn)行統(tǒng)計，統(tǒng)計結(jié)果分別如圖27(a)和圖27(b)所示。當(dāng)存在強(qiáng)迫擾動時，小尺寸液滴所占的比例提高，大尺寸液滴所占的比例減小，表明在強(qiáng)迫擾動的作用下，觀察窗內(nèi)的液滴平均尺寸有減小的趨勢。對一段時間內(nèi)觀察窗內(nèi)的液滴數(shù)目與液滴粒徑進(jìn)行時間平均，統(tǒng)計結(jié)果如表1所示。當(dāng)存在周期性擾動時，觀察窗內(nèi)的液滴數(shù)目顯著增多，時間平均之后的Sauter平均直徑顯著減小，最大Sauter平均直徑也相應(yīng)減小。與自然霧化相比，強(qiáng)迫擾動霧化的2個觀察窗內(nèi)液滴數(shù)目分別提高了6.0%與13.3%， 時間平均的Sauter平均直徑分別減小了8.8%與5.5%。由數(shù)值模擬結(jié)果可以看出，當(dāng)存在前端擾動時，霧場空間局部區(qū)域生成的液滴數(shù)目顯著增多，并且液滴粒徑有減小的趨勢。由Qin等[7]的研究成果可以得知，燃燒不穩(wěn)定的激發(fā)與液滴粒徑有強(qiáng)的相關(guān)關(guān)系，液滴粒徑的減小與液滴數(shù)目的增多更有利于激發(fā)燃燒不穩(wěn)定，強(qiáng)迫擾動下的霧化有可能成為燃燒不穩(wěn)定激發(fā)的重要中間環(huán)節(jié)。

圖25 自然霧化觀察窗內(nèi)的液滴數(shù)與SMD隨時間的變化Fig.25 Time varying droplet number and SMD of interrogation windows of natural atomization field

圖26 強(qiáng)迫擾動霧化觀察窗內(nèi)的液滴數(shù)與SMD隨時間的變化Fig.26 Time varying droplet number and SMD of interrogation window of atomization field coupled with forced perturbations

圖27 觀察窗內(nèi)的液滴粒徑分布Fig.27 Droplet diameter distribution in interrogation windows

表1 自然霧化與強(qiáng)迫擾動霧化觀察窗內(nèi)的液滴粒徑信息Table 1 Droplet information of two interrogation windows of natural atomization and forced atomization

3.2 擾動頻率對霧化特性的影響

這一部分將分析擾動頻率對于撞擊霧化的影響，通過改變電機(jī)轉(zhuǎn)速實現(xiàn)對頻率的調(diào)節(jié)，圖28給出了不同頻率的噴前脈動壓力作快速傅里葉變換之后的頻譜圖，擾動的頻率從1 347 Hz一直變化到3 563 Hz。圖29為對應(yīng)的不同擾動頻率條件下的霧場圖像，從圖29可以看出在所研究的頻率范圍(1 347～3 563 Hz)內(nèi)，撞擊式噴嘴的霧化對擾動都存在響應(yīng)，在強(qiáng)迫擾動的作用下撞擊霧場都表現(xiàn)出周期性特征，表明撞擊式噴嘴的霧化在一個非常寬的頻帶范圍內(nèi)對擾動都存在響應(yīng)。這也表明撞擊式噴嘴的工作狀態(tài)易受到外界的影響，并不是一種工作特別穩(wěn)定的噴注霧化單元，這種結(jié)構(gòu)形式的噴嘴非常容易產(chǎn)生燃燒不穩(wěn)定[20]。

圖28 不同脈動壓力的頻譜分布Fig.28 Spectral distribution of different pressure fluctuations

從宏觀上看，隨著擾動頻率的增大，空間局部聚集的弓形液滴群之間的間距在減小，表明液膜下游弓形液絲脫落的頻率也在提高，與強(qiáng)迫擾動的頻率變化保持一致。

圖29 不同頻率下的強(qiáng)迫擾動噴霧場Fig.29 Forced atomization field with different pressure fluctuation frequencies

圖30 同步測量的灰度與脈動壓力信號Fig.30 Synchronous signals of grayscale values and fluctuating pressures

(16)

圖31 壓力振蕩與非穩(wěn)態(tài)釋熱之間的時滯Fig.31 Time delay between pressure fluctuations and unsteady heat release

假定相鄰2個脈動壓力波形之間的相位差為360°，統(tǒng)計得到的不同頻率條件下霧場灰度與脈動壓力之間的相位差如圖32所示。從圖中可以看出擾動頻率對相位差影響很大，霧場相同位置處，擾動頻率越高，灰度信息與擾動壓力之間的相位差值也就越大。在相同的距離范圍內(nèi)，擾動頻率越高，統(tǒng)計得到的相位差的變化區(qū)間也就越大，表明在相同的距離內(nèi)所包含的灰度變化周期更多，這也與從圖29中觀察到的弓形液滴群之間距離減小的現(xiàn)象一致。文中得到的是噴前脈動壓力與霧場灰度之間的相位關(guān)系，通過系統(tǒng)動力學(xué)分析可以得到產(chǎn)生相同流量脈動的條件下，后端擾動與前端擾動之間的相位關(guān)系，這樣就可以間接得到后端壓力擾動與脈動灰度之間的相位關(guān)系，考慮其他子過程之后由圖31可以得到最終的相位關(guān)系。未來可進(jìn)一步考慮燃燒的影響，建立起基于周期性霧化的燃燒穩(wěn)定性分析理論體系。

圖32 不同頻率灰度值與振蕩壓力之間的相位角Fig.32 Phase angle between grayscale values and pressure fluctuations at different frequencies

3.3 擾動幅值對霧化特性的影響

液體火箭發(fā)動機(jī)推力室內(nèi)的室壓通常都存在小幅值的波動，只有當(dāng)發(fā)生燃燒不穩(wěn)定時，擾動壓力的幅值才會超出穩(wěn)態(tài)室壓的10%以上，擾動壓力的幅值是甄別燃燒不穩(wěn)定是否發(fā)生的一個重要參數(shù)，這一部分將采用數(shù)值模擬的方法研究擾動幅值對撞擊霧化的影響。當(dāng)推力室內(nèi)出現(xiàn)縱向振蕩時，噴嘴處于壓力波腹的位置，反壓振蕩影響霧化的主要作用機(jī)制是通過改變噴注壓降繼而改變噴射速度，從而對霧化產(chǎn)生影響，也就是本文所研究的作用機(jī)制。選取與某型推力室一階縱向不穩(wěn)定頻率接近的頻率作為參照頻率來研究擾動幅值對于撞擊霧化的影響。根據(jù)式(14)，設(shè)定擾動頻率為2 230 Hz，擾動幅值的大小u′依次為0.01，0.05，0.1，0.2與0.5，計算得到不同擾動幅值條件下的噴霧場分別如圖33所示。由圖33(a)可以看出，當(dāng)擾動幅值很小時，霧場空間并沒有形成局部聚集的液滴群，強(qiáng)迫擾動的噴霧場與自然噴霧場相比并沒有顯著差別，可以認(rèn)為噴嘴此時工作在準(zhǔn)穩(wěn)態(tài)條件下。只有當(dāng)擾動幅值高于某一臨界值之后，霧場才會發(fā)生顯著改變。由圖33(b)～圖33(e)可以看出，隨著擾動幅值的增大，振蕩壓力場對于霧化的影響加劇，空間局部聚集的液滴群更加明顯，并且液膜的變形更加劇烈。圖34為不同擾動幅值條件下，霧場XY截面體積分?jǐn)?shù)分布灰度圖。從圖中可以看出，擾動幅值的大小對于撞擊前射流的狀態(tài)影響較大。擾動幅值增大，射流中液體堆積形成的“鼓包”更加明顯，表明霧場的Klystron效應(yīng)更加顯著。并且擾動幅值對于射流撞擊之后的液膜波動狀態(tài)、液膜的厚度分布以及破碎長度都有影響，液膜的破碎機(jī)制也發(fā)生變化，破碎過程由撞擊波主導(dǎo)開始轉(zhuǎn)變?yōu)樽矒舨ㄅcKlystron效應(yīng)共同主導(dǎo)。

圖33 不同擾動幅值下的強(qiáng)迫擾動噴霧場Fig.33 Forced atomization field with different perturbation amplitudes

圖35 不同擾動幅值下液膜的無量綱破碎長度Fig.35 Dimensionless breakup lengths of liquid sheet with different perturbation amplitudes

圖36 不同擾動幅值下撞擊點下游3d位置的速度分量Fig.36 Velocity component on point 3d downstream impingement point with different perturbation amplitudes

4 結(jié) 論

燃燒不穩(wěn)定的顯著特征是推力室內(nèi)周期性的背壓振蕩，背壓振蕩會引起噴射速度的周期性變化，從而調(diào)制出與自然霧化顯著不同的霧化特性。本文采用數(shù)值模擬結(jié)合試驗的方法研究了強(qiáng)迫擾動條件下的射流撞擊霧化特性，主要目的就是全面把握噴嘴的工作特性，進(jìn)一步認(rèn)識自然霧化與強(qiáng)迫擾動霧化之間的差別以及霧化在燃燒不穩(wěn)定中所起的作用，得出的主要結(jié)論如下：

1)基于樹形自適應(yīng)網(wǎng)格與分段線性的VOF方法建立的霧化數(shù)值模擬方案可以細(xì)致刻畫出強(qiáng)迫擾動條件下的撞擊霧化過程，模擬得到了強(qiáng)迫擾動作用下的周期性噴霧場，可以為燃燒不穩(wěn)定正反饋機(jī)制的研究提供重要的參考，并為非穩(wěn)態(tài)燃燒的數(shù)值模擬提供準(zhǔn)確的初邊值條件。

2)強(qiáng)迫擾動的噴霧場與自然噴霧場之間存在顯著差別，自然噴霧場的液滴在空間分布比較均勻，強(qiáng)迫擾動條件下的射流撞擊噴霧場出現(xiàn)了弓形液滴群局部聚集的現(xiàn)象，并且在時間上表現(xiàn)出周期性特征，霧化頻率與強(qiáng)迫擾動的頻率一致。在本文研究的頻率范圍(1 257～3 563 Hz)內(nèi)，撞擊式噴嘴的霧化對擾動都有響應(yīng)，表明撞擊式噴嘴的霧化在一個非常寬的頻帶范圍內(nèi)對擾動都存在響應(yīng)。

3)擾動頻率主要影響霧場空間相鄰弓形液滴群之間的間距，以及霧場與振蕩壓力場之間的相位關(guān)系。擾動頻率越高，霧場空間某一固定位置與擾動壓力之間的相位差越大，并且在相同的距離范圍內(nèi)，相位差的變化區(qū)間范圍也越大。

4)擾動幅值決定了霧化Klystron效應(yīng)的強(qiáng)度，只有當(dāng)擾動幅值高于某一臨界值時，霧場才會表現(xiàn)出周期性特征。擾動幅值增大，液膜的破碎長度減小，撞擊點下游的流量波形由正弦波向陡峭前緣波轉(zhuǎn)變，表明流量由線性變化向非線性變化轉(zhuǎn)變。

5)強(qiáng)迫擾動下噴霧場的許多特性都更加有利于激發(fā)燃燒不穩(wěn)定，比如霧場表現(xiàn)出的周期性特征，霧場出現(xiàn)局部聚集的液滴群，以及流量特性由線性向非線性的轉(zhuǎn)變等等，這些特性的出現(xiàn)有可能在正反饋機(jī)制中起到推波助瀾的作用，需要在非穩(wěn)態(tài)霧化的基礎(chǔ)上開展燃燒的研究來進(jìn)一步明確這些特征對燃燒特性的影響。

下一步工作：采用光學(xué)手段直接測量噴嘴出口周期性變化噴射速度，為數(shù)值模擬提供更加準(zhǔn)確的邊界條件；進(jìn)一步測量霧場的微觀信息，采用微距鏡頭結(jié)合圖像處理獲得瞬態(tài)液滴粒徑隨時間的變化規(guī)律；在霧化研究的基礎(chǔ)上進(jìn)一步建立燃燒響應(yīng)模型，從而構(gòu)建起基于周期性霧化的燃燒不穩(wěn)定分析理論體系。

致 謝

本項工作在國家超級計算天津中心的“天河一號”超級計算機(jī)上完成，感謝天津超算中心的大力支持。