王 兵, 范文琦, 徐 勝, 高 瞻

(清華大學(xué)航天航空學(xué)院, 北京 100084)

引 言

航空航天、 機(jī)械工程、 化工生產(chǎn)和生物醫(yī)學(xué)等領(lǐng)域廣泛存在著可壓縮兩相界面流動和反應(yīng)流動過程, 流動過程中往往存在著復(fù)雜的非定常波系： 例如超燃沖壓發(fā)動機(jī)燃燒室噴霧燃燒過程中存在著燃油霧化、 液滴破碎、 液滴撞擊壁面和液滴汽化燃燒等現(xiàn)象, 及其與非定常激波系的強(qiáng)烈相互作用; 渦輪泵誘導(dǎo)輪和螺旋槳葉片等位置產(chǎn)生的空化初生、 氣泡潰滅和空泡射流生成等現(xiàn)象, 以及潰滅誘導(dǎo)的波系等。這些高瞬變兩相流與反應(yīng)流的流場組織方法與調(diào)控手段是工程應(yīng)用的關(guān)鍵技術(shù), 取決于對基本物理過程的機(jī)理認(rèn)知及關(guān)鍵特征的定量表征。

高瞬變兩相流與反應(yīng)流研究中的核心問題是在極高速、 強(qiáng)可壓等極端條件下的質(zhì)量、 動量、 熱量傳遞, 以及化學(xué)反應(yīng)機(jī)理動力學(xué)規(guī)律, 即“三傳一反”機(jī)理的闡述和數(shù)理模型的建立。其難點在于存在激波和相界面等多種強(qiáng)間斷, 以及多種間斷的復(fù)雜作用。激波是一種強(qiáng)擾動的傳播, 一般厚度約為10 μm, 而在無混溶情況下的相界面是不同流體的分界面, 沒有厚度。目前, 盡管人們做了大量的研究工作, 但針對高瞬變兩相流動過程的數(shù)理模型依然不完備, 計算方法精度低, 對高瞬變兩相流與反應(yīng)流機(jī)理與表征方法認(rèn)識不清晰, 對高瞬變兩相流動過程控制不精準(zhǔn), 兩相界面流與反應(yīng)流中非定常波系過程也尚未得到系統(tǒng)研究。

本文基于實驗室近十年的研究成果, 試圖呈現(xiàn)在極端條件下兩相流與反應(yīng)流的研究工作與進(jìn)展。第1章闡釋了高瞬變兩相流的數(shù)理模型與計算方法, 第2章介紹了兩相流體界面約束波系解析方法及運(yùn)用, 第3章針對發(fā)動機(jī)燃燒室內(nèi)的強(qiáng)可壓縮兩相噴霧反應(yīng)流動, 通過數(shù)值模擬研究了斜激波增混與強(qiáng)化燃燒機(jī)理, 探究了可壓縮反應(yīng)流動的多物理過程, 提出了定量表征燃燒模式的第三Damk?hler數(shù)DaⅢ。

1 高瞬變兩相流數(shù)理模型與計算方法

一般采用激波捕捉法或激波裝配法進(jìn)行數(shù)值處理復(fù)雜波系, 采用界面捕捉法或界面追蹤法進(jìn)行數(shù)值處理復(fù)雜相界面。例如, 國內(nèi)外學(xué)者對激波/液滴相互作用[1-5]、 激波/液柱相互作用[6-7]、 空泡潰滅形成射流[8-10]等過程展開了系列研究。數(shù)值格式研究方面, Saurel等[11]系統(tǒng)介紹了針對可壓縮兩相流體流動的擴(kuò)散界面捕捉方法, Allaire等[12]提出了針對可壓縮流體界面仿真的五方程模型。這些都是本文研究的基礎(chǔ)。

1.1 方法的基本思想

高瞬變兩相流與反應(yīng)流中同時存在多種間斷(接觸間斷、 激波、 物質(zhì)界面等), 界面兩側(cè)流體密度比大(如103及以上)、 存在快速相變。通常來講, 在兩相流動研究中, 一般采用不同模型描述物質(zhì)界面間斷和激波間斷, 這樣做, 一是難以實現(xiàn)數(shù)值格式的高精度, 二是容易出現(xiàn)非物理解, 在高瞬變復(fù)雜波系中相變求解過程可能導(dǎo)致數(shù)值不守恒。因此, 實驗室畢業(yè)生張文斌[13]、 項高明[14]和吳汪霞[15]針對該問題建立了一套高瞬變強(qiáng)可壓縮兩相流與反應(yīng)流系統(tǒng)描述的新模型。

針對可壓縮多相或多組分流動, 守恒型Euler方程可表示為

將物質(zhì)界面看作間斷, 可以遷移使用激波捕捉格式捕捉物質(zhì)界面, 因此定義界面量輸運(yùn)方程

式中,φ定義為

同時, 引入剛性氣體狀態(tài)方程, 把兩相變單相

式中,E為內(nèi)能,p為壓強(qiáng),ρ為密度,u為速度。

通過將剛性氣體狀態(tài)方程中表征物理屬性的變量當(dāng)作界面變量, 使得在數(shù)值上可以用單相算法求解兩相流動問題, 將物質(zhì)界面當(dāng)作間斷進(jìn)行求解。

將物質(zhì)界面轉(zhuǎn)化為間斷后, 可用高精度激波捕捉方法同時捕捉激波和物質(zhì)界面等間斷, 實現(xiàn)非定常多波系和相界面的捕捉。因此, 這一做法將高瞬變兩相流動數(shù)值模型問題轉(zhuǎn)化為更加成熟的高精度、 多波系激波捕捉問題。

1.2 復(fù)雜多波系與兩相界面同時捕捉格式

當(dāng)流場中存在間距較小的雙重或多重間斷時, 經(jīng)典WENO格式(WENO-JS)不再適用。Wang等[16]提出了變模板點WENO重構(gòu)格式(WENO-IS), 該格式能夠處理多重間斷, 特別是激波與物質(zhì)界面相互作用時的多重間斷。在該格式中, 在多重間斷處將模板點拆分為多段組合的模板點, 例如對于5點的迎風(fēng)模板, 采用兩個3點模板和兩個兩點模板進(jìn)行重構(gòu), 如圖1所示。在計算xi+1/2半點值時, 迎風(fēng)的5階WENO-JS格式[17]采用3個3點模板進(jìn)行重構(gòu), 而WENO-IS將5階WENO-JS格式中的一個3點模板S12拆分為兩個兩點模板S1和S2。 當(dāng)所有候選模板都穿過間斷(激波或界面)時, WENO-IS 格式仍可以選擇模板S1或S2進(jìn)行計算, 實現(xiàn)多重間斷處理。

圖1 WENO-IS格式模板點示意圖[14]Fig. 1 Stencils of WENO-IS scheme[14]

在xi+1/2半點處值為

式中,ωk為非線性權(quán)重系數(shù), WENO-IS格式中定義為

式中,αk定義為

式中,βk為候選模板的光滑指數(shù),τ5為全局光滑指數(shù),dk為候選模板的最佳權(quán)重,ε為一小量(參考5階WENO-JS格式, 此處將ε取值為10-6)。

參照5階WENO-JS格式, 定義WENO-IS格式中光滑指數(shù)βk為

同時, 定義全局光滑因子τ5為

式中, 系數(shù)cm值見表1。最終, 通過WENO-IS格式改進(jìn)WENO算法的耗散特性, 使格式具有捕捉多間斷的特性。Wang等[16]檢驗了該格式的精度以及收斂性, 此處不再贅述。

表1 系數(shù)cm值

1.3 基于化學(xué)勢弛豫的瞬變相變模型

為了實現(xiàn)兩相流場中高瞬變相變過程的數(shù)值描述, Wu等[18]進(jìn)行了相變空化初生數(shù)值模擬, 給出了判斷條件的壓強(qiáng)函數(shù)。針對化學(xué)勢弛豫過程, 采用了基于相平衡條件的化學(xué)勢松弛平衡算子, 基于熱力學(xué)平衡條件進(jìn)行迭代求解。相變程序的求解過程如圖2所示, 在求解時間步tn時, 若當(dāng)?shù)貭顟B(tài)達(dá)到相變條件, 則進(jìn)入化學(xué)勢弛豫過程進(jìn)行迭代求解, 直到氣相化學(xué)勢μv與液相Gibbs自由能gl小于某極小值ε時, 完成相變過程的求解, 更新當(dāng)?shù)氐臓顟B(tài)變量。

圖2 相變程序求解流程圖[18]Fig. 2 Flow chart of the phase-transition procedure[18]

圖2液相相變判據(jù)中,p為當(dāng)?shù)貕簭?qiáng),pthreshold為壓強(qiáng)閾值, 若當(dāng)?shù)貕簭?qiáng)p小于壓強(qiáng)閾值pthreshold, 則液相發(fā)生空化相變, 因此選擇合理、 合適的空化相變判據(jù)是該程序的核心問題。由于純凈液體各向同性, 可以承受極大的拉力, 其空化壓強(qiáng)與溫度基本呈線性關(guān)系, 因此有均質(zhì)空化模型

