安博 孟欣雨 桑為民,2)
* (西北工業大學航空學院,西安 710072)
? (中國空氣動力研究與發展中心結冰與防除冰重點實驗室,四川綿陽 621000)
** (翼型、葉柵空氣動力學國家重點實驗室,西安 710072)
過渡流臨界特性作為流體力學中的經典問題長期受到學界的關注.Hof 等[1]在Science上撰文強調,過渡流臨界特性的研究對揭示流動物理本質起到了至關重要的作用,同時對流場演化模式起到了決定性的作用.他們認為,流場演化是個極為復雜的過程,需要細節化的分類探索來準確把握諸如湍流等復雜流場的物理特性.Avila 等[2]也在Science上撰文闡述,過渡流研究對認識流場的重要作用,對比經典的Landau-Ruelle-Takens 流場演化模式,他們認為混沌的空間擴散是湍流特性的決定性過程和固有本質.同時Graham[3]在Nature上也介紹了過渡流的概念,文中再次強調了過渡流臨界特性研究的重要性,同時指出過渡流的相關研究為解釋更為復雜的流動現象鋪平了道路.此外諸多學者[4-13]在Annual Review of Fluid Mechanics和Journal of Fluid Mechanics上先后強調了過渡流臨界特性研究在流體力學研究中的應用價值.雖然他們就不同視角從不同層面探討了不同的內容,但是學者們一致肯定了過渡流研究的重要性和必要性,在正確認識物理本質的同時不斷推進本學科的蓬勃發展,為解決更為復雜的流動問題打下了堅實的基礎.
作為過渡流臨界特性的核心研究內容,流動分岔點(flow bifurcations)表征的是流場不同流動狀態(定常流動、非定常周期性流動、非定常準周期性流動、湍流) 的物理特性臨界點.常見的流動分岔點,比如Hopf 流動分岔點,它的出現意味著此時的流動隨雷諾數增加已經從定常狀態演化至非定常周期性狀態,流動會隨時間產生周期的循環變化.根本原因是流場穩定性被破壞了,取而代之的是流場周期性特征.隨著雷諾數進一步增加,流動的非定常周期特性也逐漸遭到破壞,此時的流動雖然繼承了部分周期性流動的物理特性,但逐漸出現了準周期性流動的典型特征,即基于龐加萊映射的雙環形結構.此時的流動不再以單一頻率振蕩,流動表現出長周期和短周期的雙頻振蕩.而Neimark-Sacker 流動分岔點正是非定常周期性流動和非定常準周期性流動的臨界點;Period-doubling 流動分岔點出現后,標志著流場從非定常周期性流動躍變至湍流;再比如定常狀態下,區別不同拓撲結構流場解的Saddle-node流動分岔點和在研究流場拓撲結構時捕捉到的Pitchfork 流動分岔點.這樣的基礎研究,對準確認識流場的物理特性意義巨大.
作者在之前的研究工作[14-15]中針對三角腔頂蓋驅動內流,經典方腔頂蓋驅動內流以及上下邊雙驅動方腔內流開展了全面的流場穩定性分析研究.較為完整地揭示了以上流動問題的流場過渡流臨界特性.同時發現流場拓撲結構的 π 旋轉對稱性與流場穩定性密切相關,對流場的演化也有較大影響.根據之前的研究經驗發現對稱驅動條件對于腔體內流的過渡流臨界特性影響非常大.當驅動條件布置在相向的兩個邊界時,會導致流場內的鏡像和 π 旋轉對稱性,而這些流場對稱性又與流場的穩定性息息相關,對更好地理解流場物理規律有重要意義.而單邊鏡像對稱驅動內流的相對缺乏,沒有形成完整的體系化的研究.本文獨特的單邊鏡像驅動條件使得頂蓋的剪應力分布相較于常規邊界驅動內流發生了巨大變化,導致了流場具有更強的失穩性,其流動機理也更為復雜.為了進一步明晰鏡像對稱性與流場穩定性的關系及對流場演化的影響規律,針對鏡像對稱頂蓋驅動內流開展過渡流臨界特性研究.
本文的研究目的在于揭示鏡像對稱頂蓋驅動方腔內流的過渡流臨界特性,數值模擬工作以低馬赫數和低雷諾數為主.所涉及馬赫數為Ma=0.173 2,雷諾數為Re≤20 000,因此選取傳統格子Boltzmann方法[16-19]中的經典單松弛碰撞遷移模型作為本文數值模擬的計算模型,且使用了目前應用最廣泛的LBGK D2Q9 模型[20].其中平衡態分布函數的構造為

