高馬赫數下激波液滴相互作用的數值模擬研究1)

宋家喜 潘書誠

(西北工業大學航空學院,西安 710072)

引言

液滴在高速氣流作用下發生變形和破碎的研究可追溯至20 世紀50 年代,主要應用于大規模化學武器的高速投遞、超聲速飛行中的雨滴損傷評估以及火箭發動機的燃油噴射.例如,超聲速飛行器在穿越大氣層時會遇到暴雨,其雨滴在高速氣流下的撞擊和破碎會對飛行器表面材料造成嚴重損傷[1].早在1949 年,Hinze[2]通過實驗確定了韋伯數(Weber number)為表征液滴變形和破碎過程的主要參數.此后,Hanson 等[3]以及Ranger 等[4]通過大量激波管實驗發現液滴的變形破碎與周圍剪切氣流形成的邊界層有很大關系,由此建立了液滴破碎的邊界層剝離(boundary layer stripping)模型.Patel 等[5]通過實驗研究了強表面張力下的液滴破碎行為,這種情形下Rayleigh-Taylor (RT)不穩定性波的穿刺(piercing)現象在液滴破碎中起著重要作用.另外,Pilch 等[6]對前人的研究工作進行了總結.根據不同的韋伯數,他們將液滴的破碎模態劃分為經典的五種破碎模態,分別為振蕩(vibrational) 破碎、袋狀(bag) 破碎、袋蕊(bag-and-stamen) 破碎、剪切剝離[7](stripping)破碎和毀滅(catastrophic)破碎.為了彌補傳統紋影顯示技術的不足,Theofanous 等[8-9]采用激光誘導熒光技術來觀測液滴破碎.其研究發現在高韋伯數下液滴的破碎是剪切誘導界面剝離的結果,而并非RT 波穿刺所導致.他們發現傳統實驗所觀察到的毀滅破碎機制并不存在,只是液霧遮蔽效應的一種誤導,從而對液滴破碎模態進行了重新劃分.根據氣動力相對于黏性力和表面張力的大小,他們將液滴破碎模態劃分為低韋伯數下的RTP (Rayleigh-Taylor piercing)破碎機制和高韋伯數下的剪切誘導剝離(shear-induced entrainment,SIE)破碎機制.當氣動力顯著小于黏性力和表面張力時,RTP 是主要的破碎機制.而隨著氣動力的增加,SIE 將成為主導的破碎機制.

相比于實驗測量,數值模擬能夠捕捉到更詳細的液滴界面和流場的演化細節,目前已成為研究液滴變形破碎的一種重要手段.但受限于目前計算能力,高分辨率的三維激波液滴相互作用數值模擬依然較難實現.前人的研究主要集中于激波沖擊液滴的二維數值模擬[10-11]以及三維對稱模擬[12-13].Kaiser 等[11]采用高分辨率的數值方法對二維液滴在高韋伯數下的變形破碎過程進行了數值模擬.結果表明,只有在足夠精細的分辨率下液滴界面演化過程中迎風面的帽狀結構才能得到充分解析.而黏性力和表面張力對于高韋伯數下液滴變形破碎的影響可被忽略.全三維的數值模擬在近幾年也有出現,Meng 等[14]對液滴在高韋伯數下的SIE 破碎機制進行了數值模擬,定性描述了液滴的剪切剝離破碎機制和周圍流場的演化等,包括前期液滴在激波沖擊下的扁平化過程以及中后期剪切層破碎的形成.他們對液滴質心位移、速度和加速度的統計揭示了液滴在激波沖擊下的加速推進[15]特征.另外,他們分析了液滴表面不穩定性,并對速度場進行傅里葉分解,試圖找出液滴破碎后期非對稱性形成的成因.隨后,Dorschner 等[16]對中等韋伯數和無窮大韋伯數下的液滴破碎進行了全三維數值模擬和實驗研究.他們發現液片的形成和脫落是一個反復循環的過程,Sharma 等[17]的實驗也發現了這一點.通過對流場的方位角傅里葉變換,他們觀察到了方位模態的發展,表面液滴破碎存在橫向方位調節機制.但由于他們采用的是耗散界面方法,無法獲取精確的液滴界面.因此,其數值模擬中觀察到橫向不穩定性的準確性還有待驗證.