式中,phomogenous為均質(zhì)空化壓強(qiáng),Tref,0和Tref,1以及pref,0和pref,1分別為兩個參考狀態(tài)的溫度和壓強(qiáng), 本文選擇的兩個參考狀態(tài)為Tref,0=273.15 K,pref,0=-115.0 MPa 和Tref,1=288.15 K,pref,1=-100.0 MPa。

在數(shù)值計算過程中, 當(dāng)流體壓強(qiáng)p

但若液體與壁面接觸, 受壁面濕潤性影響, 壁面處液體承壓能力迅速下降, 更易達(dá)到空化條件, 因此有異質(zhì)空化模型

式中,pheterogenous為異質(zhì)空化壓強(qiáng),p′為飽和蒸氣壓,ξ為壁面與液滴之間的靜態(tài)接觸角。

當(dāng)流體狀態(tài)滿足相變條件時, 則觸發(fā)數(shù)值模型中的相變過程, 利用基于化學(xué)勢的瞬間相變模型求解相變源項

1.4 小結(jié)

本節(jié)總結(jié)了實驗室針對高瞬變強(qiáng)可壓縮兩相流與反應(yīng)流系統(tǒng)所建立的數(shù)理模型與數(shù)值計算方法, 包括適用于多間斷捕捉的高精度格式, 以及基于化學(xué)弛豫的瞬變相變(空化)模型。這些新的模型與算法為強(qiáng)可壓縮兩相流過程的數(shù)值模擬提供了新方法。

2 兩相流體界面約束波系解析方法及運(yùn)用

激波與兩相界面的相互作用在激波沖擊液滴、 激波沖擊氣泡、 高速液滴撞擊壁面等過程中普遍存在。激波和兩相流體界面的相互作用過程更加復(fù)雜, 其與激波和固體壁面的相互作用不同, 除激波反射外, 還存在透射現(xiàn)象, 且透射波的傳播與流體介質(zhì)相關(guān)。理解激波與兩相界面的相互作用機(jī)理, 厘清激波在兩相界面的反射與透射類型, 及其發(fā)生條件, 是詳細(xì)解析復(fù)雜波系與兩相界面時空演化的前提。

本節(jié)重點介紹兩相流體界面約束波系的解析方法, 以及運(yùn)用該方法分析激波沖擊液滴以及液滴撞擊壁面的過程。

2.1 兩相界面的激波反射與透射類型

項高明等[14,19]研究分析了激波在氣液兩相界面反射與透射類型, 如圖3所示, 左側(cè)為氣相區(qū)域, 右側(cè)為液相區(qū)域, 紅色斜線為氣液相界面, 黑色線為入射激波, 定義入射激波與相界面之間的夾角為入射激波角β。

圖3 激波與相界面相互作用示意圖Fig. 3 Interaction between shock waves and phase interfaces

利用解析理論分析激波與氣液兩相流體界面相互作用的反射與透射過程, 可將其分為“常規(guī)反射+非前導(dǎo)型透射”“常規(guī)反射+前導(dǎo)型透射”和“Mach反射+前導(dǎo)型透射”3種類型, 如圖4所示。在相同入射激波強(qiáng)度條件下, 隨入射激波角β增大, 3種類型依次出現(xiàn)。

以激波Mach數(shù)Ma=2.4為例, 激波與相界面相互作用仿真結(jié)果見圖5。當(dāng)入射角β=30°時發(fā)生常規(guī)反射+非前導(dǎo)型透射, 透射激波為直線, 此時透射激波沿兩相界面的速度分量與入射激波沿兩相界面的速度分量一致, 因此透射激波、 入射激波和反射激波相交于同一點, 見圖5(a); 當(dāng)入射角β=40°時發(fā)生常規(guī)反射+前導(dǎo)型透射, 此時透射激波沿兩相界面的速度分量超過入射激波沿兩相界面的速度分量, 此時透射激波為前導(dǎo)透射激波, 見圖5(b); 當(dāng)入射角β=50°～60°時發(fā)生Mach反射+前導(dǎo)型透射, 見圖5(c)和(d)。

(a) Regular reflection + non-leading transmission

(b) Regular reflection + leading transmission

(c) Mach reflection + leading transmission圖4 激波與相界面相互作用類型示意圖Fig. 4 Schematic of wave patterns for interaction between shock waves and phase interfaces

(a) β=30° (b) β=40° (c) β=50° (d) β=60°圖5 不同入射角度激波與相界面相互作用數(shù)值仿真結(jié)果[14]Fig. 5 Numerical results of interaction between shock waves and phase interfaces at different incident angles[14]

(a) Critical angle for transition from regular reflection to Mach reflection

(b) Critical angle for transition from non-leading transmission to leading transmission圖6 激波Mach數(shù)與Mach反射臨界角關(guān)系圖[14]Fig. 6 Relationship between Mach number and critical reflection angle[14]

由圖6可知: 1)氣液兩相界面情況中Mach反射臨界角高于固壁情況, 且臨界角隨激波Mach數(shù)提高而趨于定值; 2)激波Mach數(shù)越高, 透射激波由非前導(dǎo)型向前導(dǎo)型轉(zhuǎn)變的臨界激波角越大。

2.2 高速液滴/壁面相互作用過程分析

燃燒室霧化燃料液滴碰撞壁面過程影響噴霧特性以及液膜在壁面的鋪展規(guī)律。研究發(fā)現(xiàn), 在液滴高速碰撞壁面過程中液滴內(nèi)部存在著水錘激波產(chǎn)生、 激波與流體界面作用、 液滴內(nèi)空化等現(xiàn)象, 需要探究界面約束水錘激波產(chǎn)生機(jī)制、 運(yùn)動激波與流體界面作用規(guī)律和液滴內(nèi)空化初生機(jī)理等問題。Wu等[18,20]、 Xiang等[19]研究了液柱撞擊剛性水平壁面的過程, 以液柱高速撞擊剛性水平壁面為例, 見圖7, 其中點C為液柱中心, 點O為液柱與壁面的初始撞擊點, 液柱半徑R0=5 mm, 液柱初始速度為V0=110 m/s。

圖7 液柱撞擊剛性水平壁面示意圖Fig. 7 Schematic of the high-speed liquid column impinging on rigid flat surface

Wu等[18]建立了與實驗[21]相對應(yīng)的數(shù)值仿真模型, 并利用前文所介紹的數(shù)值模擬方法獲得撞擊過程的數(shù)值結(jié)果。根據(jù)流動特征以及對典型流動現(xiàn)象的分析, 撞擊過程可分為3個階段: (a)水錘激波產(chǎn)生與脫離階段, (b)復(fù)雜非定常波系演化階段, (c)反射波匯聚與空化泡產(chǎn)生階段。下面詳細(xì)介紹各個階段的特征。

2.2.1 水錘激波產(chǎn)生與脫離階段

初始階段開始于撞擊接觸時刻, 結(jié)束于水錘激波脫離壁面的時刻。

液柱撞擊剛壁面后, 接觸區(qū)端點A將沿碰撞平面以速度VA移動, 見圖8所示, 定義接觸區(qū)端點A處液柱切線與剛性壁面之間的夾角為表面夾角θ。

由圖8可知, 隨表面夾角θ增大, 接觸區(qū)端點A的橫向速度VA逐漸減小。由Rein[22]的研究結(jié)論可知, 激波離開壁面前的臨界角度計算公式為

式中,V0為液柱初始速度,VA,c為臨界條件下點A的速度值,Vs為受限激波初始波速。計算結(jié)果表明在液柱初始速度V0=110 m/s的情況下, 臨界角θc≈4°。

進(jìn)一步, 參考Heymann[23]建議的估算方法, 計算獲得在初始撞擊點O處產(chǎn)生的水錘激波壓強(qiáng)值為

ph=ρ0V0Vs～100 MPa

式中,ph為水錘激波壓強(qiáng),ρ0為液柱初始密度。

圖8 表面夾角和接觸點橫向速度示意圖Fig. 8 Schematic diagram of surface angle and lateral velocity of contact point

利用小波理論和Newton射線分析法, 可獲得水錘激波在液柱內(nèi)的演化特性, 進(jìn)而構(gòu)建界面約束波系解析理論。在液柱與壁面碰撞過程中, 任意時刻沿壁面不斷鋪展的液柱在與壁面接觸的端點A處會產(chǎn)生一個以當(dāng)?shù)芈曀傧蛲膺\(yùn)動的弧形小波, 小波內(nèi)區(qū)域為受撞擊影響區(qū)域, 小波外區(qū)域為未受擾動區(qū)域。在碰撞的最開始階段小波運(yùn)動的速度(即當(dāng)?shù)芈曀?有限而VA更大, 因此VA的值大于小波的運(yùn)動速度, 即小波的運(yùn)動速度無法超越小波的產(chǎn)生速度, 小波被限制在液柱內(nèi)部。這些小波的包絡(luò)線即構(gòu)成了水錘激波的波面St。

臨界時刻后瞬間, 水錘激波波面St和小波波系脫離, 此時的激波波面位置和壓縮小波的波系結(jié)構(gòu)見圖9所示。

圖9 臨界時刻后瞬態(tài)波系示意圖Fig. 9 Transient wave structure after the critical time

圖10 水錘激波復(fù)雜非定常波系初始演化局部分析示意圖Fig. 10 Local analysis of initial evolution of complex unsteady wave structure with water hammer shock wave

