基于POD方法的二維方柱低雷諾數繞流流場分析研究

王掩剛, 陳俊旭, 先松川

(西北工業大學 動力與能源學院, 陜西 西安 710072)

在計算流體力學的應用中,通常求解數學模型方程組的計算規模龐大、維數較高,對計算機的容量和速度提出了極高的要求。因此如何在保證足夠精度的條件下對高維或無窮維流體動力系統進行降維建模成為了計算流體力學領域研究的熱點問題,國內外研究學者為此進行了廣泛研究。

本征正交分解(proper orthogonal Deco-Mposition,POD)方法作為一種有效的降維手段,被廣泛運用于流體動力系統的降維建模中。在國外,Favier等[1]通過POD降階模型對圓柱繞流尾跡,以及某翼型的表面分離流動進行了分析研究,結果表明該方法可以很好地模擬所研究的流場區域;Siegel等[2]從數值模擬中提取出圓柱繞流瞬態流場在短時間內的快照(snapshots)集合,運用POD方法對該流場進行了分析研究,證實了該方法可以準確的模擬和評估瞬態流場結構;Scarano等[3]則通過PIV技術和POD方法研究了入射角對二維的方柱繞流流場的影響,結果表明該方法可以很好的展現流場的旋渦脫落過程以及入射角的影響。在國內也有學者對POD方法的進行了一定的相關研究。張震宇[4]利用POD原理設計了一種針對風力機翼型動態失速時變過程的辨識方法, 結果表明該降階模型方法能夠以明顯降低的計算量精確辨識翼型的淺失速情況;倪振華、江棹榮等[5]基于POD近似的時間與空間分解來預測未知區域的風壓時間序列,研究表明該方法能夠有效合理地反映出風壓場的時間序列。

綜合國內在POD方法的相關研究來看,在非定常流場方面的應用研究相對較少,因此本文基于POD方法對二維方柱尾緣特征區域的非定常流場進行了研究,對求得的POD模態以及時間系數進行了初步分析,并用POD模態對原始流場進行了重構,證實了該方法的可行性,為后續POD方法在流場方面的研究應用提供了一定的參考。

1 計算模型及POD模型的構建

1.1 計算模型

計算采用二維方柱為研究對象,方柱的邊長D

=1 m,計算域為11D×30D,上游邊界距方柱中心4.5D,下游邊界距方柱中心25.5D,上、下側邊界分別距方柱中心5.5D。在大量網格校驗的基礎上,計算網格采用H型網格,網格密度選取370×280,并在方柱周圍以及下游尾跡區域加密,以準確捕捉其流動細節。

計算工質為空氣,雷諾數Re=100,采用層流模型;選用壓力基求解器對各控制方程進行求解,壓力速度耦合算法選用SIMPLE算法,梯度項選用least squares cell based方法;動量項選用二階迎風格式。上游邊界指定為速度入口邊界條件u=0.001 46 m/s,v=0 m/s;下游邊界指定為壓力出口邊界條件;方柱表面指定為壁面無滑移邊界;上下邊界指定滑移邊界條件u=0.001 46 m/s,v=0 m/s。

1.2 POD簡介及模型的構建

本征正交分解(POD)方法由LumIey[6]于1967年首先引入到湍流相干結構研究中。該方法的基本思想是把原來在時間域及空間域上連續物理量的場,寫成一個只與時間相關的函數和一個只與空間相關的函數的展開式序列,且它們在均方意義上是最優的,在展開式中只需要少量的項數就可較準確地描述該物理過程,從而可以提供具有足夠高精度而自由度又較少的低維模型,以降低計算維數、節省計算時間和內存。本文將以上述方柱繞流的非定常計算模型為基礎,運用POD方法來構建模型并進行研究分析。

圖1中所示的點劃線矩形區域即為本文進行POD分析的特征區域,區域范圍為X:-9.5～-4.5;Y:-2～2。選取該區域的原因在于它是方柱繞流流動的主要特征區域，能夠反映出尾跡脫落渦的整個演化過程。

圖1 POD分析的特征區域

針對以上特征區域,在流場穩定周期性波動的某段時間內,本文選取了該區域中N=60個時刻的瞬態速度場快照(snapshots)集合U(x,y,ti)(i=1,2,…,N)用于建立降維的POD空間。定義這組樣本的平均值和脈動量分別為:

(1)

(2)

根據POD方法的基本思想,目標是將脈動量V(x,y,ti)分解為空間模態φi(x,y)與時間系數ai(t)的函數即:

(3)

實際上,求解空間模態φi(x)等價于求解以下最大值問題:

(4)

且滿足

(φ,φ)=‖φ‖2=1

(5)

式中：(·,·)和‖·‖分別表示內積空間Ω上的L2-內積和L2-范數。

利用變分法,上述最大值問題可轉化為以下特征值問題:

(6)

式中:

i=1,2,…N為Vi的自相關函數,也稱POD的核;可以利用原有函數空間快照脈動量的線性組合來表示空間模態,即

(7)

將上(7)式代入(6)式中得以下特征值問題:

(8)

Φ=span(φ1,φ2,…,φN)

(9)

基于上面所求的M個POD模態,原速度場的降維模式可以表示為

(10)

以上便得到了由M個POD基函數擴展成的降維空間,使其在降維空間里求解。用Galerkin逼近將模式方程投影到POD基函數擴張成的降維空間中,得到求解時間系數ai(t)的常微分方程組,便可求得時間系數ai(t),最終得到POD的低階模型。

2 數值模擬結果及POD分析

2.1 數值模擬結果

圖2和表1是通過上述計算模型得到的方柱繞流非定常數值模擬的結果。

圖2 升阻力系數隨無量綱化時間的變化曲線

表1 方柱的計算結果與參考文獻[8-9]的對比

2.2 POD結果分析

圖3 各階POD模態的能量分布

由圖3可以看出,所計算出的特征值均成對出現,而且一對特征值所占的能量基本相同。其中第一對POD模態所占總波動能量的比例分別為51.62%與44.08%,而剩下的POD模態分別所占總能量比例均不大于2%。第二對POD模態則與相對較低能量的高次諧波相關[7],而更高階的POD模態則是包含了流體運動中的更高次諧波。

根據之前E(k)的定義可以得出前2階POD模態占總能量的比例為95.7%,一大部分的流動動能均包含在前2階模態當中;前4階POD模態所占總能量的比例已高達99.4%。由此可以看出前4階POD模態已經完全可以抓住流場流動的主要特征,因此本文僅給出前4階POD模態的速度基函數圖像,并對前4階POD的相關結果進行了初步分析以及對非定常流場的重構。

圖4分別給出了POD分析中所占能量比例較高的前2個模對相應的時間系數(a1,a2)、(a3,a4)隨無量綱化時間的變化曲線。從圖中可以看出:時間系數基本都是呈正余弦曲線的變化趨勢,且每個模對中2個時間系數的頻率及振幅基本相同,其中時間系數a1、a2的變化周期為旋渦脫落周期的1/2;由圖5a)～圖5c)還可以得到每個模對的2個時間系數相位差為π/2,幅值隨著模對階數提高而減小,但第一個模對的時間系數幅值占主導地位。據此,可以推測:對于方柱繞流的尾跡旋渦流動,隨著旋渦的生成、發展和消亡,在微尺度上存在著多個相同頻率和幅值、但旋轉相位有差異的對旋渦對(counter-rotating vortex)。能量最大的第一階渦對決定了尾跡流動的波動頻率和幅值,低能量的高階渦對則影響著流場的細微結構。

圖4 前4階時間系數隨時間t′的變化

為了清楚地看到不同模對時間系數的相互關系,圖5給出了POD分析中前3個模對相應的時間系數ai在60個時刻點相關的利薩如圖形。不同的時間系數合成軌跡為封閉的圖形,且頻率比滿足如下關系式:

式中：fx、fy分別為2個信號的頻率,nx、ny分別為利薩如圖形的外切水平線和垂直線的切點數。從圖5d)～圖5f)可看出:a1、a4的頻率比值為1∶2;a1、a5的頻率比值為1∶3;a4、a6的頻率比值為2∶3。應用利薩如圖形在不同頻率信號疊加時其形狀與相位關系可以得出:時間系數a1與a4、a5的等效相位差均為nπ/2(n=1,3,5…),而時間系數a4與a5的等效相位差nπ/4(n=1,3,5,…)。

通過以上分析,得到了前3個模對中的6個時間系數之間的頻率及其相位關系。對于一個復雜的旋渦流場來說,各階旋渦波動共存,并且相互干涉。通過獲取各階波動的關聯關系,如果能夠對相對較低能量的波動實施干擾,控制流動細微結構,從而實現旋渦運動的宏觀改善,這對于復雜流動的主動控制技術有一定的參考價值。

圖5 時間系數的利薩如圖形

2.3 基于POD結果的流場重構