式中,ωi和ei分別對應離散時空模型中不同離散方向的權系數和離散速度,即

其中,c=Δx/Δt=1 為格子速度,Δx為網格步長,Δt為時間步長.格子Boltzmann 控制方程為

其中,fi為碰撞遷移前的分布函數,Fi為碰撞遷移前的分布函數,τ 為松弛時間. ρ 和u為流體粒子的宏觀密度和速度,即

均勻直角網格具有網格質量好,魯棒性高等特點,其網格結構天然契合格子Boltzmann 方法的碰撞遷移理論,因此在LBM 數值模擬中得到了廣泛應用.本文使用的均勻直角網格,網格分辨率為1024×1024,網格步長 Δx=1/1024,在之前的研究工作中[15]我們已經開展了相關的網格獨立性驗證,研究結果表明這個網格尺度足以保證計算結果的精確性和可靠性.尤其對于腔體內流流場在Re≤20 000 的計算狀態下,能夠保證y+<0.5,具備捕捉邊界層流動信息的能力.如圖1 所示,計算域特征長度L=1.0 是方腔的邊長.頂蓋鏡像對稱驅動條件為

本文在過渡流臨界特性分析研究中設計了兩個數值模擬信息采集點(如圖1 所示)Plm(x=0.25,y=0.5)和Prm(x=0.75,y=0.5) 用于記錄流場中局部速度隨時間的變化曲線.同時,為了研究流場拓撲結構的鏡像對稱性,設計了對稱性參數 ξ,定義為

圖1 計算域Fig.1 Computational domain

其中,ulm和vlm是點Plm(x=L/4,y=L/2)水平和垂直速度分量,而urm和vrm是對應的點Prm(x=3L/4,y=L/2)的水平和垂直速度分量.
本文數值模擬的邊界條件均為平直邊界(方腔四個邊).其中頂蓋為驅動邊界,其余三邊均為物面邊界.為此本文采用了經典的非平衡態外推格式[21].該格式的基本思想是,將邊界節點上的分布函數分為平衡態和非平衡態兩部分.其中,平衡態部分由平衡態分布函數的定義近似獲得,而非平衡態部分則用非平衡態外推法求解.
如圖2 所示點D,E和F為遠場邊界點,根據LBM 的演化(碰撞遷移)原理可知,在每次演化之前需要求解每個點的分布函數,對于E點其分布函數可看作兩部分: 平衡態分布函數和非平衡態分布函數,即

圖2 平直物面邊界Fig.2 Straight wall boundary


因此,可以近似求解E點的分布函數

若考慮邊界點的碰撞過程,則邊界點E的分布函數為

綜上,可以確定壁面邊界和驅動邊界點的分布函數,各邊界點宏觀物理量構造如下
根據流動的演化機理,隨著雷諾數的增加流動會從定常狀態演化為非定常狀態[22-23].本文選取Re=1000時的數值模擬結果作為鏡像對稱頂蓋驅動內流定常結果的代表.如圖3(a)和圖3(b)所示,分別介紹了該雷諾數下的流線圖和渦量圖.由流場拓撲結構可見,由于對稱驅動的作用,流動在此時保持了非常好的對稱性,兩個對稱主渦幾乎占據了整個計算域,主渦的周邊存在著兩對對稱次級渦.從主渦和次級渦的演化和分布規律可以觀察到腔體內流的典型特征.此外,本文數值模擬結果中其余的所有渦量圖都使用了如圖3(b)所示的統一色條.
圖3(c)描繪了信息采集點Plm和Prm的水平速度分量及對稱性參數隨時間的變化曲線,其中t是無量綱流場演化步數.從三個曲線中可以看出此時的流動收斂于一個定常狀態且進一步從數值方面證明,此時的流動是嚴格對稱的(ξ ≡0).相較于之前的研究工作[14-15],沒有發現類似三角腔頂蓋驅動內流的Saddle-node 流動分岔點,說明對于定常流動,只存在一種流場解,這與經典的頂蓋驅動方腔內流研究結論一致.