2.2.2 復(fù)雜非定常波系演化階段

受限激波脫離固壁后, 繼續(xù)向遠(yuǎn)離固體壁面方向運(yùn)動, 但由于受限激波轉(zhuǎn)變?yōu)閿U(kuò)展型激波, 強(qiáng)度隨波面?zhèn)鞑ザ鴾p弱, 波后壓強(qiáng)隨激波運(yùn)動而逐漸下降。該過程伴隨著液柱內(nèi)部形成的更為復(fù)雜的非定常波系演化過程。

由t0時刻產(chǎn)生小波后, 在tc時刻所對應(yīng)的射線簇長度為

式中,c為當(dāng)?shù)芈曀佟?jù)此, 可以獲得未發(fā)生反射射線的發(fā)射角α取值范圍為

式中,R為液柱半徑。發(fā)生N次反射后, 射線發(fā)射角α取值范圍為

當(dāng)反射激波運(yùn)動到液柱上半部分時, 圖11所示為t*=0.81(t*=t2) 時的波系結(jié)構(gòu)示意圖。圖11(a)為數(shù)值計算獲得的紋影圖, 圖11(b)為波系的解析示意圖。通過射線分析, 可以得到一次反射稀疏波波頭的形狀和位置。

(a) Numerical schlieren contour

(b) Schematic of ray analysis圖11 反射激波在液柱上半部分時波系示意圖Fig. 11 Schematic of initial evolution of reflected shock waves at upper part of the cylindrical water

圖11(b)給出了部分分別對應(yīng)于未被反射的受限激波的射線、 一次反射稀疏波的射線和二次反射壓縮波的射線?？梢钥闯? 這些射線可以根據(jù)發(fā)射角α的取值區(qū)間劃分為5個區(qū)域, 包括: (Ⅰ)0<α<α0, (Ⅱ)α0<α<α1, (Ⅲ)α1<α<α2, (Ⅳ)α2<α<α3和(Ⅴ)α3<α<α4。 其中, 區(qū)域(Ⅰ)射線端點對應(yīng)受限激波波面(St2)位置, 其發(fā)射角α上限值為α0=arccos[c(t2-tc)/(2R)]; 區(qū)域(Ⅱ)和區(qū)域(Ⅲ)對應(yīng)于發(fā)生了一次反射的射線, 區(qū)域(Ⅱ)射線端點對應(yīng)一次反射稀疏波上支(見圖12藍(lán)色射線)的波頭位置, 發(fā)射角α的上限角α1介于α0和α2之間, 區(qū)域(Ⅲ)射線端點對應(yīng)一次反射稀疏波的下支(如圖12藍(lán)色射線)的波頭位置, 發(fā)射角α上限值為α2=arccos[c(t2-tc)/(4R)]; 區(qū)域(Ⅳ)和區(qū)域(Ⅴ)對應(yīng)于發(fā)生了二次反射的射線, 區(qū)域(Ⅳ)射線端點對應(yīng)二次反射壓縮波上支(如圖12紅色射線)的波頭位置, 發(fā)射角α的上限值α3介于α2和α4之間, 區(qū)域(Ⅴ)射線端點對應(yīng)二次反射壓縮波下支(如圖12紅色射線)的波頭位置, 發(fā)射角α的上限值為α4=arccos[c(t2-tc)/(6R)]。 圖12為數(shù)值密度紋影圖和射線理論分析結(jié)果的對比圖(中)和壓強(qiáng)等值線云圖的數(shù)值結(jié)果(右)。通過射線分析得到的一次反射稀疏波的波頭位置與數(shù)值紋影結(jié)果吻合良好。二次反射壓縮波波頭相對較弱, 在數(shù)值紋影結(jié)果中較難分辨。在壓強(qiáng)等值線結(jié)果中, 所顯示的局部壓強(qiáng)極大值區(qū)域即對應(yīng)了二次反射壓縮波波頭的位置, 其形狀和位置也與射線分析結(jié)果相吻合。

圖12 反射波系解析圖Fig. 12 Analytical diagram of reflected wave system

2.2.3 反射波匯聚與空化泡產(chǎn)生階段波系

通過上述解析方法求解反射波匯聚過程和波系誘導(dǎo)空化區(qū)域, 見圖13; 圖13左側(cè)為解析方法獲得的反射波波系圖, 右側(cè)為數(shù)值紋影圖。

圖13 空化區(qū)域分析圖Fig. 13 Schematic of cavitation zone

通過解析理論獲得極限交點Pdli與液柱上極點TP間的極限距離εmax為

式中,D0為液柱的初始直徑。

對比1989年Cavendish實驗室Field等[21]實驗觀察到的空化現(xiàn)象, 上述分析過程系統(tǒng)給出了空化初生位置理論解, 與實驗觀測結(jié)果相符合。

2.2.4 壁面曲率效應(yīng)

考慮壁面曲率影響, Wu等[24]進(jìn)一步研究了高速液柱撞擊曲面壁面的動態(tài)過程, 采用含流體相變的三組分可壓縮多相流動模型, 獲得了液柱高速撞擊(V0=150 m/s)過程中的復(fù)雜波系結(jié)構(gòu)和空化過程, 見圖14, 圖14(a)為高速液柱撞擊凹壁面情況, 圖14(b)為高速液滴撞擊凸壁面情況。

(a) Concave surface

(b) Convex surface圖14 液柱撞擊剛性曲壁面示意圖Fig. 14 Schematic of cylindrical liquid column impacting on rigid curved surface

對比凹壁面、 平面和凸壁面3種壁面條件下, 液滴撞擊過程中流場中最大壓強(qiáng)變化曲線, 見圖15。

圖15 3種壁面條件液柱撞擊過程壓強(qiáng)曲線Fig. 15 Pressure curves of cylindrical liquid columnimpacting on three types of surface

由圖15仿真獲得的壓強(qiáng)曲線結(jié)果可知, 在3種壁面條件下, 初始撞擊瞬間的水錘壓強(qiáng)pinitial_impinging基本一致, 約為250 MPa。對于液柱撞擊凹壁面而言, 液滴撞擊瞬間, 由于初始強(qiáng)水錘激波的作用, 壓強(qiáng)升高, 之后, 由于均質(zhì)聚焦空腔而產(chǎn)生壓強(qiáng)脈沖pfocus_cavity, 由于近壁面空腔而產(chǎn)生壓強(qiáng)脈沖pfirst_surface_cavity和psecond_surface_cavity等。但在液柱撞擊凸壁面時, 液柱撞擊后壓強(qiáng)顯著降低, 沒有較強(qiáng)的二次壓強(qiáng)脈沖產(chǎn)生。

針對液柱撞擊過程中壁面曲率效應(yīng)更詳細(xì)的分析請參見文獻(xiàn)[24]。

2.3 激波/液滴相互作用過程分析

激波和液滴相互作用過程中同樣存在著液滴內(nèi)部復(fù)雜波系演化過程, 同時還存在著激波在液滴兩相界面的演化過程和液滴中的射流形成過程。

2.3.1 波系結(jié)構(gòu)及界面演化過程分析

Xiang等[25]通過數(shù)值方法研究了激波撞擊液柱過程中的波系結(jié)構(gòu)與波系演化過程。以激波(Mas=2.4)沖擊初始直徑為d0的液滴為例, 定義無量綱時間t*為

式中,t為真實時間,u′2為入射激波波后初始速度。

將數(shù)值計算結(jié)果與Sembian等的實驗結(jié)果[6]進(jìn)行對比, 數(shù)值結(jié)果與實驗中的流場顯示結(jié)果吻合, 見圖16, 圖中左列為本課題組數(shù)值計算結(jié)果, 上半部分為數(shù)值紋影圖, 下半部分為壓強(qiáng)分布圖, 圖右列為Sembian等的實驗結(jié)果.

激波撞擊液滴過程的波系結(jié)構(gòu)及界面演化過程見圖17。圖中ISW代表入射激波, RSW代表反射激波, TW代表透射激波, REW代表反射膨脹激波, R-TW代表二次透射激波, SS代表滑移線, MS代表Mach桿。

(a) t*=0.17

(b) t*=0.31

(c) t*=0.62

(d) t*=0.75

(e) t*=1.20

(f)t*=1.62

圖17 平面激波與液柱相互作用過程波系結(jié)構(gòu)示意圖Fig. 17 Schematic diagram of the wave structures of a planar shock wave interacting with a water column

圖17結(jié)果表明激波與液滴作用前期波系結(jié)構(gòu)演化過程與實驗基本一致; 液滴內(nèi)部反射膨脹波在液滴下游壁面附近聚焦, 形成低壓區(qū)。入射激波的加速效應(yīng)和脫落渦的阻礙效應(yīng)共同作用使得液滴下游表面逐漸被壓平; 在液滴上游表面, 液體逐漸發(fā)生卷起形成一些渦結(jié)構(gòu); 在高速氣流的夾帶作用下, 液滴表面卷起的渦結(jié)構(gòu)逐漸被拉伸, 最終發(fā)生剝落流向下游。這種作用在液滴上游頂點最為明顯, 大量液體先卷起, 后被高速氣流夾帶從液滴上剝落, 最終使得液滴發(fā)生破碎。因此, 在本文算例中液滴后續(xù)出現(xiàn)的破碎, 其機(jī)制主要為剪切誘導(dǎo)夾帶機(jī)制。