另外,近幾年來有關其他參數對液滴變形破碎影響的研究也逐漸增多.Wang 等[18]實驗研究了在恒定韋伯數情況下,通過不同材料流體的相互組合研究馬赫數和雷諾數對液滴破碎的影響.他們的實驗表明,在其他參數不變的情況下,隨著來流馬赫數從0.3 增加到1.2,液滴的破碎樣式發生明顯變化.他們還通過實驗首次研究了兩個串聯液滴在激波沖擊下的變形破碎行為,比較了不同韋伯數和串聯液滴間距對液滴破碎行為的影響[19].Leung 等[20]設計了一款新型的激波管來研究高馬赫數下激波和液滴的相互作用,通過對比數值模擬結果研究了2 到5 馬赫數下激波沖擊液滴的變形破碎過程.Nykteri 等[21]用一種多尺度的兩相流方法數值模擬了軸對稱液滴在1.6 至2.64 的馬赫數區間的SIE 破碎機制,描述了液滴從早期到晚期的變形和破碎行為.Garcia-Magarino 等[22]通過對乙醇液滴破碎的實驗研究發現了一種新的破碎模態,這在以前的類水液滴破碎研究中還沒有報道過.然而,由于實驗條件的限制,極高馬赫數(比如馬赫數大于5 的高超聲速條件)下激波和液滴相互作用相關實驗的實現比較困難,與其對應的數值模擬研究也沒有出現.因此,有必要利用數值模擬的優勢開展極高馬赫數下激波和液滴相互作用的相關研究.

國內學者也開展了許多有關液滴變形和破碎的研究.陸守香等[23]在2000 年對激波沖擊下的液滴變形破碎進行了數值模擬,并提出了初始霧化時間的概念,分析了液滴的變形破碎特征.耿繼輝等[24]對液滴的變形破碎行為進行了較為系統的實驗方面的研究.實驗結果表明低韋伯數下液滴的變形破碎時間明顯更長,初始液滴形狀對液滴的變形破碎有著顯著影響.樓建鋒等[25]將VOF (volume of fluid)方法和湍流模型組合,數值模擬了相關實驗,計算結果表明韋伯數對液滴變形破碎起促進作用.楊威等[26]采用簡化的CLSVOF (coupled level-set and VOF)算法以及自適應網格加密技術對液滴在氣流中的破碎過程進行了數值模擬研究.其結論解釋了低韋伯數下液滴破碎和RT 不穩定性波的關系.另外,他們首次對低韋伯數下的液滴變形和破碎進行了三維數值模擬,解釋了此時RT 波的重要作用[27].朱萬里等[28]數值模擬了不同氣流壓力對液滴在SIE 破碎模態的影響,重點關注液滴破碎后小液滴的大小分布狀態,結果發現較高的氣流壓力會顯著增加液滴的變形.施紅輝等[29]在水平激波管中進行了超聲速條件下激波沖擊亞毫米液滴變形破碎的實驗研究,得到了各種韋伯數下液滴變形破碎的紋影圖像.中國科學技術大學的丁航教授團隊[30-31]發展了高精度的守恒清晰界面可壓縮多相流數值方法,并將其應用于激波沖擊液滴數值模擬研究.另外,有關黏性液滴[32-33]、非牛頓液滴[34]以及反射激波對界面不穩定性的影響[35]最近幾年也有學者研究.就目前的研究進展來看,國內對于液滴變形和破碎的全三維數值模擬還較為缺乏,有關各個無量綱參數對液滴破碎的影響以及液滴在推進和變形過程中的定量研究也較為缺乏.

本文的主要工作是在高韋伯數的前提下,針對高馬赫數和極高馬赫數下液滴變形破碎動力學開展數值模擬研究.根據Kaiser 等[11]的高韋伯數下二維液滴變形破碎的有關結論,液滴的表面張力和黏性力可以忽略不計.本文采用基于守恒清晰界面方法的可壓縮多相流數值模擬方法[36-39]來捕捉液滴界面的變形和流場演化,并使用自適應多分辨率方法來提高模擬精度和效率.其中,單元面通量采用5 階WENO 格式[40]進行重構,時間積分使用二階Runge-Kutta 格式[41].數值計算程序為文獻[37]中開發的自適應多分辨率多相流計算程序ALIYAH.本文的結構如下: 第1 章描述所研究的激波液滴相互作用問題的物理模型;第2 章說明所使用的多相流數值模擬方法;第3 章驗證數值模擬結果的準確性并分析高馬赫數下液滴推進變形和破碎的數值模擬結果.

1 物理模型

1.1 問題描述

本文的研究對象為單個液滴在超聲速和高超聲速氣流沖擊下的變形破碎行為.在數值模擬中,液滴周圍氣體的高速流動通過一道正激波產生.激波馬赫數定義為激波速度和聲速的比值,即

其中,us代表激波速度,a為聲速.液滴的變形破碎過程受慣性力、黏性力和表面張力的共同作用.慣性力驅使液滴發生變形最終破碎,黏性力阻礙液滴的變形,表面張力促使液滴保持原始的球形形狀.在液滴的變形破碎過程中,韋伯數We和奧內佐格數Oh為描述該物理過程最重要的兩個無量綱參數.We定義為液滴受到的慣性力和表面張力的比值,即

其中,ρg和ug分別為激波后的氣流密度和速度.D0代表液滴的初始直徑,σ=0.073 N/m 為液滴的表面張力系數.Oh定義為液滴受到的黏性力和表面張力之比,即

