韓帥斌,羅勇,張樹海
(1.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室,綿陽 621000) (2.中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000)
空腔流動是流體力學的經典流動問題,包含豐富的渦結構相關的物理現象,例如剪切層中渦的卷起與撞擊、渦與剪切層相互作用、渦聲的產生與傳播等。同時空腔流動具有較強的實際應用背景。在航空航天領域,剪切層內的渦與剪切層及空腔拐角相互作用產生的強烈噪聲不僅影響身體健康,更會導致飛行器結構疲勞,危及飛行安全,是目前迫切需要解決的問題。自20世紀50年代起,空腔流動尤其是空腔流致噪聲得到了大量研究[1-3]。空腔流致噪聲與來流參數、空腔幾何形態等密切相關。J.E.Rossiter[4]通過系列實驗總結得到Rossiter半經驗公式,能夠很好地預測輻射噪聲的振蕩頻率與來流馬赫數的關系。通過H.H.Heller等[5]、C.K.W.Tam等[6]的完善,公式預測范圍和精度得到進一步提高。空腔流致噪聲的動力學過程為:渦擾動在剪切層中不斷增長,渦撞擊空腔后緣產生聲波,聲波擾動向上游傳播,聲波擾動到達前緣誘導產生新的渦擾動[3]。該動力學過程的一個關鍵問題是渦量擾動和壓力脈動之間如何相互轉化。M.V.Morkovin等[7]研究了壓力信號轉化為渦量擾動的物理機制;Y.P.Tang等[8]對渦-邊緣相互作用的形式進行了區分,對渦波轉化壓力波的物理機理作出了相應的闡釋;萬振華[9]借助渦量及擬渦能的時空演化分析了空腔中渦結構撞擊后拐角產生壓力脈動的動力學過程。空腔中的渦結構作為噪聲的主要來源,其動力學過程對空腔流動特性及噪聲的產生和傳播起著重要作用,開展空腔流動中的渦動力學分析能夠為相應的流動控制和降噪提供理論依據。
識別流場中的渦結構是開展渦動力學分析的關鍵步驟。目前存在諸多渦判據[10-11]可用來識別流場中的渦結構,但仍沒有公認的渦的嚴格定義。既往的渦動力學研究大多采用歐拉框架下的分析工具例如閉合流線、渦量、Q判據、 準則等識別流場中的渦,然而這些方法通常僅使用某瞬時流場的速度或速度梯度提取渦結構,揭示的是流動的瞬時狀態,缺失了動力學系統中與時間相關的信息,無法準確反映流動結構的歷史累積效應以及相應的動力學特性。另外歐拉框架下的判據通常不具有客觀性,即流場結果會隨參考系改變而改變。
近年來,拉格朗日方法在復雜流動系統的動力學分析及渦識別方面取得了較大進展[12]。拉格朗日方法在時空相空間內跟蹤流體粒子運動,能夠客觀地揭示非線性動力系統的動力學特性和流動機理。基于拉格朗日擬序結構(Lagrangian Coherent Structures,簡稱LCS)識別流場中的渦結構并研究其時空演化規律,能夠揭示開式空腔中的拉格朗日渦的動力學特性。
本文采用拉格朗日擬序結構,并結合Q判據和渦量通量分析,對開式空腔流動中渦的生成、脫落、對流以及撞擊破裂過程展開研究,分析其動力學特性,為流動控制及降噪提供理論依據。
通過求解非定常可壓縮Navier-Stokes方程(簡記為N-S方程)對空腔流動進行直接數值模擬,守恒形式的無量綱N-S方程為
?tρ+?i(ρvi)=0
(1)
(2)
(3)
其中,
(4)
(5)
(6)
粘性系數的Sutherland公式為
(7)
N-S方程的對流項和粘性項分別采用五階WENO格式和六階中心差分格式進行離散。時間項采用三階TVD Runge-Kutta格式求解獲得高精度非定常無量綱流場數據。
拉格朗日擬序結構是基于追蹤流體粒子運動而提取的流場中主導流動的吸引性、排斥性或者強剪切擬序結構,是流動的“骨架”[12],目前已被廣泛應用于流動分離[13]、渦動力學[14]、生物流體力學[15]等。諸多方法可以用來提取流場中的LCS,其中有限時間李雅普諾夫指數(Finite Time Lyapunov Exponent,簡稱FTLE)是被廣泛應用的一種方法。對于在非自治系統

