空腔流動的拉格朗日渦動力學分析

韓帥斌，羅勇，張樹海

(1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000) (2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000)

0 引 言

空腔流動是流體力學的經(jīng)典流動問題，包含豐富的渦結(jié)構(gòu)相關(guān)的物理現(xiàn)象，例如剪切層中渦的卷起與撞擊、渦與剪切層相互作用、渦聲的產(chǎn)生與傳播等。同時空腔流動具有較強的實際應(yīng)用背景。在航空航天領(lǐng)域，剪切層內(nèi)的渦與剪切層及空腔拐角相互作用產(chǎn)生的強烈噪聲不僅影響身體健康，更會導致飛行器結(jié)構(gòu)疲勞，危及飛行安全，是目前迫切需要解決的問題。自20世紀50年代起，空腔流動尤其是空腔流致噪聲得到了大量研究[1-3]?？涨涣髦略肼暸c來流參數(shù)、空腔幾何形態(tài)等密切相關(guān)。J.E.Rossiter[4]通過系列實驗總結(jié)得到Rossiter半經(jīng)驗公式，能夠很好地預測輻射噪聲的振蕩頻率與來流馬赫數(shù)的關(guān)系。通過H.H.Heller等[5]、C.K.W.Tam等[6]的完善，公式預測范圍和精度得到進一步提高。空腔流致噪聲的動力學過程為：渦擾動在剪切層中不斷增長，渦撞擊空腔后緣產(chǎn)生聲波，聲波擾動向上游傳播，聲波擾動到達前緣誘導產(chǎn)生新的渦擾動[3]。該動力學過程的一個關(guān)鍵問題是渦量擾動和壓力脈動之間如何相互轉(zhuǎn)化。M.V.Morkovin等[7]研究了壓力信號轉(zhuǎn)化為渦量擾動的物理機制；Y.P.Tang等[8]對渦-邊緣相互作用的形式進行了區(qū)分，對渦波轉(zhuǎn)化壓力波的物理機理作出了相應(yīng)的闡釋；萬振華[9]借助渦量及擬渦能的時空演化分析了空腔中渦結(jié)構(gòu)撞擊后拐角產(chǎn)生壓力脈動的動力學過程?？涨恢械臏u結(jié)構(gòu)作為噪聲的主要來源，其動力學過程對空腔流動特性及噪聲的產(chǎn)生和傳播起著重要作用，開展空腔流動中的渦動力學分析能夠為相應(yīng)的流動控制和降噪提供理論依據(jù)。

識別流場中的渦結(jié)構(gòu)是開展渦動力學分析的關(guān)鍵步驟。目前存在諸多渦判據(jù)[10-11]可用來識別流場中的渦結(jié)構(gòu)，但仍沒有公認的渦的嚴格定義。既往的渦動力學研究大多采用歐拉框架下的分析工具例如閉合流線、渦量、Q判據(jù)、 準則等識別流場中的渦，然而這些方法通常僅使用某瞬時流場的速度或速度梯度提取渦結(jié)構(gòu)，揭示的是流動的瞬時狀態(tài)，缺失了動力學系統(tǒng)中與時間相關(guān)的信息，無法準確反映流動結(jié)構(gòu)的歷史累積效應(yīng)以及相應(yīng)的動力學特性。另外歐拉框架下的判據(jù)通常不具有客觀性，即流場結(jié)果會隨參考系改變而改變。

近年來，拉格朗日方法在復雜流動系統(tǒng)的動力學分析及渦識別方面取得了較大進展[12]。拉格朗日方法在時空相空間內(nèi)跟蹤流體粒子運動，能夠客觀地揭示非線性動力系統(tǒng)的動力學特性和流動機理?；诶窭嗜諗M序結(jié)構(gòu)(Lagrangian Coherent Structures，簡稱LCS)識別流場中的渦結(jié)構(gòu)并研究其時空演化規(guī)律，能夠揭示開式空腔中的拉格朗日渦的動力學特性。

本文采用拉格朗日擬序結(jié)構(gòu)，并結(jié)合Q判據(jù)和渦量通量分析，對開式空腔流動中渦的生成、脫落、對流以及撞擊破裂過程展開研究，分析其動力學特性，為流動控制及降噪提供理論依據(jù)。

1 數(shù)值分析方法

1.1 高精度流場計算方法