式中,μl和 ρl分別代表液滴的黏性系數和密度,D0代表液滴的初始直徑.其中,本文算例中的液滴黏性系數為1 mPa·s,密度為1 t/m3.

物理計算域如圖1 所示,沿流向取6 倍的液滴初始直徑,沿展向和垂向都為4 倍的液滴初始直徑.這樣的計算域大小既避免了液滴受到壁面反射的影響過大,又不至于因為計算域太大而導致過高的計算量.首先模擬384×256,768×512 和1536 ×1024 三種有效網格分辨率下的二維算例.在充分平衡計算所需代價和數值模擬效果之后,最終選擇使用分辨率為768×512×512 來進行全三維液滴破碎的數值模擬,這部分的內容會在3.1 節進行詳細分析.與Meng 等[42]的研究類似,這樣的分辨率相當于沿液滴直徑布置128 個網格.如圖1所示,為了真實模擬實驗中的激波管壁,計算域前后上下面均設為對稱邊界條件,而左右界面分別為入流和出流邊界條件.液滴的初始形狀為直徑3.5 mm的球形.其中心位于矩形計算域的中軸線距左邊界5.5 mm 處,而激波初始位置位于液滴中心前2.0 mm處.液滴初始密度 ρl=1 t/m3,激波前空氣密度為1.2 kg/m3,激波前空氣壓為100 kPa.根據Rankine-Hugoniot 激波關系式計算得出不同激波馬赫數Ms的流動參數如表1 所示.

圖1 激波驅動液滴破碎數值模擬的計算域Fig.1 Computational domain for shock-driven droplet breakup

表1 不同激波馬赫數下的激波后流動狀態Table 1 Flow conditions behind various Mach air shocks

1.2 控制方程

本文研究液滴破碎所使用的控制方程為

其中U=(ρ,ρu,E)T為守恒量,F=(uρ,ρu?u+pI,u(E+p))T為對流通量,Fv=(0,T,T·u)T為黏性通量,S代表考慮表面張力的源項.t為物理時間,ρ 代表密度,u為速度矢量,T為黏性應力張量.E為總能量,p代表壓強.控制方程需要加入狀態方程得到封閉,對于空氣,使用理想氣體狀態方程

而水的狀態方程采用Tait 狀態方程

其中,γ=7.15,A=100 kPa,B=331 MPa,ρ0=1 t/m3.

2 數值方法

2.1 基于守恒清晰界面模型的可壓縮多相流數值方法

在笛卡爾網格下將控制方程在計算域 Ω 內使用有限體積方法進行離散,計算域 Ω 被隨時間演化的兩相界面劃分為液相域 Ω1(t) 和氣相域 Ω2(t).圖2 給出了二維笛卡爾網格下的清晰界面方法示意圖,其中 ΔΓi,j(藍色直線) 代表每個網格單元內對精確界面(紅色曲線)的分段線性近似.對控制方程在每個子域 Ωn(t) 和網格單元 Δi,j的交集 Ωn(t)∩Δi,j內進行積分并應用高斯定理可得,

圖2 二維切割網格下的守恒離散示意圖Fig.2 Two-dimensional schematic of the conservative discretization in a cut cell

其中 Γ(t) 為每個網格單元隨時間變化的由level-set函數[43]確定的分段界面. ?Δi,j為每個網格單元的邊界.n為網格單元邊(2D)或面(3D)的單位法向量,而nΓ代表界面的單位法向量.Ωn(t)∩Δi,j的體積為αi,jΔi,j,其中αi,j為體積分數,? (Ωn(t)∩Δi,j) 在網格(i,j)內可以用近似,其中A如圖2 所示為?Δi,j被ΔΓi,j切割后的剩余量.這樣式(7)可表達為

將上式擴展到三維可得

其中,Δt為時間步長,代表n時刻網格單元(i,j,k)的單元平均狀態量,Fi,j,k為單元的界面通量,X(ΔΓi,j,k)為氣液界面交換項.對式(9)在單相流場上進行求和可得

右手邊的第二項為界面交換項,只存在于被切割的網格單元.對整個計算域的兩相流體來說,界面交換項的符號相反.對整個兩相流體進行求和可得

這說明對兩相流體來說,本方法在全局上是守恒的.更多關于守恒性的驗證可參考本數值方法的文獻[36].最大時間步長根據CFL 數來確定,本文中使用的CFL 數均為0.5.兩相界面通過水平集函數φ(x)來描述,即 φ(x)=0 代表兩相界面,φ(x)＜0 代表流體1,φ(x)＞0 代表流體2,|φ(x)| 為單元網格x的中心坐標距離界面的距離.水平集函數的演化用以下的對流方程描述