2.3.2 入射激波穿過液滴后演化規(guī)律

進(jìn)一步對液滴入射激波的調(diào)制規(guī)律進(jìn)行了研究, 主要分析入射激波穿過液滴后的形態(tài)和強(qiáng)度演化情況, 見圖18, 仿真結(jié)果表明, 入射激波穿過液滴后, 形成上下對稱的4個三波點(TP1, TP2, TP3和TP4), 其中TP1與TP3之間、 TP2與TP4之間是Mach桿, TP3和TP4之間為Mach桿相交之后形成的新激波, 稱之為恢復(fù)激波。

在不同Mach數(shù)下, 圖18中三波點TP1和TP3在不同Mach數(shù)下的運(yùn)動軌跡見圖19, TP1的運(yùn)動軌跡幾乎為一條直線, 并且受入射激波強(qiáng)度影響較小, 而TP3的運(yùn)動軌跡在運(yùn)動后期也接近直線。

圖18 三波點運(yùn)動軌跡(Ma=2.4)Fig. 18 Trajectory of triple points at Ma=2.4

不同時刻入射激波波面的形狀變化情況見圖20, 隨三波點運(yùn)動, 恢復(fù)激波寬度逐漸變大, 恢復(fù)激波波面與入射激波波面間距不斷縮小, 恢復(fù)激波逐漸追趕達(dá)到入射激波。

圖19 不同Mach數(shù)條件下三波點運(yùn)動軌跡Fig. 19 Trajectory of triple points at different Mach numbers

圖20 入射激波波面形狀變化情況Fig. 20 Shape of incident shock wave

分析上述現(xiàn)象, 繪制激波穿過液滴后的壓強(qiáng)云圖, 見圖21, 可以看出恢復(fù)激波波后壓強(qiáng)高于入射激波壓強(qiáng), 因此恢復(fù)激波具有“追趕”效應(yīng)。

圖21 激波穿過液滴后的壓強(qiáng)云圖Fig. 21 Pressure filed after shock wave passing through droplets

進(jìn)一步, 對比不同入射激波強(qiáng)度(Ma=1.47, 2.0, 2.4)條件下恢復(fù)激波中點位置波后壓強(qiáng)隨時間的變化(見圖22), 發(fā)現(xiàn)隨時間推移, 恢復(fù)激波的壓強(qiáng)逐漸下降, 最終接近入射激波波后壓強(qiáng)。

圖22 恢復(fù)激波中點位置波后壓強(qiáng)變化曲線Fig. 22 Pressure curves at the midpoint of recovery shock wave

對于3種不同入射強(qiáng)度的激波, 對比恢復(fù)激波壓強(qiáng)恢復(fù)到1.04, 1.05和1.06倍入射激波波后壓強(qiáng)p2所需的壓強(qiáng)恢復(fù)時間, 見圖23。發(fā)現(xiàn)入射激波強(qiáng)度越強(qiáng), 恢復(fù)激波恢復(fù)到入射激波波后壓強(qiáng)所需時間越長。

2.3.3 液滴運(yùn)動特性分析

研究激波作用下液滴的運(yùn)動特性, 獲得液滴質(zhì)心位移(見圖24)、 液滴質(zhì)心加速度(見圖25)、 液滴中心線流向?qū)挾?見圖26)和液滴展向高度(見圖27)的變化曲線。

圖23 恢復(fù)激波壓強(qiáng)恢復(fù)時間Fig. 23 Pressure recovery time of recovery shock wave

圖24 液滴質(zhì)心位移Fig. 24 Drift of droplet centroid

圖25 液滴無量綱質(zhì)心加速度Fig. 25 Non-dimensional acceleration of droplet centroid

圖26 液滴中心線流向?qū)挾菷ig. 26 Flow width of droplet centerline

圖27 液滴展向高度Fig. 27 Spanwise height of droplets

2.3.4 運(yùn)動波系誘導(dǎo)液滴空化理論

運(yùn)動波系與液滴相互作用過程中, 由于反射膨脹波(REW)在液滴內(nèi)部匯聚形成負(fù)壓區(qū), 一旦流體內(nèi)部局部壓強(qiáng)小于當(dāng)?shù)仫柡驼羝麎? 就會誘導(dǎo)發(fā)生空化, 見圖28。

(a) t*=0.27

(b)t*=0.36

(c)t*=0.65

(d)t*=1.03

(e)t*=1.18

(f)t*=1.22

圖28 運(yùn)動激波誘導(dǎo)流體空化過程仿真結(jié)果Fig. 28 Dynamic shock wave induced liquid cavitation process

利用小波分析法獲得了反射膨脹波在液滴內(nèi)部形成的負(fù)壓區(qū)軌跡, 見圖29所示, 分析表明反射膨脹波最終在Fe點匯聚, 導(dǎo)致此處壓強(qiáng)降低并發(fā)生空化現(xiàn)象。

(a) t*=0.34

(b) t*=0.21

(c) t*=1.03

(d) Numerical schlieren圖29 反射膨脹波在液滴內(nèi)部形成負(fù)壓區(qū)軌跡圖Fig. 29 Trajectory of reflection expansion waves forming negative pressure in droplet

進(jìn)一步分析表明匯聚點Fe位置與臨界入射角βc具有如下關(guān)系

2.4 激波/含空泡液柱射流產(chǎn)生機(jī)制

Xiang等[25]針對激波與含空泡液滴相互作用, 建立了數(shù)值仿真模型, 見圖30。圖中深藍(lán)色區(qū)域為液柱, 中心為空泡(空氣泡或蒸汽泡), 入射激波向右方傳播。

研究激波與含空泡液滴相互作用過程中的射流形態(tài)變化規(guī)律。在Mas=2.4激波作用下, 含空泡液柱界面形態(tài)演化過程見圖31所示。

圖30 激波與含空泡液柱相互作用簡化模型Fig. 30 Simplified model of interaction between shock wave and liquid column with cavity

圖31 激波作用下含空泡液柱界面形態(tài)演化過程Fig. 31 Interface evolution of the water ring interacting with planar shock wave

繪制不同半徑比下, 不同強(qiáng)度入射激波(分別為Mas=1.47,2.0,2.4)作用下, 射流形成瞬間的兩相界面形狀, 見圖32。

圖32 不同半徑比和激波強(qiáng)度作用下射流形成瞬間界面形狀Fig. 32 Interface shape at the moment of jet formation at several radius ratios and shock wave intensities

圖31計算結(jié)果表明: 1)入射激波越強(qiáng), 液滴內(nèi)部氣泡壓縮效應(yīng)越明顯, 射流寬度越寬; 2)半徑比為0.25時, 入射激波強(qiáng)度越弱, 射流刺穿氣泡所需要的時間越長; 3)半徑比為0.5和0.75時, 不同入射激波強(qiáng)度下射流刺穿氣泡時液環(huán)的寬度變化不大。

進(jìn)一步, 還研究了不同類型空泡(空氣泡和蒸氣泡)對射流產(chǎn)生過程的影響規(guī)律[20]?？张萆淞餍螤顚Ρ纫妶D33, 撞擊后空泡被分割為一個在中心軸上的小空泡和兩個較大的軸對稱分布在中心軸兩側(cè)的空泡翼, 在空泡潰滅后, 剩余的大空泡翼被進(jìn)一步分割成一系列獨(dú)立的小空氣泡。

圖33 不同類型氣泡中射流形成形狀Fig. 33 Shapes of jet formation in different types of bubbles

對比不同類型氣泡潰滅后剩余氣泡的演化過程[20], 見圖34。

圖34仿真結(jié)果表明, 潰滅后, 空氣泡被分割, 兩側(cè)軸對稱分布的空氣泡翼為尖頭圓尾形狀, 中心軸上的殘余空氣泡則非常小。對于兩側(cè)較大的空氣泡翼而言, 在主縱向射流和界面碰撞產(chǎn)生的水錘壓強(qiáng)波的作用下, 空氣泡翼尖端被持續(xù)壓縮, 空氣泡內(nèi)的氣體在被壓縮的過程中高速排向空泡鈍頭側(cè)(光滑的尾端), 并形成排氣射流。排氣射流的最高速度達(dá)到1 500 m/s。與此同時, 空氣泡翼由被壓縮的尖頭向鈍頭尾端橫向潰縮, 使得液體橫向流動產(chǎn)生橫向射流, 液體橫向射流速度可達(dá)600 m/s。

圖34 潰滅后不同類型氣泡液滴剩余氣泡的演化過程Fig. 34 Evolution of remaining cavities of the droplet with different cavities after collapsing

2.5 小結(jié)

本節(jié)介紹了兩相流體界面約束波系解析的理論方法, 通過理論解析和數(shù)值方法研究了高速液滴/壁面、 激波/液滴、 激波/含空泡液柱相互作用過程中的多波系演化過程, 解析獲得了波系誘導(dǎo)空化初生位置的理論解, 實現(xiàn)了界面約束波系和波系誘導(dǎo)空化過程的高精度數(shù)值解析, 其結(jié)果與文獻(xiàn)中實驗結(jié)果一致。