通過求解非定?？蓧嚎sNavier-Stokes方程(簡記為N-S方程)對空腔流動進行直接數(shù)值模擬，守恒形式的無量綱N-S方程為

?tρ+?i(ρvi)=0

(1)

(2)

(3)

其中，

(4)

(5)

(6)

粘性系數(shù)的Sutherland公式為

(7)

N-S方程的對流項和粘性項分別采用五階WENO格式和六階中心差分格式進行離散。時間項采用三階TVD Runge-Kutta格式求解獲得高精度非定常無量綱流場數(shù)據(jù)。

1.2 拉格朗日擬序結(jié)構(gòu)

拉格朗日擬序結(jié)構(gòu)是基于追蹤流體粒子運動而提取的流場中主導流動的吸引性、排斥性或者強剪切擬序結(jié)構(gòu)，是流動的“骨架”[12]，目前已被廣泛應(yīng)用于流動分離[13]、渦動力學[14]、生物流體力學[15]等。諸多方法可以用來提取流場中的LCS，其中有限時間李雅普諾夫指數(shù)(Finite Time Lyapunov Exponent，簡稱FTLE)是被廣泛應(yīng)用的一種方法。對于在非自治系統(tǒng)

(8)

(9)

式中：λmax為矩陣的最大特征值。

在實際計算中，F(xiàn)TLE既可以前向時間積分也可以逆向時間積分，對應(yīng)的FTLE的嵴分別捕捉流場中的前向LCS(pLCS)和逆向LCS(nLCS)，pLCS通常是流場中的排斥性結(jié)構(gòu)，nLCS則通常是流場中的吸引性結(jié)構(gòu)。H.Yangzi等[16]采用pLCS和nLCS的邊界識別拉格朗日渦邊界并追蹤渦運動，本文采用該方法對空腔中的拉格朗日渦結(jié)構(gòu)進行動力學分析。

1.3 Q判據(jù)

歐拉框架下常用來識別渦的一種方法是J.C.R.Hunt等[17]提出的Q判據(jù)。速度梯度張量可以分解為對稱的應(yīng)變率張量S和反對稱的渦張量Ω，即

(10)

Q值的定義即為

(11)

在流場中Q值為正的地方，也即渦張量對流體變形的貢獻大于應(yīng)變率張量貢獻的區(qū)域被識別為渦。Q判據(jù)具有伽利略不變性，但在旋轉(zhuǎn)或加速參考系中則不能正確進行渦識別。對于本文中的開式空腔流動，參考系不存在旋轉(zhuǎn)或加速，Q判據(jù)能夠給出渦的正確描述，因而被用來與LCS識別的拉格朗日渦對比驗證，并描述相應(yīng)的渦動力學。

2 數(shù)值計算及驗證

2.1 計算模型及設(shè)置

本文直接數(shù)值模擬的空腔流動計算域及空腔附近網(wǎng)格如圖1所示，其中空腔長深比為2∶1，各項無量綱長度參數(shù)為

L=2D,Ll=7D,Lr=20D,Ly=20D

(12)

腔外網(wǎng)格為810×300，腔內(nèi)網(wǎng)格為240×120，并在腔壁附近加密以捕捉邊界層?？涨涣鲃訒a(chǎn)生強烈的聲波，為避免聲波在邊界的反射污染流場，在計算域邊界設(shè)置海綿層[18]吸收聲波，從而獲得準確的高精度流場數(shù)據(jù)。海綿層的長度參數(shù)為

Li=3D,Lo=10D,Lt=10D

(13)

本文針對亞聲速流動開展研究，分析所用的流場物理量均為無量綱量。空腔的流動參數(shù)設(shè)置為Ma=0.8，Re=2 500，Pr=0.7，其中雷諾數(shù)基于空腔深度D。

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

(b) 腔體附近網(wǎng)格分布圖1 二維空腔流動計算域及網(wǎng)格Fig.1 Computational domain and mesh of the open cavity

2.2 計算結(jié)果及驗證

本文直接數(shù)值模擬所用代碼已在參考文獻[13，19]中得到很好驗證。針對本算例，首先進行網(wǎng)格收斂性驗證。考慮一套加密網(wǎng)格，其中腔外網(wǎng)格為1 245×470，腔內(nèi)網(wǎng)格為400×200，并比較兩套網(wǎng)格下空腔底部(0.5D，0)和空腔外部(0,2D)處的壓力隨時間變化。兩套網(wǎng)格下壓力變化基本一致，如圖2所示，可以看出：幅值和相位差別均很小，說明本計算所用網(wǎng)格密度已經(jīng)足夠。

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