其中,uφ代表水平集對流速度.在被界面切割的單元里面,uφ等于界面速度uΓ,uΓ可以通過在相界面求解一個兩種物質的黎曼問題來得到[36].平均曲率κ為界面法向矢量的散度,而界面法向矢量nΓ可以通過以下求解得到

在未被界面切割的單元中,水平集對流速度等于外推得到的界面速度可由擴展方程(extension equation)的穩態解得到

為了維持水平集函數的符號距離特性 |?φ|=1,需要對水平集函數進行重新初始化[44]

這里 φ0指重新初始化開始之前的水平集場.為了獲取界面附近網格點的流體狀態量,通過虛擬流體法(ghost fluid method)[45]將界面一側流體狀態量外推到另一側

其中,q代表需要外推的流體狀態量,而N≡(Nx,Ny,Nz)是根據水平集函數確定的界面法向量.另外有關單元孔隙A、單元體積分數 α,以及界面交換項X的推導細節可以參考圖2 和文獻[36-39].

為了統計液滴推進和變形,本文用以下公式計算液滴質心位移xc、速度uc和加速度ac

其中,αl,x,ρl,V,u分別為網格單元的液相體積分數、流向坐標值、密度、單元格體積和速度.液滴質心位移、速度和加速度分別使用下式進行無量綱化

式中,x0為液滴初始位置的球心坐標,D0為液滴的初始直徑,ug為激波后的氣流速度.另外,本文所提到的液滴破碎時間均為無量綱時間,使用文獻[4]中提出的公式對液滴破碎時間進行無量綱化,其中

式中,ρg和 ρl分別為激波后的空氣密度和液滴密度.

2.2 基于塊結構的自適應多分辨率網格加密方法

在液滴破碎的過程中,為了精確捕捉到液滴界面的變形,需要在液滴界面和激波附近布置足夠精細的網格,而其他地方使用較粗的網格以節省計算量.為此,應用一種基于塊結構的自適應多分辨率方法[37,46],這種方法由動態數據結構和自適應更新過程組成.在自適應更新過程中,網格的粗化和細化過程通過投影算子Pl+1→l和預測算子Pl→l+1來實現,以一維結構網格為例,兩者分別定義為

這樣,使用投影算子Pl+1→l從l+1 層網格單元計算得到l層網格單元的過程是精確的,而利用預測算子Pl→l+1根據l層網格單元來預測l+1 層網格單元的值會產生一定的預測誤差.l層第i個網格單元的預測誤差可以定義為在本文的工作中,采用五階插值來定義預測算子,即m=2,相應的系數為

當預測誤差大于某一個閾值時,相應的單元或塊結構將會被細化,這個閾值定義為

其中,D0為空間維度,Lmax為自適應數據結構的最大層數,ε 為參考誤差值.另外,始終使用最密的網格捕捉界面.這樣,如圖1 所示,在界面和流場劇烈變化(如激波)區域加密網格.

3 結果分析

3.1 網格分辨率的影響

完全解析激波和液滴相互作用問題中所有的物理細節需要的計算量對于目前計算機很難承受.根據Popinet[47]的研究,完全解析液滴表面張力影響所需的網格尺度可按照單元格韋伯數來進行估算

按照這樣的估算,僅僅解析Ms=3 下的相關算例就需要沿液滴直徑布置2.2×106個網格,這樣的全三維數值模擬所需的網格量達到了驚人的 1021,這遠遠超出了目前超級計算機的計算能力.此外,根據Chang 等[12]的研究,為了捕捉液滴破碎過程中的KH 不穩定波,他們使用了八層自適應網格,相當沿液滴直徑布置800 個網格.對于黏性邊界層的捕捉,另外還需在液滴表面沿直徑布置4000 個網格,這樣的計算量也超出了絕大多數計算機的計算能力.因此,為了使計算量可以承受,本文把研究的主要關注點放在液滴界面的變形、液滴加速過程的量化以及流場結構的演化上.首先對不同網格分辨率下的結果進行比較以驗證計算結果的收斂性.隨著計算網格的不斷細化,數值模擬可以捕捉到更加精細的流場結構和界面變形的細節,但隨之增加的計算代價是巨大的.因此,需要選擇一個合適的網格分辨率,既能捕捉到相對完整的界面變形特征,又能盡可能地節省計算量.

如圖3 所示,在激波馬赫數為3 的情況下三種不同網格分辨率的二維液滴破碎數值紋影圖表明液滴的界面變形和流場結構具有相似的演化趨勢.顯然,分辨率越高,捕捉到的液滴界面和流場結構越精細.在最粗的384×256 分辨率下,盡管液滴的整體變形特征與768×512,1536×1024 分辨率下類似,但界面結構特征的細節還不夠完善.液滴上游的帽狀結構沒有捕捉到,赤道區域的剪切層也未捕捉到.對比768×512 和1536×1024 下的界面變形和流場結構可以發現,768×512 分辨率已經能夠捕捉到相對完善的界面特征和比較精細的流場結構.