圖3 Re=1000 時的定常計算結果Fig.3 Steady state atRe=1000
本文選取了三個計算狀態分別揭示周期性流動(Re=1700)、準周期性流動(Re=1735) 和湍流(Re=20 000)的流場特性.
圖4 展示了不同雷諾數信息采集點Prm的水平速度分量及通過傅里葉變換之后的速度頻譜曲線.可以看出當雷諾數增加至1700 時流動已經演化為非定常周期性流動,此時流場穩定性已被破壞,取而代之的是以頻率為f=1/T=0.343 的周期性振蕩.當雷諾數進一步增加至1735 時,流動演化為非定常準周期性流動,雖然保留了周期性流動的部分特征,但是此時的流動不在以周期性單一頻率振蕩演化.從圖中可以看到代表準周期性流動典型特征的雙頻率,其中f1=0.338 是從周期性流動繼承而來的基本頻率,而f2=0.169 是伴隨周期性流動出現的調制頻率,其他頻峰則是基本頻率和調制頻率的線性組合.當雷諾數增加至20 000 時,流場已完全被湍流取代,但是從速度頻譜圖中仍可以觀察到之前準周期性流動的雙頻率,而圍繞雙頻率頻峰的基本都是寬帶噪音.此時從速度曲線圖也可以看到流動演化的無序性和隨機性.圖5 給出了對應雷諾數的速度相圖,其中閉合的單一曲線代表了周期性流動.而準周期性流動的相圖不再是單一閉合的曲線,也不是混沌狀態的一團亂麻(湍流),而是介于二者之間的過渡狀態,其物理特性既繼承了周期性流動的部分特征同時預示了可能出現的湍流.

圖4 速度頻譜圖Fig.4 Velocity spectrum

圖5 速度相圖Fig.5 Velocity phase map
結合圖6,展示了周期性流動在一個完整周期內不同時刻的流場拓撲結構.圖6(a) 給出了周期性流動完整周期內水平速度隨時間的變化曲線以及不同時刻的選取方式.此時,流動仍舊繼承了流動定常解拓撲結構的主要特征,可以觀察到兩個主渦及其附近的次級渦依舊存在,只不過整個流場以f=0.343的頻率循環振蕩.

圖6 完整周期內不同時刻渦量圖(Re=1700)Fig.6 Vorticity snapshots at different time steps within a full period T (Re=1700)
圖7 描述了準周期性流動(Re=1735)在不同龐加萊交叉點的流場拓撲結構,給出了準周期性流動龐加萊交叉點的選取方式.圖中urm代表準周期性解信息采集點Prm處沿x方向的速度分量代表準周期性解調制頻率f2對應的長周期在一個整周期內的速度隨時間t的曲線.龐加萊交叉點的選取滿足urm=-0.003 46 和的條件.圖7(b)~圖7(d)分別展示了準周期性解在不同時刻(龐加萊交叉點)的瞬時渦量圖,雖然整體流動趨勢基本一致,但是流動細節有不同的呈現(見圖7(d)所示的對應不同時刻的流函數).這是因為此時流場已經演化為準周期性流動,盡管龐加萊交叉點的選取方式一致,但是對應的是準周期性解調制頻率f2所對應的長周期在不同時刻的瞬時渦量圖.可以想象,如果此時的流動仍為周期性流動,那么這三個時刻對應的計算結果將完全一致.這也進一步證實了流動此時已經演化為準周期性流動的事實.從流場拓撲結構來看,此時的準周期性解與周期性解的差別不是特別明顯,需要從流動細節觀察區分.因為流動剛從周期性演化至準周期性不久,其準周期性特征不是特別明顯,單從流場拓撲結構很難區分.

圖7 準周期性解(Re=1735)Fig.7 A quasi-periodic solution (Re=1735)
結合圖8 分析了湍流流動狀態的流場特性,圖8(a)揭示了湍流狀態下的能量頻譜圖,可以看到湍流慣性子區的斜率為 -5/3,跟文獻[24-26]的結論一致.伴隨能量級串現象的出現,大尺度的渦結構已破碎變成細小的渦結構,而能量也依次傳遞至小尺度的渦結構直至Kolmogorov 尺度的能量耗散現象.