3 可壓縮反應(yīng)流動多物理過程

3.1 強(qiáng)可壓氣液兩相含化學(xué)反應(yīng)流動的數(shù)值方法

超聲速噴霧燃燒過程中存在湍流與激波作用下液霧的彌散與蒸發(fā)、 兩相混合、 著火與火焰?zhèn)鞑サ戎T多復(fù)雜物理化學(xué)過程。在噴霧燃燒中, 霧化的燃料液滴直徑通常為微米量級, 一般遠(yuǎn)小于網(wǎng)格尺寸, 使用Euler-Lagrange方法可以較好地捕獲流體和液滴的運(yùn)動過程。

因此, 清華大學(xué)噴霧燃燒與推進(jìn)實驗室整合近十年的噴霧燃燒仿真經(jīng)驗積累[26-35], 發(fā)展了一套高精度數(shù)值模擬軟件TURFsim[36-37], 該軟件可實現(xiàn)全速域高精度大渦模擬噴霧燃燒求解。該軟件采用Cartesian網(wǎng)格, 使用有限差分法進(jìn)行方程離散, 應(yīng)用浸沒邊界法處理復(fù)雜幾何構(gòu)型, 采用亞格子模型封閉大渦模擬湍流項, 使用Euler-Lagrange方法求解噴霧兩相流動, 采用簡化或者詳細(xì)化學(xué)反應(yīng)機(jī)理描述燃料的化學(xué)反應(yīng)過程。

在數(shù)值格式方面, TURFsim軟件集成了1～6階精度的數(shù)值格式, 同時集成了一部分高精度混合數(shù)值格式。本文主要介紹針對強(qiáng)可壓氣液兩相流動采用的數(shù)值格式, 強(qiáng)可壓縮連續(xù)相流動中包含多尺度旋渦, 存在著激波-旋渦的相互作用, 要求數(shù)值格式能夠捕捉流動中的復(fù)雜激波波系, 同時需要分辨湍流流動。任兆欣[27]采用Hu等[38]提出的6階精度自適應(yīng)中心-迎風(fēng)WENO混合格式離散氣相控制方程中的對流項, 能夠使用WENO格式捕捉激波, 而在光滑區(qū)(湍流)使用中心-迎風(fēng)格式, 采用6階精度對稱型緊致差分格式離散氣相控制方程中的黏性擴(kuò)散項, 采用3階顯式Runge-Kutta積分方法進(jìn)行時間推進(jìn)。

在TURFsim軟件的Euler-Lagrange方法中, 液滴離散相與氣體連續(xù)相之間的質(zhì)量、 動量和能量交換采用雙向耦合模型。針對液滴彌散相, 采用Lagrange控制方程描述, 將液滴近似為點源模型(遵從點源模型假設(shè)), 在Lagrange框架下求解描述液滴物理量的常微分方程

(1)

(2)

(3)

式中,xd,i和ud,i分別為液滴在i方向的位置和速度,Tp和md分別為液滴溫度和質(zhì)量,ui是i方向氣相速度,Fd,i是液滴在i方向的受力,cp是氣相混合物的定壓比熱容,cL是燃料液體的比熱,LV是液體在標(biāo)準(zhǔn)沸點溫度條件下的蒸發(fā)潛熱,BM是傳質(zhì)系數(shù)。

基于氣液兩相滑移速度定義當(dāng)?shù)匾旱蜶eynolds數(shù)Red

式中,ρ是氣相混合物的密度,μ是氣相混合物的黏性系數(shù),dd是液滴直徑。

式(1)中f1是對Stokes阻力定律的修正, 經(jīng)驗公式為

液滴運(yùn)動的特征弛豫時間τd定義為

其中,ρd是液滴的密度。

式(2)中, 描述對流效應(yīng)對傳熱過程影響的Nusselt數(shù)Nu為

式中,Pr為氣相Prandtl數(shù), 定義為

式中,λ是氣相混合物的熱傳導(dǎo)系數(shù)。

式(2)中，f2是液滴蒸發(fā)對傳熱過程的修正系數(shù), 定義為

式中, 無量綱蒸發(fā)參數(shù)β定義為

式(3)中描述對流效應(yīng)對傳質(zhì)過程影響的Sherwood數(shù)Sh為

式中,Sc為氣相Schmidt數(shù), 定義為

通過上述Lagrange形式的液滴控制方程, 獲得液滴的質(zhì)量、 速度、 溫度、 質(zhì)量變化和運(yùn)動軌跡。

本文采用稀疏相假設(shè), 使用雙向耦合模型, 即忽略液滴之間的相互作用, 描述液滴蒸發(fā)過程與周圍氣相的質(zhì)量、 動量和能量交換, 由氣相網(wǎng)格節(jié)點數(shù)據(jù)線性插值獲得液滴位置處的流場物理信息(如氣體流動速度、 溫度、 壓強(qiáng)等), 通過插值方法將液滴對周圍氣相的質(zhì)量、 動量和能量的傳遞賦予液滴周圍的氣相網(wǎng)格節(jié)點。點源模型采用了可壓縮修正, 詳細(xì)方法可參見文獻(xiàn)[27]。

3.2 斜激波增強(qiáng)混合機(jī)理分析

超燃沖壓發(fā)動機(jī)燃燒室內(nèi)氣體流動時間僅為毫秒量級, 因此增強(qiáng)燃料(噴霧燃油)和氧化劑(空氣)的混合, 有利于提高發(fā)動機(jī)整體效率。利用超燃沖壓發(fā)動機(jī)內(nèi)流道中的斜激波等波系促進(jìn)燃料和氧化劑混合成為一種經(jīng)濟(jì)可行的技術(shù)途徑。

薛淑艷[39]將該問題抽象為斜激波對超聲速平面混合層流動的混合增強(qiáng)問題, 物理模型見圖35。

圖35 斜激波作用下超聲速混合層示意圖Fig. 35 Mixing layer with interaction with oblique shock wave

超聲速平面混合層的顯著特征是剪切導(dǎo)致的大渦擬序結(jié)構(gòu), 大渦擬序結(jié)構(gòu)控制著混合層中的質(zhì)量、 動量和能量的傳遞過程。當(dāng)斜激波與超聲速混合層相互作用后, 混合層中大渦擬序結(jié)構(gòu)產(chǎn)生響應(yīng), 同時, 激波入射混合層后發(fā)生折射和反射, 形成復(fù)雜波系結(jié)構(gòu)。研究斜激波對超聲速混合層增混機(jī)理的關(guān)鍵在于揭示激波與旋渦之間的相互作用機(jī)制。

3.2.1 旋渦穿過激波時激波的演化過程

在旋渦的作用下, 激波會發(fā)生Mach反射或正規(guī)反射。以混合層中旋渦穿過激波過程中激波演化過程為例, 激波變形數(shù)值紋影圖見圖36。該數(shù)值紋影圖中使用變量q1對流場中的激波和弱波進(jìn)行表征， 其計算公式為

式中,qmax為全場中q的最大值,q1為無量綱數(shù), 取值范圍為[0.36,0.8], 激波強(qiáng)度越強(qiáng), 密度梯度越大, 則q1越小; 反之, 流場越光滑, 則q1越大。

追蹤圖36(a) 中旋渦A與激波的相互作用, 圖中旋渦A沿流向方向(x方向)運(yùn)動。當(dāng)渦邊緣與激波接觸時, 二者相互作用很弱, 因此激波基本沒有變形, 見圖36(a)。當(dāng)渦心運(yùn)動到激波附近時, 渦心處形成折射激波, 渦兩側(cè)與激波運(yùn)動方向相反的位置出現(xiàn)了慢速衍射激波(SD), 此處的激波強(qiáng)度增加, 與激波運(yùn)動方向相同的位置出現(xiàn)了快速衍射激波(FD), 此處的激波強(qiáng)度削弱。激波在渦的作用下呈現(xiàn)出S形結(jié)構(gòu), 見圖36(b)。隨著渦繼續(xù)運(yùn)動, 穿過激波, S形結(jié)構(gòu)更加顯著, 見圖36(c)和(d)。當(dāng)渦心遠(yuǎn)離激波時, 在激波強(qiáng)度增強(qiáng)的慢速衍射激波(SD)處出現(xiàn)了Mach反射, 形成Mach桿結(jié)構(gòu)(MS), 見圖36(e)。在Mach桿的壓強(qiáng)梯度作用下反射形成二次激波MR1和MR2, 見圖36(f)。

(a) Initial interaction

(b) Vortex center at shock wave

(c) Vortex center passing through shock wave

(d) Vortex center passed through shock wave

(e) Mach stem formation

(f) Secondary shock wave formation

3.2.2 旋渦穿過激波時的演化過程

在激波作用下, 旋渦受到壓縮, 渦量增強(qiáng), 見圖37。在未穿過激波時, 混合層中旋渦渦心處的渦量處于-39 000～-31 000 s-1之間, 而有激波作用后, 混合層旋渦渦心處渦量顯著增強(qiáng), 處于-54 000～ -40 000 s-1之間。

(a) Case without shock wave

(b) Case with shock wave圖37 有無激波情況時的渦量分布圖Fig. 37 Vorticity distribution with/without shock wave