(8)

(9)
式中:λmax為矩陣的最大特征值。
在實際計算中,FTLE既可以前向時間積分也可以逆向時間積分,對應的FTLE的嵴分別捕捉流場中的前向LCS(pLCS)和逆向LCS(nLCS),pLCS通常是流場中的排斥性結構,nLCS則通常是流場中的吸引性結構。H.Yangzi等[16]采用pLCS和nLCS的邊界識別拉格朗日渦邊界并追蹤渦運動,本文采用該方法對空腔中的拉格朗日渦結構進行動力學分析。
歐拉框架下常用來識別渦的一種方法是J.C.R.Hunt等[17]提出的Q判據。速度梯度張量可以分解為對稱的應變率張量S和反對稱的渦張量Ω,即

(10)

Q值的定義即為
(11)
在流場中Q值為正的地方,也即渦張量對流體變形的貢獻大于應變率張量貢獻的區域被識別為渦。Q判據具有伽利略不變性,但在旋轉或加速參考系中則不能正確進行渦識別。對于本文中的開式空腔流動,參考系不存在旋轉或加速,Q判據能夠給出渦的正確描述,因而被用來與LCS識別的拉格朗日渦對比驗證,并描述相應的渦動力學。
本文直接數值模擬的空腔流動計算域及空腔附近網格如圖1所示,其中空腔長深比為2∶1,各項無量綱長度參數為
L=2D,Ll=7D,Lr=20D,Ly=20D
(12)
腔外網格為810×300,腔內網格為240×120,并在腔壁附近加密以捕捉邊界層。空腔流動會產生強烈的聲波,為避免聲波在邊界的反射污染流場,在計算域邊界設置海綿層[18]吸收聲波,從而獲得準確的高精度流場數據。海綿層的長度參數為
Li=3D,Lo=10D,Lt=10D
(13)
本文針對亞聲速流動開展研究,分析所用的流場物理量均為無量綱量。空腔的流動參數設置為Ma=0.8,Re=2 500,Pr=0.7,其中雷諾數基于空腔深度D。

(a) 二維空腔流動計算區域示意圖

(b) 腔體附近網格分布圖1 二維空腔流動計算域及網格Fig.1 Computational domain and mesh of the open cavity
本文直接數值模擬所用代碼已在參考文獻[13,19]中得到很好驗證。針對本算例,首先進行網格收斂性驗證。考慮一套加密網格,其中腔外網格為1 245×470,腔內網格為400×200,并比較兩套網格下空腔底部(0.5D,0)和空腔外部(0,2D)處的壓力隨時間變化。兩套網格下壓力變化基本一致,如圖2所示,可以看出:幅值和相位差別均很小,說明本計算所用網格密度已經足夠。

(a) 空腔底部x=0.5D處壓力隨時間變化

(b) 空腔外x=0,y=2D處壓力隨時間變化圖2 網格收斂性驗證Fig.2 Validation of grid convergency
對空腔(0.995D,D)處的速度和壓力脈動采樣分析,結果如圖3所示,流動呈現出強烈的周期性,且無量綱時間周期為T=3.75。對包含5個流動周期的468個流場片段進行快速傅里葉變換(FFT),采樣點的頻譜如圖4所示。

圖3 速度及壓力時間演化Fig.3 Time evolution of velocity and pressure

