壓縮性對渦環物理特征及其傳播速度的影響規律

林海燕， 向 陽， 張 斌， 劉 洪

(上海交通大學 航空航天學院， 上海 200240)

渦結構廣泛地存在于自然界及工程應用中，由于渦環是3維流場中最基本的渦系結構，眾多學者開展了對其特征性能的研究.Shariff等[1]指出研究渦環可以揭示渦動力學的一些物理特性及基本原理.針對不可壓縮渦環，研究者從渦環動力學的特征[2-3]、極限生長原理[4-5]、多渦環相互作用[6]等多個方面開展了大量研究.但是，相比于不可壓縮渦環，可壓縮渦環的研究十分有限且大多局限于研究可壓縮渦環的流動結構特征，如針對渦環夾止和尾跡射流[7]及由激波氣泡相互作用產生的渦環的生長演化過程[8]的研究，幾乎均未揭示可壓縮性對渦環物理特征的影響.對于可壓縮流動，一般而言可以分為亞聲速、跨聲速和超聲速3個區域，其流場存在不同的壓縮性特點.不同的壓縮性會使可壓縮渦環產生不同的結構特征.

一般采用激波管產生軸對稱的可壓縮渦環，而通過調節激波管驅動段和被驅動段之間的壓力比可以改變入射激波的馬赫數(Mashock)，即激波到達激波管出口時的馬赫數，進而得到具有不同可壓縮性的渦環.Elder等[9]在1952年首次通過紋影實驗觀測到了由激波管產生的可壓縮渦環向下游傳播的過程.Arakeri等[10]利用粒子圖像測速(PIV)技術研究當Mashock=1.1，1.2，1.3 時的可壓縮渦環，發現渦環完全形成后的傳播速度是激波到達激波管出口時波后速度的0.7倍.Dora等[11]同樣通過PIV技術研究Mashock=1.27，1.37時的可壓縮渦環，獲得了渦環結構參數的定量數據，并驗證了由實驗測得的渦環傳播速度與理論計算公式求得的結果相符.此外，部分學者還關注了可壓縮渦環結構中存在的激波.Baird[12]通過差分干涉測量法對可壓縮渦環進行研究的結果表明，當Mashock=1.5時，渦環的徑向區域存在嵌入后向激波，但該激波在強烈黏性剪切的作用下終止，并不能傳播至渦核中心.Brouillette 等[13]通過電花影圖和紋影拍攝技術研究不同條件下激波管產生的可壓縮渦環，發現能夠產生嵌入激波的Mashock的最小值為1.34.進而，Brouillette等[14]根據激波的結構特征對可壓縮渦環進行了更細致的分類：利用紋影拍攝技術發現在臨界驅動長度條件下，當Mashock<1.43時，將會產生沒有嵌入激波的渦環；當Mashock≥1.43時，則產生帶有嵌入激波的渦環；當Mashock=1.6時，嵌入激波的逆壓梯度將產生二次渦，這種二次渦被稱為反向旋轉渦環.隨后，又有學者針對反向旋轉渦環開展了更為細致的研究，如：Kontis 等[15]通過高速紋影拍攝技術觀察到了當Mashock=1.63時，主渦環下游存在多重反向旋轉渦環；Thangadurai等[16]通過石蠟油煙霧顯示技術同樣觀測到了反向旋轉渦環并詳細地研究了其形成機理及演變過程.

上述研究結果均表明，對于可壓縮渦環，可壓縮性的強弱對渦環流動結構的影響有很大的不同.然而，有關渦環自身的可壓縮性描述及可壓縮性對渦環自身結構參數影響的較為系統的研究尚未見報道.激波結構的存在與否是一個可對渦環結構的穩定性產生重要影響的因素，而渦環自身的結構參數(如渦核結構，渦環半徑，傳播速度)對于研究由渦環結構引起的夾帶、混合、燃燒等現象均十分重要.本文通過有限體積法對激波管產生的可壓縮渦環進行數值計算，以研究渦環自身的可壓縮性及其在不同可壓縮性條件下對渦環半徑、渦核半徑、渦環傳播速度以及渦核內渦量分布的影響.