薛淑艷[39]通過分析二維、 無黏、 無體積力的渦量動力學(xué)方程, 研究了激波增混的斜壓效應(yīng)。首先通過對動量方程求旋度, 獲得渦量動力學(xué)方程

3.2.3 激波作用下的混合特性

定義混合層厚度δ, 用其表征混合程度。在激波作用下, 受到壓縮作用, 混合層厚度先減小, 但混合層厚度增長率增大, 最終, 混合層厚度超過無激波作用時的厚度, 見圖38。

圖38 有無激波條件下混合層厚度變化曲線Fig. 38 Mixing layer thickness curve with/without shock wave

沿流向選取4個不同位置截面, 繪制混合層的湍動能, 見圖39, 圖39(a)為x=0.6 m, 為激波前位置; 圖39(b)為x=0.77 m, 為激波與混合層作用區(qū); 圖39(c)為x=0.8 m, 為剛離開激波的混合層區(qū); 圖39(d)為x=1.0 m, 為遠(yuǎn)離激波的混合層區(qū).

(a) x=0.6 m

(b) x=0.77 m

(c) x=0.8 m

(d) x=1.0 m

由圖39可知, 激波與混合層作用過程中, 混合層湍動能峰值迅速升高, 而且受激波壓縮作用, 湍動能分布范圍變窄; 當(dāng)混合層中心離開激波時, 湍動能呈雙峰結(jié)構(gòu), 左側(cè)尖銳峰為激波所處位置, 右側(cè)光滑峰為混合層中心所處位置; 當(dāng)混合層遠(yuǎn)離激波后, 混合層的湍動能較作用前增強(qiáng)。

同時, 繪制對應(yīng)位置的橫向脈動強(qiáng)度分布圖, 如圖40所示。

(a) x=0.6 m

(b) x=0.77 m

(c) x=0.8 m

(d) x=1.0 m

由圖40可知, 經(jīng)過激波作用, 混合層的橫向脈動強(qiáng)度也得到了一定的增強(qiáng)。

結(jié)合以上分析, 經(jīng)過激波作用, 混合層首先被壓縮, 但混合層厚度增長率變大, 因此最終混合層的厚度增大, 混合層的湍動能和流動脈動強(qiáng)度增強(qiáng), 混合層整體混合程度得到增強(qiáng)。

3.3 斜激波增混與強(qiáng)化燃燒

上述斜激波增強(qiáng)混合機(jī)理基于氣相無反應(yīng)過程進(jìn)行分析, 相關(guān)結(jié)論可以遷移到含燃料液滴噴霧燃燒過程, 但液霧兩相燃燒的燃燒效率取決于湍流流動、 液霧彌散、 蒸發(fā)與燃燒之間的關(guān)聯(lián), 因此需要考察斜激波增混和強(qiáng)化燃燒過程所涉及的旋渦對液霧彌散和蒸發(fā)的影響、 斜激波和化學(xué)反應(yīng)過程的相互影響。

3.3.1 強(qiáng)可壓縮剪切流動中液霧彌散與蒸發(fā)

在強(qiáng)可壓縮流動中存在復(fù)雜的膨脹波與壓縮波, 甚至小激波, 導(dǎo)致旋渦流場中交替出現(xiàn)局部的高溫高壓區(qū)與低溫低壓區(qū), 典型的密度梯度模量和無量綱溫度分布圖見圖41[27]。圖41(a)中綠色虛線為無量綱展向渦量等值線, 明亮條帶狀結(jié)構(gòu)即為小激波, 可以看出小激波存在于旋渦之間的強(qiáng)拉伸區(qū)域; 圖41(b)表明高速強(qiáng)可壓效應(yīng)引起高溫區(qū)(流動壓縮區(qū))和低溫區(qū)(流動膨脹區(qū))的空間交替分布結(jié)構(gòu)。

(a) Density gradient modulus

(b) Non-dimensional temperature圖41 密度梯度模量和無量綱溫度分布圖Fig. 41 Distribution of density gradient modulus and non-dimensional temperature

自高速氣流入口隨機(jī)釋放液霧, 獲得液霧的瞬時分布圖見圖42, 在大尺度旋渦卷吸作用下, 液滴發(fā)生彌散, 由于展向旋渦順時針運(yùn)動, 液滴也在旋渦作用下順時針旋轉(zhuǎn)。由圖42(a)可知, 液滴更多聚集于旋渦之間的強(qiáng)拉伸區(qū)域以及旋渦外緣, 而且由于小激波的存在, 強(qiáng)拉伸區(qū)域中產(chǎn)生條帶狀的液滴聚集區(qū)域, 形狀與小激波類似; 由圖42(b)可知, 由于旋渦間局部高溫區(qū)的存在, 液滴進(jìn)入到高溫區(qū)后被當(dāng)?shù)貧饬骷訜? 溫度升高, 同時旋渦間的高壓環(huán)境提高了液滴沸點, 液滴的濕球溫度隨之升高, 導(dǎo)致液滴能達(dá)到更高溫度。此外, 旋渦內(nèi)的局部低溫氣流會對在其中運(yùn)動的液滴降溫。

(a) Dispersed liquid mist

(b) Non-dimensional temperature of droplets圖42 彌散液霧分布和液滴無量綱溫度分布圖Fig. 42 Distribution of dispersed liquid mist and non-dimensional temperature of droplets

觀察燃料液滴在流動過程中的蒸發(fā)過程, 得到無量綱燃料蒸氣質(zhì)量分?jǐn)?shù)分布圖, 見圖43。分析發(fā)現(xiàn)燃料蒸氣質(zhì)量分?jǐn)?shù)分布直接受液滴的彌散分布影響, 燃料蒸氣富集于旋渦間的強(qiáng)拉伸區(qū)域以及旋渦外緣, 在旋渦外呈條帶狀富集。因此, 在高速強(qiáng)可壓縮剪切流動中, 液滴的傾向性彌散更加顯著, 直接影響當(dāng)?shù)厝剂险魵夥植? 進(jìn)而控制著燃料-氧化劑的局部摻混以及著火與燃燒特性。

圖43 無量綱燃料蒸氣質(zhì)量分?jǐn)?shù)瞬時分布圖Fig. 43 Instantaneous distribution of non-dimensional fuel vapor mass fraction

3.3.2 斜激波與含可燃預(yù)混氣旋渦相互作用

對斜激波增混與強(qiáng)化燃燒進(jìn)行數(shù)值模擬[40-41]。對比研究了無斜激波、 斜激波角β=27°、 斜激波角β=37°的3種工況。3種工況超聲速空氣氣流在入口的溫度T0=1 500 K, 壓強(qiáng)P0=0.1 MPa, 噴霧當(dāng)量比Φ0=1.2。 同時, 定義無量綱氣體溫度

T*=T/T0

式中,T代表氣體溫度,T0代表氣流入口初始溫度。

計算獲得的無量綱氣體溫度分布和燃料蒸氣質(zhì)量分?jǐn)?shù)分布結(jié)果見圖44, 圖中色標(biāo)表示無量綱溫度值, 灰色圓點表示燃料顆粒, 藍(lán)色虛線表示燃料質(zhì)量分?jǐn)?shù)YF=0.05, 藍(lán)色虛線內(nèi)的區(qū)域表示燃料質(zhì)量分?jǐn)?shù)YF>0.05。3 種工況下, 整個流場中無量綱氣體溫度T*在關(guān)于混合分?jǐn)?shù)Z中的散點分布見圖45, 圖中藍(lán)色虛線表征絕熱火焰溫度(對有激波存在的工況, 基于波后氣體溫度計算絕熱火焰溫度)。在無激波作用時, 燃燒釋熱和蒸發(fā)冷卻的競爭機(jī)制, 富燃混合氣中液霧蒸發(fā)冷卻作用明顯(低于環(huán)境溫度T0), 無量綱氣體溫度峰值出現(xiàn)于貧燃混合氣, 見圖45(a); 在斜激波角β=27°工況中, 斜激波作用導(dǎo)致化學(xué)反應(yīng)加強(qiáng), 增強(qiáng)了釋熱速率, 因此無量綱氣體溫度峰值右移到接近化學(xué)計量條件, 見圖45(b); 斜激波角度β=37°時, 激波增強(qiáng)燃燒效應(yīng)更加明顯, 激波前后的溫度散點分布發(fā)生分離, 見圖45(c), 由于斜激波增強(qiáng)作用, 波后流場中存在部分流體微團(tuán)高于波后氣體平均溫度, 所以其火焰溫度高于基于波后氣體平均溫度計算的絕熱火焰溫度。

(a) No oblique shock wave

(b) Oblique shock wave angle β=27°

(c) Oblique shock wave angle β=37°圖44 無量綱氣體溫度T*與燃料蒸氣質(zhì)量分?jǐn)?shù)YF分布圖Fig. 44 Distribution of non-dimensional gas temperature T* and fuel vapor mass fraction YF

(a) No oblique shock wave

(b) β=27°

(c) β=37°圖45 無量綱溫度T*關(guān)于混合分?jǐn)?shù)Z的散點分布圖Fig. 45 Scatter distribution of non-dimensional temperature T* over mixing fraction Z