(b) 空腔外x=0，y=2D處壓力隨時間變化圖2 網(wǎng)格收斂性驗證Fig.2 Validation of grid convergency

對空腔(0.995D，D)處的速度和壓力脈動采樣分析，結(jié)果如圖3所示，流動呈現(xiàn)出強烈的周期性，且無量綱時間周期為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給出的半經(jīng)驗公式[20]

(14)

可見，圖4得到的主頻與式(14)中的第2模態(tài)即St2=0.686吻合較好，而且與K.Krishnamurty[21]的實驗結(jié)果頻率St=0.656吻合良好。Krishnamurty的實驗紋影與直接數(shù)值模擬獲得的數(shù)值紋影對比如圖5所示。

(a) Krishnamurty實驗紋影

(b) 數(shù)值紋影圖5 實驗紋影圖與DNS計算所得數(shù)值紋影對比Fig.5 Comparison between experimental Schlieren and Numerical Schlieren from DNS

從圖5可以看出：兩者在定性上吻合良好，波系結(jié)構(gòu)較為一致。

3 渦動力學分析

3.1 流場基本結(jié)構(gòu)

一個周期內(nèi)的流線結(jié)構(gòu)如圖6所示，空腔內(nèi)部主要有四個回流區(qū)，它們構(gòu)成了流場的主要拓撲結(jié)構(gòu)。在t=0.2T時刻，空腔尾緣點處渦量等值面撕裂，在t=0.8T時刻在前緣開始形成閉合流線結(jié)構(gòu)，這些特征在歐拉框架下通常被看作渦破裂和渦生成。然而流線和渦量都不具有客觀性，所揭示的流場結(jié)果會隨參考系的變化而變化。另外瞬時流態(tài)不能揭示渦結(jié)構(gòu)的動態(tài)演化過程，其動力學特性也無法揭示。因此接下來本文采用拉格朗日方法對流場渦結(jié)構(gòu)進行動力學特性分析。

(a) t=0

(b) t=0.2T

(f) t=T圖6 空腔流動的流線及渦量分布Fig.6 Streamlines and vorticity distribution in the cavity

3.2 拉格朗日渦動力學分析

對空腔流動的流場進行FTLE的計算，積分時長是影響計算結(jié)果的重要參數(shù)，通常隨著積分時長的增加，F(xiàn)TLE的嵴所提取的LCS結(jié)構(gòu)更加銳利。因此在本文計算中，積分時長選取為一個時間周期T=3.75，使得FTLE的嵴能夠準確捕捉到流場中的LCS。一個周期內(nèi)不同時刻的pLCS(紅色)和nLCS(藍色)如圖7所示，pLCS為排斥性結(jié)構(gòu)，附近的流體粒子被其排斥從而流向動力學特性不同的區(qū)域；nLCS則為吸引性結(jié)構(gòu)，在時空演化過程中吸引周圍來自不同動力學特性區(qū)域的流體粒子向其聚集靠攏，因此pLCS和nLCS是不同動力學特性區(qū)域的邊界。渦結(jié)構(gòu)作為流場中的擬序結(jié)構(gòu)，在渦邊界內(nèi)外，流動的特性存在差異，因而渦邊界可用LCS進行識別。將pLCS和nLCS包絡(luò)的閉合區(qū)域作為渦結(jié)構(gòu)，兩者交叉處的鞍點作為追蹤渦結(jié)構(gòu)運動的特征點來分析拉格朗日渦的動力學特性。