圖8 湍流狀態計算結果Fig.8 Results for chaos
圖8(b)給出了某一特定時刻的流場拓撲結構,其中實線代表逆時針旋轉的渦,虛線代表順時針旋轉的渦,雖然此時流場已演化至湍流,流動已變得隨機和無序,但是部分特征仍舊明顯,如流場內大尺度的渦結構都發生了破壞,此時主導流場拓撲結構的基本都是小尺度的細碎渦結構.并且每個尺度的渦結構都有對應的振蕩頻率,如圖8(b)所示,頻峰f1所對應的頻率為主渦的振蕩頻率,相較于定常結果(見圖3),主渦的結構發生了明顯破壞,且尺度變小.頻峰f2,f3,和f4分別對應不同位置的次級渦的振蕩頻率.頻峰f5和f6對應邊角處次級渦的振蕩頻率.除了這些較大尺度的渦系結構所對應的振蕩頻率,從能量頻譜曲線中還可以觀察到其他頻峰對應的附著在主渦和次級渦周圍的細碎渦的頻率.
隨著雷諾數的增加,流動會從定常狀態演化至非定常狀態,本文針對鏡像對稱頂蓋驅動方腔內流,根據研究不同雷諾數下的擾動衰減系數(Lyapunov指數) ε=ln(urm-Uˉ),發現流場穩定性最初的破壞伴隨Hopf 流動分岔點的出現而發生,這與我們之前研究工作[14-15]中的結論一致.如圖9 所示,當雷諾數從1500 增至1692 時,擾動衰減系數的斜率不斷增大,從一個負值逐漸趨向于0.當擾動衰減系數的斜率變為0 時,說明此時流動已經演化為非定常周期性流動,意味著Hopf 流動分岔點出現在雷諾數等于1691 和1692 之間.對比經典的頂蓋驅動方腔內流[14-15,27-28](ReH=8025±25),我們發現,該流場的穩定性非常差,很容易失穩.說明頂蓋鏡像驅動對于腔體內流有很強的失穩作用.

圖9 擾動衰減系數Fig.9 Perturbation decay rate
隨著雷諾數的進一步增加,從1725 增至1735,流動由非定常周期性流動演化為非定常準周期性流動.說明Neimark-Sacker 流動分岔點出現在雷諾數等于1725 和1735 之間.圖10(a)和圖10(b)分別展示了速度頻譜圖和相圖,其中黑色和紅色曲線分別代表雷諾數為1725 和1735 的計算結果.結合圖10(a),可以觀察到周期性解的振蕩頻率為f=0.338,其他頻峰均為f的整數倍,為基本頻率的諧振頻率.而準周期性解則有兩個頻率,分別為基本頻率f1=0.333,調制頻率f2=0.169,其他頻峰則是f1和f2的線性組合,如f3=f1+f2=0.502 ,f4=f1+2f2=0.671.如圖10(b)所示,周期性解的相圖是一個閉合的單一曲線,而準周期性解的相圖則由一組曲線族構成.類似頂蓋驅動和頂底雙邊驅動方腔內流[14-15],研究發現流動非定常周期性的破壞通常都是伴隨Neimark-Sacker 流動分岔點的出現而發生.沒有發現類似三角腔頂蓋驅動和四邊驅動方腔內流時所出現的Period-doubling流動分岔點.并且類似其他腔體內流流動,流動的非定常周期性并不能長期保持,隨著雷諾數增加,伴隨準周期性流動典型雙環形結構的出現,流動很快會進一步演化為準周期性流動.

圖10 周期性和準周期性計算結果對比Fig.10 Comparison between periodic and quasi-periodic solutions
經歷了Neimark-Sacker 流動分岔點之后,當雷雷諾數增加至1750 時,流場演化逐漸表現出無序性和隨機性,說明湍流出現在雷諾數等于1735 和1750之間.圖11(a)和圖11(b)分別展示了速度頻譜圖和速度相圖,其中黑色和紅色曲線分別代表雷諾數為1735和1750 的計算結果.從速度頻譜圖(圖11(a))可以明顯觀察到準周期性解的兩個頻率,分別為基本頻率f1=0.333,調制頻率f2=0.169.當雷諾數增加至1750 時,雖然仍舊可以觀察到從準周期性流動中繼承來的兩個振蕩頻率,但是被一系列寬頻噪音所包圍,說明此時流動已然變為湍流.如圖11(b)所示,準周期性解的速度相圖是閉合的曲線族(圖11(b)局部放大圖),而湍流的速度相圖明顯是一個無序、隨機的混亂系統,更進一步證實此時流動特性主要表現為湍流.至此,流動的演化路徑已基本明晰,這與頂蓋和頂底雙邊驅動內流的研究結論保持一致,即流動先從定常演化為非定常周期性,再演化為準周期性流動,最終演化為湍流.

圖11 準周期性和湍流計算結果對比Fig.11 Comparison between quasi-periodic and chaotic solutions
由于本文的驅動條件為嚴格的頂蓋鏡像對稱,所以流動在初期表現出了嚴格的對稱性,且對稱性參數一直為0(見圖3).但是當我們觀察對稱性參數(圖12)時,可以看到當Hopf 流動分岔點出現時,對稱性參數不再是0,且隨著雷諾數增大而逐漸增加.這就說明當Hopf 流動分岔點出現時,伴隨著流場穩定性的破壞,流場鏡像對稱性也發生了破壞,這與我們之前研究[15]中觀察到的結論類似,即流場 π 旋轉對稱性與流場穩定性同時喪失.同時對比與該流動較為類似的Taylor-Couette 流動,我們發現了相同的結論,鏡像對稱性的破壞往往伴隨著流場穩定性的喪失,流動不可能不經歷對稱性破壞而直接演化為湍流[29-30].