圖4 空腔壓力及速度脈動頻譜分布Fig.4 Frequency spectrum of velocity and pressure perturbation in the open cavity
從圖4可以看出:主頻St=0.668。
Rossiter給出的半經驗公式[20]
(14)
可見,圖4得到的主頻與式(14)中的第2模態即St2=0.686吻合較好,而且與K.Krishnamurty[21]的實驗結果頻率St=0.656吻合良好。Krishnamurty的實驗紋影與直接數值模擬獲得的數值紋影對比如圖5所示。

(a) Krishnamurty實驗紋影

(b) 數值紋影圖5 實驗紋影圖與DNS計算所得數值紋影對比Fig.5 Comparison between experimental Schlieren and Numerical Schlieren from DNS
從圖5可以看出:兩者在定性上吻合良好,波系結構較為一致。
一個周期內的流線結構如圖6所示,空腔內部主要有四個回流區,它們構成了流場的主要拓撲結構。在t=0.2T時刻,空腔尾緣點處渦量等值面撕裂,在t=0.8T時刻在前緣開始形成閉合流線結構,這些特征在歐拉框架下通常被看作渦破裂和渦生成。然而流線和渦量都不具有客觀性,所揭示的流場結果會隨參考系的變化而變化。另外瞬時流態不能揭示渦結構的動態演化過程,其動力學特性也無法揭示。因此接下來本文采用拉格朗日方法對流場渦結構進行動力學特性分析。

(a) t=0

(b) t=0.2T

(f) t=T圖6 空腔流動的流線及渦量分布Fig.6 Streamlines and vorticity distribution in the cavity
對空腔流動的流場進行FTLE的計算,積分時長是影響計算結果的重要參數,通常隨著積分時長的增加,FTLE的嵴所提取的LCS結構更加銳利。因此在本文計算中,積分時長選取為一個時間周期T=3.75,使得FTLE的嵴能夠準確捕捉到流場中的LCS。一個周期內不同時刻的pLCS(紅色)和nLCS(藍色)如圖7所示,pLCS為排斥性結構,附近的流體粒子被其排斥從而流向動力學特性不同的區域;nLCS則為吸引性結構,在時空演化過程中吸引周圍來自不同動力學特性區域的流體粒子向其聚集靠攏,因此pLCS和nLCS是不同動力學特性區域的邊界。渦結構作為流場中的擬序結構,在渦邊界內外,流動的特性存在差異,因而渦邊界可用LCS進行識別。將pLCS和nLCS包絡的閉合區域作為渦結構,兩者交叉處的鞍點作為追蹤渦結構運動的特征點來分析拉格朗日渦的動力學特性。
pLCS和nLCS包絡區域形成的渦結構主要有三個:空腔前緣點剪切層卷起后形成的渦,空腔后緣點處即將發生撞擊和撕裂的渦以及空腔右側的主渦。這與流線所揭示的四個主要回流區有所不同,位于空腔右下拐角處的回流區Ⅱ和位于空腔左側的回流區Ⅲ未被LCS揭示出來。由于LCS揭示的是擬序結構,能在一定生命周期內保持動力學及結構特征穩定,而回流區Ⅲ即使從流線觀察不能保持穩定的形態且無清晰邊界,因此未被LCS識別出來;回流區Ⅱ位于空腔右側底部拐角處,速度較低,其靠近壁面的邊界由于無滑移邊界速度幾乎為零,吸引性或排斥性較弱,因而只有與主渦相交的邊界被LCS揭示出來,而完整的渦結構未被識別出來。流場中Q值大于零的區域如圖7中灰黑色區域所示。Q值揭示的渦區域也主要有三個,即空腔右側的主渦和剪切層中卷起和向下游發展的兩個渦,與LCS包絡所揭示的渦結構覆蓋的區域基本一致。空腔右下拐角處的角渦盡管被Q值所揭示,但強度較弱,且形態和邊界不夠清晰穩定。流線揭示的回流區Ⅲ同樣未被Q所識別如圖6所示,表明空腔左側盡管存在瞬時封閉流線,但不能保持動力學穩定性及清晰邊界,該回流區不能識別為渦結構。綜上,LCS與Q值都清晰地揭示了空腔流動的流場中存在的主要渦結構,而空腔右側的角渦則由于吸引性或排斥性較弱,未被LCS識別出來。
空腔右側的主渦在整個流動周期中形態基本維持不變,繞渦心順時針旋轉,周期與流場周期一致。在運動過程中,空腔后緣點處剪切層中的渦撕裂后,一部分流體被nLCS及pLCS包絡并隨主渦運動進入到主渦中;同時主渦靠近剪切層一側流體被pLCS吸引靠近剪切層進入主流,隨主流對流至下游。
剪切層中的渦結構則經歷了復雜的生成、增大、脫落、對流、破裂等過程,其完整周期為流場中壓力和速度脈動周期的兩倍,而且剪切層內兩個渦結構的完整運動動力學過程和特性是相同的。在圖7中,t=0.4T時刻如圖7(c)所示,剪切層中的渦結構在空腔前緣點開始生成并逐漸增大直至t=T時刻如圖7(c)~圖7(f)所示,從前緣點脫落進入到剪切層中,并隨主流一起向下游對流。在對流過程中如圖7(a)~圖7(e)所示,渦結構發生形變扭曲,并在t=T+0.8T時刻如圖7(e)所示,開始撞擊空腔后緣點,整個撞擊過程中如圖7(f)、圖7(a)、圖7(b)、圖7(c)所示,渦結構不斷被空腔后緣點撕裂,最終在t=2T+0.4T如圖7(c)所示,被撕裂為兩部分,一部分隨主流在空腔外向下游流動,另一部分則進入主渦部分隨主渦在空腔內部流動。