圖3 不同網格分辨率下的數值紋影圖Fig.3 Numerical schlieren images for various grid resolutions

圖4 統計了三種網格分辨率下液滴質心的無量綱位移和無量綱速度隨時間演化規律.384×256 分辨率下的液滴質心運動曲線和另外兩者在后期有顯著差異,其可能原因是由于加密后的網格能夠解析出從原始液滴剝離形成的薄液片和小液滴,而小液滴的存在使得液滴整體的質心變得靠后,導致384 ×256 分辨率下液滴后期的位移和速度明顯低于其余兩者.由于768×512 和1536×1024 下質心運動曲線之間的差異已經相當的微小,為了節省計算資源,后續采用768×512×512 的分辨率來進行液滴破碎的全三維數值模擬.這樣的網格分辨率足以捕捉到相對完整的界面演化狀態和液滴的加速運動規律,相對而言計算量也在可接受的范圍之內.

圖4 不同網格分辨率下的無量綱質心位移和速度演化Fig.4 Evolution of the dimensionless center-of-mass drift and velocity for various grid resolutions

3.2 液滴表面張力和黏性的影響

一方面,在高We和低Oh下,氣動力為驅使液滴變形破碎的主導機制.此時液滴的表面張力和黏性力在抵抗液滴破碎的過程中所起到的作用相當微弱,表面張力和黏性力在量級上是可以忽略不計的.另外一方面,目前的網格分辨率根本達不到完整捕捉液滴表面張力造成的毛細效應的影響.

為了驗證表面張力和黏性作用,在768×512 的網格分辨率下,分別模擬了在忽略表面張力和黏性力以及考慮表面張力和黏性力2 種條件下的二維液滴破碎行為.如圖5 所示,在激波馬赫數為3 的條件下,表面張力和黏性力對液滴界面變形的影響甚微,僅僅從尾渦的充分發展區域才能發現微小的差異.另外,如圖6 所示,兩種情形下的液滴質心位移、速度和加速度曲線幾乎完全一致.唯有液滴破碎中后期的加速度曲線有所出入,這一出入更多是數值誤差引入的振蕩,并非由表面張力和黏性力導致.經過如上對比,在高We、高Ms和低Oh下,表面張力和黏性模型可以被忽略,故在進行液滴的全三維數值模擬過程中忽略表面張力和黏性力的影響,這也與Kaiser 等[11]的做法一致.

圖5 表面張力和黏性力對激波沖擊液滴數值紋影圖的影響Fig.5 Effects of capillary and viscous forces on Numerical schlieren images for shock-droplet simulations

圖6 表面張力和黏性力對液滴無量綱質心位移、速度和加速度演化的影響Fig.6 Effects of capillary and viscous forces on evolution of the dimensionless center-of-mass drift and velocity and acceleration for shock-droplet simulations

3.3 和實驗結果的比對

圖7 和圖8 為不同視角下Ms=3 的激波沖擊液滴變形破碎的界面演化圖,上面一行為數值模擬結果,下面一行為相應無量綱時間下的實驗結果.實驗數據來源于Theofanous 等[9]的實驗,其激光誘導熒光可視化技術顯示了液滴在激波沖擊下的變形破碎過程.實驗所使用的液滴直徑和數值算例均為3.5 mm,并且激波馬赫數均為3.盡管實驗使用的液滴黏性和表面張力與水不同,但其SIE 機制是相似的.值得注意的是,原始實驗數據的激波運動方向是從右向左的,而本文所有數值模擬結果的激波運動方向均為從左至右.為了方便和數值模擬結果做比較,將實驗結果的圖片旋轉 180°.

圖7 數值模擬結果(上)和實驗可視化(下)側視圖對比Fig.7 Comparison of numerical results (upper) and experimental visualizations (lower) for side view

圖8 數值模擬結果(上)和實驗可視化(下)前側30°視圖對比Fig.8 Comparison of numerical results (upper) and experimental visualizations (lower) for 30° front side view

從圖中可以看出,數值模擬得出的液滴破碎形態、破碎時間都和實驗結果吻合,表明該方法可以有效捕捉液滴破碎過程中的界面變形.激波沖擊過后,液滴在高速氣流的作用下發生迅速的扁平化.在扁平化的過程中液滴迎風面的前駐點區域一直保持光滑,臨近赤道區域出現凸起.同時液滴背風面逐漸變平,赤道區域在氣流剪切作用下拉伸出薄的液層.在實驗和數值結果中均可觀察到這些SIE 破碎機制的顯著特征.當然,目前所允許的分辨率所計算出的結果還不足以捕捉到SIE 破碎機制的全部特征,比如實驗觀察到的液滴背風面在扁平化的過程中會出現帶狀條紋.據Chang 等[12]的研究,這是KH 不穩定性在液滴表面不斷發展的結果,當然目前的分辨率還不足以捕捉到液滴表面的KH 波.而液滴赤道區域被拉伸的液層在早期就已經出現破碎形成小的細霧,這些破碎特征在數值模擬中可能需要更高的分辨率才可以捕捉到.