pLCS和nLCS包絡(luò)區(qū)域形成的渦結(jié)構(gòu)主要有三個：空腔前緣點剪切層卷起后形成的渦，空腔后緣點處即將發(fā)生撞擊和撕裂的渦以及空腔右側(cè)的主渦。這與流線所揭示的四個主要回流區(qū)有所不同，位于空腔右下拐角處的回流區(qū)Ⅱ和位于空腔左側(cè)的回流區(qū)Ⅲ未被LCS揭示出來。由于LCS揭示的是擬序結(jié)構(gòu)，能在一定生命周期內(nèi)保持動力學及結(jié)構(gòu)特征穩(wěn)定，而回流區(qū)Ⅲ即使從流線觀察不能保持穩(wěn)定的形態(tài)且無清晰邊界，因此未被LCS識別出來；回流區(qū)Ⅱ位于空腔右側(cè)底部拐角處，速度較低，其靠近壁面的邊界由于無滑移邊界速度幾乎為零，吸引性或排斥性較弱，因而只有與主渦相交的邊界被LCS揭示出來，而完整的渦結(jié)構(gòu)未被識別出來。流場中Q值大于零的區(qū)域如圖7中灰黑色區(qū)域所示。Q值揭示的渦區(qū)域也主要有三個，即空腔右側(cè)的主渦和剪切層中卷起和向下游發(fā)展的兩個渦，與LCS包絡(luò)所揭示的渦結(jié)構(gòu)覆蓋的區(qū)域基本一致?？涨挥蚁鹿战翘幍慕菧u盡管被Q值所揭示，但強度較弱，且形態(tài)和邊界不夠清晰穩(wěn)定。流線揭示的回流區(qū)Ⅲ同樣未被Q所識別如圖6所示，表明空腔左側(cè)盡管存在瞬時封閉流線，但不能保持動力學穩(wěn)定性及清晰邊界，該回流區(qū)不能識別為渦結(jié)構(gòu)。綜上，LCS與Q值都清晰地揭示了空腔流動的流場中存在的主要渦結(jié)構(gòu)，而空腔右側(cè)的角渦則由于吸引性或排斥性較弱，未被LCS識別出來。

空腔右側(cè)的主渦在整個流動周期中形態(tài)基本維持不變，繞渦心順時針旋轉(zhuǎn)，周期與流場周期一致。在運動過程中，空腔后緣點處剪切層中的渦撕裂后，一部分流體被nLCS及pLCS包絡(luò)并隨主渦運動進入到主渦中；同時主渦靠近剪切層一側(cè)流體被pLCS吸引靠近剪切層進入主流，隨主流對流至下游。

剪切層中的渦結(jié)構(gòu)則經(jīng)歷了復雜的生成、增大、脫落、對流、破裂等過程，其完整周期為流場中壓力和速度脈動周期的兩倍，而且剪切層內(nèi)兩個渦結(jié)構(gòu)的完整運動動力學過程和特性是相同的。在圖7中，t=0.4T時刻如圖7(c)所示，剪切層中的渦結(jié)構(gòu)在空腔前緣點開始生成并逐漸增大直至t=T時刻如圖7(c)～圖7(f)所示，從前緣點脫落進入到剪切層中，并隨主流一起向下游對流。在對流過程中如圖7(a)～圖7(e)所示，渦結(jié)構(gòu)發(fā)生形變扭曲，并在t=T+0.8T時刻如圖7(e)所示，開始撞擊空腔后緣點，整個撞擊過程中如圖7(f)、圖7(a)、圖7(b)、圖7(c)所示，渦結(jié)構(gòu)不斷被空腔后緣點撕裂，最終在t=2T+0.4T如圖7(c)所示，被撕裂為兩部分，一部分隨主流在空腔外向下游流動，另一部分則進入主渦部分隨主渦在空腔內(nèi)部流動。

(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

前緣點和后緣點附近區(qū)域的渦通量為

(15)

式中：ω為渦量。

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

圖8 渦通量時間演化Fig.8 Time evolution of vorticity flux

4 結(jié) 論

(1) 吸引性和排斥性拉格朗日擬序結(jié)構(gòu)邊界的包絡(luò)作為空腔流動中的渦邊界，可有效識別出空腔中形態(tài)基本維持不變的主渦和剪切層內(nèi)發(fā)生復雜動力學過程的渦?？涨挥覀?cè)的角渦則由于吸引性和排斥性較弱，未被識別出來。

(2)Q判據(jù)以及渦量通量的時間演化分析表明，在流動周期內(nèi)，拉格朗日擬序結(jié)構(gòu)在空腔前緣點的卷起與脫離伴隨著前緣渦的生成與脫落，并向下游發(fā)展；拉格朗日擬序結(jié)構(gòu)在后緣點的撞擊破裂伴隨著Q與渦的破裂。

(3) 剪切層中的渦結(jié)構(gòu)經(jīng)歷了復雜的生成、增大、脫落、對流、破裂等過程，空腔內(nèi)部的主渦則在整個流動周期內(nèi)都維持相對穩(wěn)定形態(tài)。