(a) t=0

(b) t=0.2T

(c) t=0.4T

(d) t=0.6T

(e) t=0.8T

(f) t=T圖7 LCS及Q值分布Fig.7 LCS and Q value distribution in the cavity
前緣點和后緣點附近區域的渦通量為
(15)
式中:ω為渦量。
剪切層中渦的生成、脫落、撞擊和破裂時刻的量化分析通過式(15)的時間演化特性分析給出,分析結果與LCS分析給出的特征時刻相互對照驗證。選取的前后緣點附近區域分別如圖7中A、B所示,兩個時間周期內的渦通量時間演化如圖8所示,注意在兩個區域內渦通量都是負值。由A區域內的渦通量時間演化可知:在t=0.35T時刻,渦通量開始逐漸積累增加,意味著渦結構開始形成;在t=0.87T時刻渦通量達到極大值并開始減小,意味著此時渦結構開始從前緣點脫落。由B區域內的渦通量時間演化可知:在t=0.76T時刻,渦通量開始迅速增加并在t=1.12T時刻達到極大值,意味著從上游而來的渦結構開始對流至后緣點處并撞擊后緣點,隨后渦通量開始減小,直至t=1.43T時刻渦通量減小至極小值,此時渦結構撕裂成為兩部分,一部分隨主流對流而下,另一部分則進入空腔中的主渦內。以上基于渦通量的渦結構動力學特性分析與LCS分析結果基本一致,印證了基于LCS分析渦動力學的準確性。

圖8 渦通量時間演化Fig.8 Time evolution of vorticity flux
(1) 吸引性和排斥性拉格朗日擬序結構邊界的包絡作為空腔流動中的渦邊界,可有效識別出空腔中形態基本維持不變的主渦和剪切層內發生復雜動力學過程的渦。空腔右側的角渦則由于吸引性和排斥性較弱,未被識別出來。
(2)Q判據以及渦量通量的時間演化分析表明,在流動周期內,拉格朗日擬序結構在空腔前緣點的卷起與脫離伴隨著前緣渦的生成與脫落,并向下游發展;拉格朗日擬序結構在后緣點的撞擊破裂伴隨著Q與渦的破裂。
(3) 剪切層中的渦結構經歷了復雜的生成、增大、脫落、對流、破裂等過程,空腔內部的主渦則在整個流動周期內都維持相對穩定形態。