3.4 界面的演化及流場分析

圖9 為激波沖擊液滴前期(t*=0.029 48) 的對稱面上的數值紋影 lg(|?ρ|+1) 云圖,顯示了流場的波系結構演化.入射激波到達液滴的上游面會形成反射和透射.如圖9 所示,由于液滴的密度和聲速遠遠大于周圍氣體,為了在界面處滿足力學平衡[48],液滴迎風面會形成一道向上游傳播的反射激波,而一部分入射激波穿透液滴表面形成透射激波.

圖9 Ms=3 的激波沖擊下流場的早期數值紋影圖Fig.9 Numerical schlieren images for Mach 3 simulations at early stage

由于液滴的可壓縮性遠小于周圍空氣,液滴內部透射激波的運動速度要快于周圍入射激波的速度.透射激波在撞擊到液滴下游界面處再次發生反射,形成膨脹波.膨脹波向液滴上游傳播,由于液滴下游界面存在曲率,膨脹波在液滴內部發生匯聚,形成一個低壓區[49].隨著入射激波和兩相界面的夾角超過臨界值,規則反射向馬赫反射轉變,形成馬赫桿[50].馬赫桿在液滴的下游駐點處匯聚,形成二次波結構.液滴赤道處脫落的微霧在一對反向渦的作用下形成了二次循環區[49].如圖10(a)所示,液滴的變形程度并不明顯,而隨著二次波系統和二次循環區的形成(見圖10(b)～圖10(f)),在周圍氣流的影響下液滴開始發生扁平化,受剪切作用,液滴赤道處液片受到拉伸變薄,隨后發生破裂.

圖10 Ms=3 的激波沖擊下流場的數值紋影圖Fig.10 Numerical schlieren images for Mach 3 simulations

圖11 中的液滴流場壓力云圖表明隨著入射激波掃過液滴的后緣,周圍氣流壓力的不均衡導致了液滴的推進和變形行為.從圖11(a)可以看到,液滴的前駐點區域壓力最大,而后駐點區域壓力較小,液滴在此壓力場的作用下開始向下游推進.相比于前后駐點處的壓力,液滴赤道區域的壓力明顯更低,這樣的壓力差導致了液滴的扁平化行為.圖12 中的渦量云圖,氣流流過液滴最大直徑處,在液滴的背風面產生了分離,分離的氣流形成了一對明顯的渦.如圖12(b)所示,隨著主渦不斷向下游擴展,背風面的流動分離處又產生了一對與主渦方向相反的二次渦,如圖12(c)和圖12(d),這兩對渦的相互作用導致了液滴尾流區域的形成,也導致了液滴赤道區域的剪切剝離過程.

圖11 Ms=3 的激波沖擊下流場的壓力云圖Fig.11 Pressure contours for Mach 3 simulations

圖12 Ms=3 的激波沖擊下流場的Z 軸方向的渦量云圖Fig.12 Z-Vorticity contours for Mach 3 simulations

3.5 液滴的循環破碎過程

液滴破碎在更高的馬赫數下更加劇烈和清晰,因此通過Ms=11 情形的模擬結果來分析液滴的破碎過程.圖13 中顯示了液滴在Ms=11 的激波沖擊下的前側 45°界面演化圖,為了方便觀察薄液膜和小液滴，采用光線追蹤技術渲染模擬結果.如圖中的箭頭和虛線所標識,液滴在激波驅動下的整個破碎過程包含了三個重復的破碎模式.液滴在激波的沖擊下從初始球形破碎為小液滴并非只經歷了一次性的SIE 破碎,而是逐階段發生的,每次都會經歷相似的破碎過程.首先,液滴的赤道處受到剪切拉伸變為薄的液片,繼續拉伸導致液片較薄的地方出現破損而形成破洞.破洞逐漸擴展合并形成韌帶,最終韌帶的變細以及其互相撞擊導致破碎為小液滴.

圖13 Ms=11 下激波沖擊液滴破碎的光線追蹤渲染(前側45°視圖).第一次到第三次破碎分別使用箭頭、紅色虛線和藍色虛線來標識Fig.13 Ray-traced rendering of shock-droplet breakup for Mach 11 simulation (45° front side view).The first to third breakup are identified by arrows,red dotted lines and blue dotted lines,respectively