在斜激波角β=27°工況中, 無量綱溫度T*的瞬時分布圖見圖46, 給出了斜激波誘導(dǎo)下旋渦中可燃混合氣點火的時間演化過程, 圖中藍(lán)色虛線表征燃料蒸氣質(zhì)量分?jǐn)?shù)。圖46表明, 剪切旋渦穿過入射斜激波后, 旋渦中的可燃混合氣沒有立即點火, 只有在激波作用區(qū)域下游分布于渦辮中的混合氣被消耗并形成渦辮火焰。

進(jìn)一步觀察斜激波角β=27°工況下波系結(jié)構(gòu)演化過程, 其數(shù)值紋影圖見圖47。由圖47可知, 在強(qiáng)旋渦和弱激波的相互作用下, 原斜激波經(jīng)歷了明顯的變形過程, 當(dāng)剪切旋渦運(yùn)動到激波面, 入射激波面被扭曲, 旋渦上半部分處產(chǎn)生衍射激波, 見圖47(a); 衍射激波與旋渦相互作用形成了穿越旋渦的折射激波, 衍射激波和折射激波的交點處壓強(qiáng)較高, 具有強(qiáng)烈的壓縮作用, 因此在下游產(chǎn)生了化學(xué)反應(yīng)核心, 見圖47(b); 然而弱折射激波不能穿過旋渦, 激波-旋渦作用產(chǎn)生的波系結(jié)構(gòu)與正規(guī)反射激波結(jié)構(gòu)相似, 旋渦中心可燃混合氣不發(fā)生點火, 見圖47(c); 但隨著旋渦穿過斜激波, 衍射激波在流向方向不斷拉伸, 在兩個相鄰旋渦之間誘導(dǎo)出局部著火核心, 見圖47(d); 在該強(qiáng)可壓縮剪切流動中形成的小激波也會誘導(dǎo)旋渦周圍產(chǎn)生化學(xué)反應(yīng)核心, 見圖47(e)和(f).通過對斜激波角β=27°工況火焰核心形成過程分析, 發(fā)現(xiàn)激波和旋渦相互作用時, 激波壓縮作用導(dǎo)致的升溫過程, 誘導(dǎo)產(chǎn)生了火焰核心, 可將該過程歸為局部“熱自燃”狀態(tài), 該過程發(fā)生在激波作用區(qū)的下游。

入射激波強(qiáng)度增大, 在斜激波角β=37°工況下無量綱溫度T*的瞬時分布圖見圖48。由圖48可知, 隨斜激波角變大, 激波強(qiáng)度提升, 激波和旋渦的相互作用更加劇烈, 導(dǎo)致波后氣體溫度和壓強(qiáng)提升, 旋渦下半部分穿越激波時, 火焰立即出現(xiàn), 此時旋渦上半部分外緣也已出現(xiàn)了邊緣火焰。

(a) τ=0

(b) τ=150Δt

(c) τ=200Δt

(d) τ=250Δt

(e) τ=300Δt

(f) τ=350Δt

(a) τ=0

(b) τ=150Δt

(c) τ=200Δt

(d) τ=250Δt

(e) τ=300Δt

(f) τ=350Δt

(a) τ=0

(b) τ=150Δt

(c) τ=200Δt

(d) τ=250Δt

(e) τ=300Δt

(f) τ=350Δt

該工況對應(yīng)的數(shù)值紋影圖見圖49。在該工況下, 斜激波和強(qiáng)可壓縮旋渦相互作用時, 旋渦上半部分周圍形成彎曲的衍射激波, 與其相交的折射激波穿過旋渦核心導(dǎo)致強(qiáng)度變?nèi)? 而旋渦下半部分產(chǎn)生了一道強(qiáng)衍射激波, 與另一道反射激波相交, 在二者交點下游生成了著火核心, 見圖49(b); 旋渦下部產(chǎn)生的強(qiáng)衍射激波向上游傳播, 與發(fā)射激波分離, 傳播過程中壓縮并點燃旋渦下部的可燃?xì)? 激波面和化學(xué)反應(yīng)面相耦合, 形成局部的類似“爆震燃燒”現(xiàn)象, 見圖48(c)和49(c); 之后, 隨著燃料氣體的耗盡, 波后的化學(xué)反應(yīng)核心逐漸消失, 旋渦上部和下部的兩道衍射激波相交, 在交匯點處產(chǎn)生壓強(qiáng)峰值, 見圖49(d), (e)和(f)。

(a) τ=0

(b) τ=150Δt

(c) τ=200Δt

(d) τ=250Δt

(e) τ=300Δt

(f) τ=350Δt

3.3.3 局部近似爆震燃燒發(fā)生機(jī)理及分析

進(jìn)一步探究斜激波對高噴霧當(dāng)量比燃燒的增混燃燒效果, 基于斜激波角β=37°工況, 改變噴霧當(dāng)量比Φ0=3.6, 分析其在τ=220Δt時的無量綱溫度瞬時分布圖, 見圖50, 其中藍(lán)色虛線內(nèi)包含的區(qū)域表示燃料質(zhì)量分?jǐn)?shù)YF=[0.05,0.12]。

圖50 噴霧當(dāng)量比Φ0=3.6工況下無量綱溫度瞬時分布圖Fig. 50 Instantaneous distribution of T* in case Φ0=3.6

分析對應(yīng)的數(shù)值紋影圖, 見圖51, 其中紅色實線內(nèi)包含的區(qū)域表示化學(xué)反應(yīng)速率RF∈[5 000,20 000]kg/(m3·s)。

由圖50和51可知, 由于噴霧當(dāng)量比提高, 該工況旋渦內(nèi)燃料含量高, 提供波后燃燒的釋能高, 因此激波后的氣體壓強(qiáng)與化學(xué)反應(yīng)速率峰值更高。因此在富燃混合分?jǐn)?shù)條件下, 仍會出現(xiàn)較高速率的化學(xué)反應(yīng)過程。而且, 化學(xué)反應(yīng)形成于高標(biāo)量耗散區(qū)域, 這與以往亞聲速弱可壓縮流動中的研究結(jié)論不同, 這表明激波作用下的壓縮效應(yīng)較當(dāng)?shù)貥?biāo)量耗散更加重要。

圖51 噴霧當(dāng)量比Φ0=3.6工況下數(shù)值紋影圖Fig. 51 Numerical schlieren in case Φ0=3.6

(a) Ma, P*, ρ*and T* curves

(b) RF and P* curves

(d) RF and χ versus Z

繪制噴霧當(dāng)量比Φ0=3.6工況下在垂直坐標(biāo)y/δ0=-70處的相同流向區(qū)域內(nèi)的剖面, 比較無量綱壓強(qiáng)P*和化學(xué)反應(yīng)速率RF隨位置變化曲線, 即爆震波形成的時序過程, 見圖53, 圖中箭頭表示流向坐標(biāo)x/δ0增大的方向。在τ=190Δt時, 發(fā)生了熱自燃, 化學(xué)反應(yīng)速率峰值位置滯后于壓強(qiáng)峰值位置, 見圖53(a); 但由于衍射激波傳播方向與來流速度方向相反, 二者相對Mach數(shù)增大使衍射激波增強(qiáng), 壓強(qiáng)峰值不斷上升, 壓強(qiáng)峰值也誘導(dǎo)產(chǎn)生相應(yīng)的化學(xué)反應(yīng)速率峰值, 見圖53(b); 在τ=210Δt時, 前導(dǎo)激波與反應(yīng)面在空間上已經(jīng)非常接近, 該激波與化學(xué)反應(yīng)的耦合現(xiàn)象可以定義為局部“近似爆震燃燒”狀態(tài), 見圖53(c); 但隨著剪切旋渦的對流運(yùn)動, 富燃混合氣向下游流動, 來流中缺乏新鮮燃料氣, 因此壓強(qiáng)峰值所處混合分?jǐn)?shù)由富燃狀態(tài)轉(zhuǎn)向貧燃狀態(tài), 見圖53(d)和(e); 當(dāng)前導(dǎo)激波開始壓縮貧燃預(yù)混氣, 激波作用區(qū)的化學(xué)反應(yīng)較弱, 激波與反應(yīng)面發(fā)生分離, 見圖53(f)。分析發(fā)現(xiàn), 該過程中局部“近似爆震燃燒”的出現(xiàn)取決于當(dāng)?shù)丶げǖ膫鞑ヒ约鞍殡S旋渦運(yùn)動過程中燃料氣體的空間非定常、 非均勻分布。

3.4 第三Damk?hler數(shù)及局部準(zhǔn)等容燃燒過程分析

上述研究結(jié)果表明, 等容燃燒和等壓燃燒過程的本質(zhì)是壓強(qiáng)波傳播特征時間尺度和化學(xué)反應(yīng)特征時間尺度的差異。為此, 提出表征不同燃燒模式的第三Damk?hler數(shù)DaⅢ

式中,τa為壓強(qiáng)波傳播時間,τf為化學(xué)反應(yīng)時間,L為特征長度,a為當(dāng)?shù)芈曀?ρ為當(dāng)?shù)貧怏w密度,Yf為當(dāng)?shù)厝剂辖M分濃度,ωf為當(dāng)?shù)鼗瘜W(xué)反應(yīng)速率。

(a) τ=190Δt

(b) τ=200Δt

(c) τ=210Δt

(d) τ=220Δt

(e) τ=230Δt