基于上述分析,對于本文所研究的低雷諾數條件下的方柱繞流,前4階POD模態表征了絕大部分的流場能量,為了做進一步驗證,下文通過前4階POD模態對流場進行了重構,并比較了不同階數的POD模態對流場的重構效果。

圖6和圖7給出了應用POD方法構造的速度分量前四階POD基函數等值線云圖。由圖中可以看到成對的相似模式和空間變化,這與方柱繞流尾部典型旋渦脫落的對流特點有關。對于流向速度u而言,前2階POD模態在橫向出現了正負交替的反對稱渦核結構,流向則是正負交替的渦核;3階和4階POD模態在橫向呈現對稱結構,流向則出現正負交替的渦核。對于橫向速度v而言,前2階POD模態顯示出橫向完全對稱,流向交替出現的渦核結構;3、4階POD的模態則在橫向表現為反對稱,流向交替出現的渦核結構。這一結果與Van Oudheusden B W等[3]、Ben Chiekh M等[7]分析結果基本一致。

圖6 流向速度u的前4階POD基函數等值線云圖

圖7 橫向速度v的前4階POD基函數等值線云圖

圖8 原始流場與不同階數POD重構結果對比

3 結 論

本文的研究結果表明:

1) 應用本文的數值分析手段,捕捉到了方柱繞流非定常流動特征,且與公開文獻結果吻合較好。

2) 對于本文的研究對象,前4階POD模態所占波動總能量為99.4%,可以用于較準確描述其非定常流場特征。

3) 各階POD模態的時間系數之間都有明確的頻率與相位關系,其中能量最大的第一對POD模態對方柱尾跡流動的波動頻率和幅值起決定作用,這對于流動分析控制有一定的參考價值。

4) 在所研究的雷諾數下,前4階POD模態就可以很好地重構流場,這對以后低階模型的建立有一定的指導意義。

參考文獻：

[1] Favier J, Cordier L, Kourta A. Accurate POD Reduced-Order Models of Separated Flows[J]. Physics of Fluids, 2007, 8(3): 259-265

[2] Siegel S, Cohen K, Seidel J, et al. Short Time Proper Orthogonal Decomposition for State Estimation of Transient Flow Fields[C]∥43rd AIAA Airspace Sciences Meeting and Exhibit, 2005: 296

[3] Van Oudheusden B W, Scarano F, Van Hinsberg N P, et al. Phase-Resolved Characterization of Vortex Shedding in the Near Wake of a Square-Section Cylinder at Incidence[J]. Experiments in Fluids, 2005, 39(1): 86-98

[4] 張震宇. 風力機翼型動態失速的 POD 模型降階方法[J]. 南京航空航天大學學報, 2011, 43(5): 577-580

Zhang Zhenyu. Reduced-Order POD Model for Dynamic Stall of Wind Turbine Airfoils [J]. Nanjing University of Aeronautics and Astronautics, 2011, 43(5): 577-580 (in Chinese)

[5] 倪振華, 江棹榮, 謝壯寧. 本征正交分解技術及其在預測屋蓋風壓場中的應用 [J]. 振動工程學報, 2007, 20(1): 1-8

Ni Zhenhua, Jiang Zhaorong, Xie Zhuangning. POD Technique and Its Application to Prediction of Wind Pressure Fields on Roof [J]. Journal of Vibration Engineering, 2007, 20(1): 1-8 (in Chinese)

[6] Lumley J L. The Structure of Inhomogeneous Turbulent Flows.∥Atmospheric Turbulence and Radio Wave Propagation[M]. Yaglom AM, Tatarski VI. Nauka, Moscow, 1967: 166-178

[7] Ben Chiekh M, Michard M, Guellouz M S, et al. POD Analysis of Momentumless Trailing Edge Wake Using Synthetic Jet Actuation[J]. Experimental Thermal and Fluid Science, 2012, 46: 89-102

[8] 余化軍. 圓柱和方柱繞流及矩形柱渦激振動的二維數值分析[D]. 天津大學, 2012

Yu Huajun. Two Dimensional Numerical Analysis of Flow over a Circular and Square Cylinder and Vortex-Induced Vibration of Rectangular Cylinder[D]. Tianjin: Tianjin University, 2012 (in Chinese)

[9] Singh A P, De A K, Carpenter V K, et al. Flow Past a Transversely Oscillating Square Cylinder in Free Stream at Low Reynolds Numbers[J]. International Journal for Numerical Methods in Fluids, 2009, 61(6): 658-682