圖13 分別對每次液滴破碎做了標記,箭頭、紅色虛線和藍色虛線分別代表了液滴第一次、第二次和第三次的破碎過程.液滴下一次的破碎并非是在上一次破碎全部完成之后才開始的,其之間存在重疊.例如,在t*=0.103 6 時刻,第二次破碎剛剛進行到破洞形成的階段,此時液滴赤道處就已經拉伸出第三次破碎前期階段的薄液片了.這種重復循環的破碎機制會一直持續到液滴完全破碎為小液滴以至于主液滴已經無法在赤道拉伸出薄液片為止.雖然這種循環破碎在本文記錄的時間范圍內只重復了三次,但其是周期發生的.在最近有關液滴破碎的相關研究中,同樣發現了這種循環破碎的現象[16-17].在Dorschner 等[16]的實驗中,他們發現液滴從破碎開始到結束一共經歷了四次循環的破碎過程.同時,他們還測量了每次破碎發生的時間,破碎時間定義為液片從主液滴開始脫落的時間,也可以理解為薄液片出現破洞的時間.在韋伯數趨于無窮大的情況下,液滴第一到第三次破碎的時間分別為t*≈0.5,t*≈1.0和t*≈1.6,對應圖13(c)、圖13(f)和圖13(h).這驗證了本文的數值模擬結果和Dorschner 等[16]的實驗結果是相吻合的.根據他們的推斷,這種逐次循環破碎行為可能和變形后液滴的渦脫落有關.

3.6 馬赫數對激波沖擊液滴的影響

一直以來,We和Oh被認為是控制液滴變形最重要的兩個無量綱參數,故以往大多數學者都會把研究重點放在不同韋伯數下液滴變形和破碎的相關研究上,而關于其他無量綱參數對于液滴變形破碎的影響卻極少被研究.為了填補這一空白,利用數值模擬相對于實驗研究更易實現高馬赫數環境這一特點,本文研究了高馬赫數情形下馬赫數對液滴變形和破碎的影響.根據Theofanous 等[9]的實驗研究,高韋伯數下SIE 破碎機制為液滴的主導破碎機制.為了驗證這一結論,對比Ms=3 ,Ms=6 和Ms=11 三種馬赫數下液滴的變形破碎過程.

如圖14 所示,使用無量綱時間將3 種馬赫數下液滴的界面演化過程縮至同一時空尺度下.結果表明不同馬赫數下液滴的界面演化形態并無二致,基本可以歸為兩個階段: (1)首先,t*=0～0.3 時,液滴在周圍氣流壓力差的作用下發生扁平化行為,液滴赤道區域的剪切剝離還未開始,此時液滴基本上可以看作一個剛體,質心不會出現過大的偏移;(2)然后,在t*=0.3～0.6 階段,隨著液滴赤道區域剪切剝離的開始,液滴扁平化和剪切剝離同時進行,液滴逐漸發生破碎行為.

圖14 不同馬赫數下的液滴界面演化過程Fig.14 Evolution of the droplet interface for various Mach air shocks

雖然三種馬赫數下這兩個階段基本一致,但液滴的局部變形形態有所差異.在液滴扁平化的過程中,Ms=3 的液滴迎風面區域基本保持光滑,而背風面在t*=0.2 時就開始變平,之后基本都處于這種狀態.隨著馬赫數的提高,液滴背風面的扁平化程度有所下降,根據圖11 的壓力云圖可以定性分析這一現象產生的原因.在Ms=3 的情況下,液滴迎風面始終處于一個壓力較高的狀態,如此較為均勻的壓力分布使得液滴的迎風面得以保持光滑.隨著氣流在液滴的赤道處發生分離,液滴的赤道偏向背風面處形成了一個低壓區.在液滴的尾流區,馬赫桿的匯聚在液滴的尾流區形成了一個高壓區,如圖11(a)所示.液滴的背風面在尾流高壓區和赤道低壓區的雙重作用下被壓扁變平.然而,在Ms=6 尤其是Ms=11 的算例中,由于激波的速度太快,馬赫桿匯聚形成的高壓區存在時間較短,液滴背風面扁平化不明顯.在剪切剝離的過程中,高馬赫數所受的剪切力更大,液滴較早發生剪切剝離,剪切剝離出的液層更薄.

接下來使用液滴流向直徑和初始直徑之比Dcro/D0以及橫向直徑和初始直徑之比Dwise/D0來對液滴的變形過程進行定量研究.如圖15 所示,液滴流向直徑隨時間以近似線性的趨勢減小,而橫向直徑的增長表現出明顯的非線性.仔細觀察,液滴的橫向膨脹可分為兩個線性階段.在液滴的扁平化前期,即t*=0～0.3階段,液滴橫向直徑緩慢增長.而到了中后期的t*=0.3～0.6 階段,由于赤道周圍受到強烈的剪切拉伸,液滴沿橫向迅速膨脹.

圖15 不同馬赫數下無量綱流向直徑和無量綱橫向直徑的演化Fig.15 The evolution of the dimensionless cross-stream and the dimensionless streamwise diameter at various Mach air shocks

圖16 對比了液滴質心的位移、速度、加速度以及拽力系數的演化曲線.液滴的質心位移、速度及加速度使用式(18)進行計算和無量綱化,而拽力系數[43]定義為