(f) τ=240Δt

當(dāng)?shù)谌鼶amk?hler數(shù)DaⅢ>>1時, 燃燒過程產(chǎn)生的壓強(qiáng)波來不及向周圍環(huán)境傳播, 導(dǎo)致當(dāng)?shù)貕簭?qiáng)和溫度大幅提升, 產(chǎn)生了“壓強(qiáng)脈沖峰”, 該過程為準(zhǔn)等容燃燒過程; 當(dāng)DaⅢ<<1時, 燃燒過程產(chǎn)生的壓強(qiáng)波迅速向周圍環(huán)境傳播, 因此當(dāng)?shù)貕簭?qiáng)基本不變, 該過程即為等壓燃燒過程; 當(dāng)DaⅢ≈1時, 燃燒波和壓強(qiáng)波耦合, 壓強(qiáng)波仍可向周圍環(huán)境傳播, 同時也引起當(dāng)?shù)厝細(xì)獾捏w積膨脹并使當(dāng)?shù)貕簭?qiáng)和溫度有一定賦值的升高, 為介于等容和等壓之間的燃燒過程。

研究表明, 準(zhǔn)等容燃燒發(fā)生時, 其局部第三Damk?hler數(shù)一般可達(dá)等壓燃燒的50倍以上, 因此, 可以通過第三Damk?hler數(shù)判別流場中的燃燒模式。在燃燒流場中, 若化學(xué)反應(yīng)過程非常迅速, 那么第三Damk?hler數(shù)會在局部反應(yīng)區(qū)域發(fā)生急劇升高, 與“壓強(qiáng)峰”所對應(yīng), 可用于燃燒不穩(wěn)定性及其激勵機(jī)制的分析過程中。

張耘隆[34]研究超聲速混合層燃燒過程時, 成功利用第三Damk?hler數(shù)標(biāo)識局部發(fā)生的準(zhǔn)等容燃燒, 發(fā)現(xiàn)在超聲速混合層燃燒過程中, 當(dāng)著火距離(小于計算域長度)大于渦脫落距離時, 著火之前混合層的大尺度渦中蘊(yùn)含了燃料與空氣的良好預(yù)混氣, 致使后續(xù)可發(fā)生準(zhǔn)等容燃燒, 這是進(jìn)一步促使燃燒不穩(wěn)定發(fā)生的物理機(jī)理。尕永婧[42]采用第三Damk?hler數(shù)成功分析了液氧煤油發(fā)動機(jī)高頻燃燒不穩(wěn)定性機(jī)理。

3.5 小結(jié)

本節(jié)針對強(qiáng)可壓縮氣液兩相反應(yīng)流動, 構(gòu)建了基于Euler-Lagrange框架的高精度數(shù)值仿真軟件TURFsim, 以斜激波與超聲速混合層相互作用為例, 研究了剪切層流動中液霧顆粒的彌散及蒸發(fā)過程, 分析了極限條件火核生成和火焰的傳播模式, 數(shù)值展現(xiàn)了液霧彌散、 小激波、 局部爆震波的生成及其時空演化, 揭示了斜激波增強(qiáng)混合與增強(qiáng)燃燒的物理規(guī)律, 提出了新的無量綱參數(shù)——第三Damk?hler數(shù)DaⅢ, 用其定量表征局部準(zhǔn)等容燃燒, 實現(xiàn)了高瞬態(tài)兩相流與反應(yīng)流關(guān)鍵特征辨識及定量分析。

4 總結(jié)與展望

4.1 總結(jié)

本文總結(jié)了清華大學(xué)噴霧燃燒與推進(jìn)實驗室近幾年針對極高速、 強(qiáng)可壓、 高瞬變等極端條件下氣液兩相界面流與液霧反應(yīng)流, 在流體物理、 數(shù)值模型與高精度算法等方面的研究成果與進(jìn)展。

在高瞬變兩相流數(shù)理模型與數(shù)值計算方法方面, 以剛性氣體狀態(tài)方程中表征兩相流體物理屬性的變量作為界面變量, 通過建立表征界面量的輸運(yùn)方程描述物質(zhì)界面, 從而將物質(zhì)界面轉(zhuǎn)化為數(shù)學(xué)間斷, 統(tǒng)一了兩相界面流動的求解方法; 構(gòu)建了能夠同時捕捉復(fù)雜波系與兩相界面的變模板點WENO格式(WENO-IS), 實現(xiàn)了光滑區(qū)的高分辨率和大密度比界面的魯棒捕捉; 提出了基于化學(xué)勢弛豫的瞬態(tài)相變模型, 建立了基于壓強(qiáng)判斷函數(shù)的空化初生準(zhǔn)則, 實現(xiàn)了對高瞬變相變過程以及相變觸發(fā)復(fù)雜波系時空演化過程的數(shù)值模擬。

結(jié)合數(shù)值模擬與理論分析, 將激波在氣液兩相界面的反射與透射現(xiàn)象分為“常規(guī)反射+非前導(dǎo)型透射”“常規(guī)反射+前導(dǎo)型透射”和“Mach反射+前導(dǎo)型透射”3種類型, 提出了入射激波臨界角等判據(jù)。針對高速液滴撞擊壁面問題和激波液滴相互作用問題, 結(jié)合小波分析方法以及射線追蹤分析方法, 建立了一套完整的曲界面約束下非定常激波動力學(xué)理論解析方法。借此, 開展了高速液柱和壁面、 激波和液滴等相互作用過程中激波系時空演化的數(shù)值模擬與理論解析, 并推導(dǎo)了上述過程中曲界面匯聚膨脹波誘導(dǎo)液滴流體空化初生位置的理論公式。

針對強(qiáng)可壓縮氣液兩相反應(yīng)流動, 構(gòu)建了基于Euler-Lagrange框架的高精度數(shù)值仿真軟件TURFsim; 以斜激波作用超聲速混合層為例, 開展了數(shù)值模擬研究, 分析了激波和旋渦相互作用過程中旋渦和激波的演化規(guī)律, 獲得了斜激波作用下混合層首先被壓縮, 混合層厚度增長率升高, 進(jìn)而使最終混合層厚度增大的結(jié)論, 發(fā)現(xiàn)斜激波作用下混合層湍動能和流動脈動強(qiáng)度增強(qiáng), 進(jìn)而導(dǎo)致混合層整體混合效率增強(qiáng)的結(jié)論; 在超聲速液霧混合層反應(yīng)流動中, 分析了極限條件火核生成和火焰的傳播模式, 數(shù)值展現(xiàn)了液霧彌散、 小激波、 局部爆震波的生成及其時空演化, 提出了第三Damk?hler數(shù)DaⅢ定量表征不同燃燒模態(tài)的分析方法, 并成功應(yīng)用于局部準(zhǔn)等容燃燒過程的辨識與演化分析。

4.2 展望

極端條件兩相界面流與反應(yīng)流機(jī)理、 模型與算法研究仍有大量科學(xué)問題尚待解決, 面向工程應(yīng)用仍有大量關(guān)鍵技術(shù)亟待突破, 鑒于作者水平有限, 難以面面俱到和闡述深入, 但有必要在如下方面開展系統(tǒng)深入的研究。

1)在高速強(qiáng)可壓縮效應(yīng)及非定常激波方面, 進(jìn)一步通過實驗測試和數(shù)值模擬等方法獲得精細(xì)化流場結(jié)構(gòu)、 激波演化及智能化辨識, 并將其應(yīng)用于高超聲速飛行器氣動優(yōu)化、 飛發(fā)一體化建模和仿真, 開展超聲速燃燒、 爆震燃燒以及內(nèi)外流耦合作用機(jī)理與規(guī)律等研究工作。

2)在多相空氣動力學(xué)及快速流體相轉(zhuǎn)變與多物理耦合機(jī)理研究方面, 可進(jìn)一步開展激波、 爆震波在復(fù)雜介質(zhì)中傳播機(jī)理分析, 針對發(fā)動機(jī)結(jié)冰和跨介質(zhì)飛行中涉及的液滴撞擊、 液滴蒸發(fā)、 燃料液滴燃燒、 復(fù)雜界面和波的相互作用等開展研究工作。

3)高速強(qiáng)可壓復(fù)雜介質(zhì)的多物理耦合數(shù)理模型及算法研究方面, 可進(jìn)一步擴(kuò)展不同間斷類型的統(tǒng)一描述模型、 剛性源項的松弛求解器和精確、 魯棒化學(xué)反應(yīng)動力學(xué)機(jī)理等數(shù)學(xué)、 物理、 化學(xué)模型, 開發(fā)高效高精度數(shù)值格式、 自適應(yīng)網(wǎng)格算法、 任意多面體網(wǎng)格高精度算法和化學(xué)反應(yīng)機(jī)理調(diào)控與加速算法等, 并針對高性能(尤其E級)計算機(jī)開展算法與程序軟件應(yīng)用, 為工程應(yīng)用提供依據(jù)和指導(dǎo)。

致謝感謝國家自然科學(xué)基金(51676111, 91952205)對實驗室開展科研工作的資金資助, 感謝國內(nèi)外同行對研究工作的關(guān)注、 支持、 批評與指正, 這是實驗室全體科研人員勇敢探索、 不斷前進(jìn)的動力。