1 研究方法

1.1 計算方法和模型

采用課題組自行研發的軟件進行數值計算；采用有限體積法對3維可壓縮Navier-Stokes方程進行空間離散；采用5階WENO(Weighted Essentially Non-Oscillatory)格式重構對流項；時間推進采用3階Runge-Kutta法；針對激波管產生渦環這個物理現象，為了方便驗證，采用與文獻[17]相同的物理模型.激波管計算模型示意圖及各組成部分的結構尺寸如圖1所示.激波管由驅動段A和被驅動段B組成，內徑為64 mm；C和D為環境區域，即渦環的形成區域及運動區域；被驅動段和環境區域的初始壓力為101.325 kPa；整個計算區域的初始溫度為300 K；改變驅動段的壓力即可得到不同驅動段與被驅動段的壓力比(Rp)，進而得到不同的Mashock.定義入射激波到達激波管出口時時間t=0；激波管出口處x=0.

通過改變渦環生長演化最主要區域(0 mm

表1 網格信息

為了驗證計算方法的正確性，在Rp=10及驅動段長度L=230 mm的條件下進行計算TS1-1，并將計算結果與文獻[18]中相同條件下算例的計算結果(TS2)進行比較，渦環生長過程中不同時刻渦環截面密度梯度云圖對比結果如圖3所示.由圖3可見：TS1-1與TS2的渦核中心位置以及嵌入激波位置均較為一致(見圖3(a))；TS1-1與TS2在對稱軸附近產生的滑移線位置也較為一致(圖3(b))；TS1同樣得到了反向旋轉渦環結構(圖3(c)).進而在Rp=7，L=165 mm的條件下進行計算(TS1-2)，并將獲得的計算結果與文獻[17]的計算結果(TS3)進行對比.當t=560 μs時，流場中當地流體沿x軸橫向穿過渦環中心線上的徑向速度(v)的分布規律如圖4所示，其中TS3-m為TS3的實測值.由圖4可見，TS1-2與TS3的實驗及計算結果無論是趨勢還是拐點數據均較為吻合，從而驗證了所提計算方法的準確性.

圖2 不同網格條件下Γ隨t的變化規律

圖3 TS1與TS2的渦環截面密度梯度云圖對比結果

為了研究不同可壓縮性條件下產生的渦環，設計5個具有不同Rp值的算例，每個算例對應的Mashock(由激波到達出口時根據正激波關系計算獲得)取值如表2所示.

圖4 徑向速度分布規律對比

表2 5個算例的激波馬赫數

1.2 流函數確定渦邊界法

因此，在渦環的隨體坐標系下，流函數ψ0為

(1)

式中：u0為渦環傳播速度，由各時刻渦量最大值點的所在位置確定.為了消除取渦量最大值所在位置時的測量誤差，對時間及位置曲線進行高斯擬合[21]，得到u0的平均值.

ψ0的計算步驟如下：

(1) 定義在對稱軸(r=0)上，ψ0=0；以ψ0=0的點為起點，根據速度場(u,v)進行積分，從而求得整個計算平面的ψ0值.

(2) 獲得整個區域ψ0值的曲線分布，可知ψ0=0(不包括對稱軸r=0)形成一個封閉曲線，該曲線即為渦環邊界.根據算例5計算得到的流函數分布如圖5所示.從圖中可以清晰地看到渦環邊界ψ0=0.鑒于渦環的對稱性，在ψ0=0的上半部分及r=0所包圍的區域計算環量.

圖5 流函數分布示意圖

1.3 渦核內渦量分布的無量綱化

渦核內的渦量分布是一個對理論分析(數學角度)及渦環的特征描述(物理角度)非常重要的參數.目前，已有多個用于描述渦核內渦量分布的模型被提出.Saffman[22]的研究表明，在真實流體中不可壓縮渦環的渦核內渦量分布滿足高斯分布.

對渦量分布和位置的無量綱計算公式如下：

ωn=ω/ω0

(2)

yn=(y-y0)/rc

(3)

式中：ω為渦量；ω0為峰值渦量；y為ω對應的徑向位置；y0為ω0對應的徑向位置；rc為渦核半徑.

1.4 可壓縮渦環傳播速度計算

可壓縮渦環由于自誘導速度而向下游傳播.Moore[23]提出了該傳播速度的理論計算公式：

(4)

式中：Γ可通過對由流函數方法確定的邊界區域的渦量進行積分得到；rr為渦環半徑；Mavortex為渦馬赫數，可用于表征渦環的可壓縮性，

(5)

c∞為無窮遠處的聲速.

2 結果與討論

2.1 可壓縮渦環的形成及其演化過程

激波離開激波管出口后向下游運動的同時發生衍射，導致激波管出口壓力減小，產生向激波管上游傳播的膨脹波，使得激波管內流體速度增大；被加速的流體在激波管出口形成剪切層后卷起形成渦環，渦環不斷生長并向下游傳播；渦環完成生長后，從激波管流出的流體會形成尾跡射流.尾跡射流的長度及持續時間取決于激波管驅動段的長度和壓力.

圖6 當Rp=3時 的演化過程云圖

圖7 當Rp=6時 的演化過程云圖

2.2 渦環可壓縮性的定量表征

可壓縮渦環根據結構特征可分為不含嵌入激波的渦環及含有嵌入激波的渦環.Kontis等[15]提出根據Mashock是否大于1.43來定量判斷是否存在嵌入激波，但Mashock只能作為判斷渦環是否會產生的一種參考值，對于表征渦環可壓縮性的相對強弱并不具有物理意義.Baird[12]指出把含有嵌入激波的渦環稱為超聲速渦環，并根據由渦環對稱軸上最大軸向速度計算得到的局部馬赫數(Malocal)定量區分亞聲速和超聲速渦環.此外，Moore[23]通過理論推導提出Mavortex也可用于表征渦環的可壓縮性.

傳統的空氣動力學把Ma<0.8稱為亞聲速，Ma=0.8～1.2稱為跨聲速，Ma>1.2則稱為超聲速.而在可壓縮渦環的研究中，并沒有系統地研究過渦環自身定量的可壓縮性.在不同Rp條件下，3種表征可壓縮性的馬赫數(Mashock，Malocal，Mavortex)的取值變化趨勢如圖9所示.由圖9可見：①Mashock與Rp近似成正比，且都大于1，表明都為可壓縮渦環.② 當Rp≥4時，Malocal均大于1.根據文獻[12]的定義，Rp≥4產生的皆為超聲速渦環，但當Rp=4時Malocal只是略大于1，結合圖8可知，此時產生的渦環中并不存在嵌入激波，說明產生的不是傳統意義上的超聲速渦環，而更趨向于跨聲速渦環.③ 當Rp=4時，Mavortex略小于1，表明此時的渦環處于跨聲速渦環的范圍；當Rp=5時，結合圖8可知該條件下的渦環存在嵌入激波，Malocal也遠大于1，表明此時產生的是超聲速渦環，但其Mavortex的值反而小于1.導致這種現象的原因可能是渦核受到了激波的影響，軸向被拉長，徑向被壓縮，使得渦核的半徑產生變化，導致Mavortex的數值受到了影響，但這種影響在Rp=6之后就不再顯著.

通過以上分析，為了使渦環可壓縮性的定量表征能夠與傳統的流動可壓縮性相對應，將渦環也區分為亞聲速，跨聲速和超聲速3種特征，通過Malocal和Mavortex可以很明顯地區分具有弱壓縮性的亞聲速特征渦環和具有強壓縮性的超聲速特征渦環.在跨聲速區域，由于從是否存在嵌入激波這一現象上無法區分可壓縮渦環的特征，只能通過Malocal和Mavortex的值作近似判斷.由于Malocal的計算較為簡便，故可通過其對渦環壓縮性的強弱進行預判；而Mavortex是基于整個渦環的結構參數計算獲得的，因此更能夠代表渦環自身的可壓縮性，同時也更具有實用意義.

圖8 不同Rp下主渦環ω及 的分布云圖

圖9 3種表征渦環壓縮性的Ma隨Rp的變化

3 可壓縮性對渦環物理特征的影響

本節分析Mavortex與可壓縮渦環的相關結構參數.描述渦環的重要參數包括ω，rc，rr及u0，下面將對這些參數進行定量分析.

圖10 不同Mavortex條件下ωn分布及其相應的高斯分布對比

3.1 渦核內渦量分布

渦核內無量綱渦量沿yn的分布(ωn)規律以及在對應條件下的高斯分布情況如圖10所示.其中，yn=0為渦核中心.由圖10可見：隨著Mavortex取值的增加，渦量的聚集程度也有所增加；可壓縮渦環中渦核內渦量的分布不滿足高斯分布，并且隨著Mavortex的增加，渦核內渦量更加集中，也愈加偏離高斯分布.

3.2 渦環半徑和渦核半徑

Norbury[24]根據渦環的無量綱半徑對渦環進行分類.可壓縮渦環直徑的定義如圖11所示.渦環直徑取渦量絕對值最大的兩點之間的距離，渦環半徑則取為該距離的1/2；渦核直徑由穿越渦核中心線(平行于y軸)上的軸向速度分布中最大值和最小值的所在位置確定(見圖12)， 渦核半徑則取為渦核直徑的1/2.

圖11 通過渦環截面渦量分布定義的渦環直徑

圖12 通過渦核中心u的分布定義的渦核直徑

渦環半徑和渦核半徑隨Mavortex的變化趨勢如圖13所示.由圖13可見：隨著Mavortex的增加，渦環半徑呈逐漸增加的趨勢；渦核半徑在Mavortex較小時呈逐漸增加的趨勢，而在Mavortex>1.2時，即出現相對較強的嵌入激波后，由于激波的拉伸作用，使得渦核徑向半徑略微變小；在Mavortex接近1的區域，Mavortex對兩個半徑的影響突然加劇，表明可壓縮性開始對渦環產生更劇烈的影響，也可以認為渦環進入超聲速特征的區域.

3.3 渦環傳播速度

圖14 u0及隨Mavortex的變化情況

4 結論

本文結合有限體積法求解3維Navier-Stockes方程，對激波管產生的可壓縮渦環進行數值模擬,計算結果與文獻的相關實驗結果及計算結果相符合；分析渦環的演化過程及其相關結構；給出了定量表征渦環可壓縮性的參數，并定量分析可壓縮性對渦環參數的重要影響；揭示渦環在可壓縮性作用下自身物理特征的變化情況.研究主要得到以下結論：

(1) 可壓縮渦環的可壓縮性可以通過Malocal及Mavortex定量表征.由于Malocal的計算方法較為簡單，所以可通過其對渦環壓縮性的強弱進行預判；由于Mavortex更能代表渦環自身的可壓縮性，所以更具有實用意義.

(2) 激波管產生的可壓縮渦環能夠表現出亞聲速、跨聲速和超聲速3種特征，根據Malocal和Mavortex的取值可以很明顯地區分具有弱壓縮性的亞聲速特征渦環和具有強壓縮性的超聲速特征渦環.在跨聲速區域，由于從嵌入激波這一現象上無法區分可壓縮渦環的特征，只能通過Malocal和Mavortex的值作近似判斷；而Mashock只能表征產生渦環的激波的強度，對于可壓縮渦環的特征判定不具備實際意義.

(3) 隨著可壓縮性的增加，渦環半徑呈上升趨勢，渦核半徑由于嵌入激波的存在使得其在超聲速特征區域略微有所減小，而渦環的傳播速度逐漸增加；在3種可壓縮性特征區域計算得到的渦環傳播速度和理論公式計算結果相符合，表明該理論計算公式可同時適用于3種特征區域，從而可以更為準確地預測渦環的傳播速度；可壓縮渦環渦核內渦量的分布不符合高斯分布，隨著可壓縮性的增加，渦量更加集中.

致謝感謝上海交通大學高性能計算中心超級計算機Л提供的計算幫助.