圖16 不同馬赫數下的無量綱質心位移、速度、加速度以及拽力系數演化Fig.16 The evolution of the dimensionless center-of-mass drift and velocity and acceleration and drag coefficient at various Mach air shocks

其中m為液滴的總質量

式中,ac為液滴質心加速度,ρg和ug分別為波后氣體密度和速度,uc為液滴質心速度,Fd為液滴沿流向受力的總和.在液滴位移x*曲線中,三種馬赫數下液滴的位移曲線的差異很小.然而,液滴的無量綱速度u*和加速度a*曲線卻有所差異.隨著激波馬赫數的增加,液滴的無量綱速度和加速度也相應增大,這樣的結果與Meng 等[43]在液滴的二維破碎模擬中所得到的結果一致.分析3 種馬赫數下u*和a*曲線的演化,可以發現,隨著馬赫數的增加,液滴的無量綱速度和加速度的差異在逐漸減小.這表明在高超聲速的情況下,馬赫數對液滴運動的影響在逐漸減弱.

觀察圖16 中液滴無量綱速度和加速度曲線的差異,結合式(17)對無量綱參數的定義,這樣的定義顯然引入了密度和速度的分布.考慮到液滴沿著推進方向的運動其實在各個地方的速度是不均勻的,式(17)的定義本質上是液滴各微團推進速度和加速度的平均值.其實,還可以直接對液滴質心的位移曲線進行一次和二次求導來定義速度和加速度.這種情況下認為液滴是一個無形狀的點.如圖17 所示,這樣求出的三種馬赫數下的無量綱速度曲線基本吻合.這樣做雖然可以得到不變性,但是相比于式(17)并不能更物理地表征液滴的變形,因為液滴并不是以相同外形和均勻速度進行推進的.另外,因為液滴質心位移曲線是根據積分表達式計算出的一系列離散的數值,而非解析表達式,對質心位移曲線求導時采用的差分離散近似會產生數值振蕩過大的問題.在Meng 等[43]的研究中也討論了這樣的做法,這樣求出的液滴速度加速度曲線存在較大的噪聲誤差(如圖17 右).采用式(17)的定義能夠將該噪聲震蕩最小化.

圖17 直接對液滴質心位移求導得出的無量綱速度和加速度演化Fig.17 The evolution of the dimensionless velocity and acceleration derived directly from the dimensionless center-of-mass drift

參考液滴的拽力系數曲線可以發現,不同馬赫數下液滴的拽力系數曲線除了在前期峰值和后期振蕩區有明顯差異,其他地方吻合得特別好.分析液滴的速度曲線發現速度的斜率并非是一個常數,即液滴并非是以一個恒定的加速度來運動的.正如前面所分析的,液滴的變形破碎過程可分為兩個階段,在t*=0～0.3的這一階段,液滴的變形不明顯,此刻液滴還可以看作一個剛體,以一個相對較小的恒定加速度運動.隨著液滴扁平化的不斷發展,液滴的迎風面積增大,所受到的氣動力不斷增加,故液滴的加速度也不斷增大.隨著剪切剝離的開始,液滴的加速度也開始出現較大的震蕩,通過液滴的無量綱加速度曲線a*可以很清晰地看到這些振蕩規律.Meng 等[14]的研究認為振蕩的產生可能是旋渦脫落導致的,具體的振蕩機制還需要采用更高分辨率進行更加深入的研究.

4 結論

在高韋伯數的前提下,本文對激波馬赫數Ms從3 增加到11 時激波和液滴的相互作用進行了全三維數值模擬研究.求解得到的流場和液滴界面演化過程表明高馬赫數(Ms=3)和極高馬赫數(Ms＞5)下的液滴變形和破碎遵循統一的SIE 破碎機制.液滴的變形和破碎可分為液滴扁平化和剪切剝離兩個階段.其中液滴扁平化是前后駐點高壓區和赤道周圍低壓區的壓力差所導致的,而剪切剝離的產生原因是流動分離導致液滴背風面形成漩渦結構.另外可以發現在激波驅動下液滴從球形主液滴破碎為小的子液滴會經歷多個循環重復的破碎階段,這種逐次的循環破碎行為可能和變形后液滴的渦脫落有關.模擬結果也表明不同馬赫數下激波沖擊液滴所導致的變形和加速特征具有相似性.液滴的變形階段都遵循具有上述特征的SIE 破碎機制,而液滴的加速過程均展現出了非恒定的運動規律.在前期液滴橫向直徑變化不大的情況下,液滴可看作剛體以一個恒定的加速度運動,隨著液滴扁平化的發生,液滴的加速度也隨之增大.基于本文工作,未來研究將考慮極高馬赫數條件下液滴周圍空氣非平衡效應(比如)對激波誘導液滴變形破碎產生的影響.