圖12 不同雷諾數的對稱性參數Fig.12 Symmetry at different Re
雖然此時對稱參數在數量級上非常小,但是足以說明流場對稱性在逐漸喪失,這樣的微小差別在肉眼觀察流場拓撲結構時很難發現.但隨著雷諾數進一步增加,流場非對稱性就可以直觀地顯現出來(見圖6).圖13 展示了流場演化之標準周期性流動時的速度曲線和對稱性參數曲線.可以清晰地看到此時的對稱參數以0.015 左右的振幅周期性振蕩.

圖13 水平速度分量及對稱性參數(Re=1700)Fig.13 Velocity and symmetry series (Re=1700)
在本文之前的研究工作中,我們發現對于頂蓋鏡像對稱驅動方腔內流這一特定流場,Hopf 流動分岔點出現在雷諾數等于1691 和1692 之間.這個結論從圖14 中也可進一步證實,當雷諾數增加至1700 時,流動從定常狀態演化至非定常周期性流動.如圖所示,紅色符號代表由初始狀態計算得到的結果,而黑色符號意味著,首先基于初始狀態計算得到周期性結果(Re=1700),然后在此基礎上計算雷諾數小于1700 的流動.其中×代表定常結果,△,□ 和 ○ 分別代表非定常周期性流動的最小值、平均值和最大值.

圖14 流動滯后現象Fig.14 Flow hysteresis
如圖所示,基于Re=1700 的周期性結果將雷諾數分別降至1600,1500 和1400,流動仍然保持周期性特性,并非之前觀察到的定常結果(由初始態計算得到),說明存在流動滯后(flow hysteresis)現象.進一步降低雷諾數至1350,此時流動才回落至定常狀態,說明流動滯后現象發生在 1350 <Re<1700 這個區間.同時也說明此前捕捉到的Hopf 流動分岔點為亞臨界(subcritical)形式.流動滯后現象的出現意味著在 1350 <Re<1700 這個區間,對應的每個雷諾數會有兩種解的可能性,一種是定常狀態,另一種是周期性狀態.并且,根據觀察得到,此周期性流動的基本規律與之前基于初始狀態計算得到的周期性解特性基本一致.如圖15 所示,展示了Re=1600 時的周期性解在一個完整周期內不同時刻的渦量圖.

圖15 完整周期內不同時刻渦量圖(Re=1600)Fig.15 Vorticity snapshots at different time steps within a full period T (Re=1600)
針對鏡像對稱頂蓋驅動方腔內流,本文開展了流動從定常流動到湍流的數值模擬和流場穩定性分析研究,捕捉并解釋各種流動現象,從物理層面揭示該流場的流動機理,具體結論如下.
(1) 流場穩定性的破壞是以Hopf 流動分岔點的出現而開始.
(2) 相較于經典頂蓋方腔驅動內流,流場穩定性更容易喪失,Hopf 流動分岔點的臨界雷諾數為ReH=1691.5±0.5.
(3) 流場穩定性被破壞的同時,也喪失了流場鏡像對稱性.
(4) 流動喪失穩定性后會迅速從非定常周期性流動演化為非定常準周期性流動,Neimark-Sacker 流動分岔點出現在ReNS=1712.5±12.5.
(5) 當雷諾數增至ReC=1762.5±12.5,湍流出現,流動變得無序隨機.
(6) 當流動進一步演化后,隨著雷諾數增大,大尺度的渦結構發生破壞變成小尺度的細碎渦結構,同時能量從大尺度的渦結構傳遞至小尺度的渦結構直至Kolmogorov 尺度的能量耗散.
(7) 流場演化遵循經典的Ruelle-Takens 模式,從定常演化為非定常周期性流動,再到準周期性流動,最后演化為湍流.
(8) 在 1350 <Re<1700 這個區間存在流動滯后現象,對于同一個雷諾數有兩種可能的解,一種是定常流動,另一種是非定常周期性流動.并且發現Hopf流動分岔點為亞臨界型.
為了更好地分析上文介紹的三種典型流動狀態(非定常周期性流動、非定常準周期性流動和湍流)的流場拓撲結構特征,本文準備了相關視頻動畫進一步展示了流場細節.同時對應流動滯后現象,還準備了對應同一個雷諾數可能的非定常周期性解的流動特性.