基于POD方法的二維方柱低雷諾數(shù)繞流流場(chǎng)分析研究

王掩剛, 陳俊旭, 先松川

(西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院, 陜西 西安 710072)

在計(jì)算流體力學(xué)的應(yīng)用中,通常求解數(shù)學(xué)模型方程組的計(jì)算規(guī)模龐大、維數(shù)較高,對(duì)計(jì)算機(jī)的容量和速度提出了極高的要求。因此如何在保證足夠精度的條件下對(duì)高維或無(wú)窮維流體動(dòng)力系統(tǒng)進(jìn)行降維建模成為了計(jì)算流體力學(xué)領(lǐng)域研究的熱點(diǎn)問(wèn)題,國(guó)內(nèi)外研究學(xué)者為此進(jìn)行了廣泛研究。

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

綜合國(guó)內(nèi)在POD方法的相關(guān)研究來(lái)看,在非定常流場(chǎng)方面的應(yīng)用研究相對(duì)較少,因此本文基于POD方法對(duì)二維方柱尾緣特征區(qū)域的非定常流場(chǎng)進(jìn)行了研究,對(duì)求得的POD模態(tài)以及時(shí)間系數(shù)進(jìn)行了初步分析,并用POD模態(tài)對(duì)原始流場(chǎng)進(jìn)行了重構(gòu),證實(shí)了該方法的可行性,為后續(xù)POD方法在流場(chǎng)方面的研究應(yīng)用提供了一定的參考。

1 計(jì)算模型及POD模型的構(gòu)建

1.1 計(jì)算模型

計(jì)算采用二維方柱為研究對(duì)象,方柱的邊長(zhǎng)D

=1 m,計(jì)算域?yàn)?1D×30D,上游邊界距方柱中心4.5D,下游邊界距方柱中心25.5D,上、下側(cè)邊界分別距方柱中心5.5D。在大量網(wǎng)格校驗(yàn)的基礎(chǔ)上,計(jì)算網(wǎng)格采用H型網(wǎng)格,網(wǎng)格密度選取370×280,并在方柱周圍以及下游尾跡區(qū)域加密,以準(zhǔn)確捕捉其流動(dòng)細(xì)節(jié)。

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

1.2 POD簡(jiǎn)介及模型的構(gòu)建

本征正交分解(POD)方法由LumIey[6]于1967年首先引入到湍流相干結(jié)構(gòu)研究中。該方法的基本思想是把原來(lái)在時(shí)間域及空間域上連續(xù)物理量的場(chǎng),寫成一個(gè)只與時(shí)間相關(guān)的函數(shù)和一個(gè)只與空間相關(guān)的函數(shù)的展開(kāi)式序列,且它們?cè)诰揭饬x上是最優(yōu)的,在展開(kāi)式中只需要少量的項(xiàng)數(shù)就可較準(zhǔn)確地描述該物理過(guò)程,從而可以提供具有足夠高精度而自由度又較少的低維模型,以降低計(jì)算維數(shù)、節(jié)省計(jì)算時(shí)間和內(nèi)存。本文將以上述方柱繞流的非定常計(jì)算模型為基礎(chǔ),運(yùn)用POD方法來(lái)構(gòu)建模型并進(jìn)行研究分析。

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

圖1 POD分析的特征區(qū)域

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

(1)

(2)

根據(jù)POD方法的基本思想,目標(biāo)是將脈動(dòng)量V(x,y,ti)分解為空間模態(tài)φi(x,y)與時(shí)間系數(shù)ai(t)的函數(shù)即:

(3)

實(shí)際上,求解空間模態(tài)φi(x)等價(jià)于求解以下最大值問(wèn)題:

(4)

且滿足

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

(5)

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

利用變分法,上述最大值問(wèn)題可轉(zhuǎn)化為以下特征值問(wèn)題:

(6)

式中:

i=1,2,…N為Vi的自相關(guān)函數(shù),也稱POD的核;可以利用原有函數(shù)空間快照脈動(dòng)量的線性組合來(lái)表示空間模態(tài),即

(7)

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

(8)

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

(9)

基于上面所求的M個(gè)POD模態(tài),原速度場(chǎng)的降維模式可以表示為

(10)

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

2 數(shù)值模擬結(jié)果及POD分析

2.1 數(shù)值模擬結(jié)果

圖2和表1是通過(guò)上述計(jì)算模型得到的方柱繞流非定常數(shù)值模擬的結(jié)果。

圖2 升阻力系數(shù)隨無(wú)量綱化時(shí)間的變化曲線

表1 方柱的計(jì)算結(jié)果與參考文獻(xiàn)[8-9]的對(duì)比

2.2 POD結(jié)果分析

圖3 各階POD模態(tài)的能量分布

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

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

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

圖4 前4階時(shí)間系數(shù)隨時(shí)間t′的變化

為了清楚地看到不同模對(duì)時(shí)間系數(shù)的相互關(guān)系,圖5給出了POD分析中前3個(gè)模對(duì)相應(yīng)的時(shí)間系數(shù)ai在60個(gè)時(shí)刻點(diǎn)相關(guān)的利薩如圖形。不同的時(shí)間系數(shù)合成軌跡為封閉的圖形,且頻率比滿足如下關(guān)系式:

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

通過(guò)以上分析,得到了前3個(gè)模對(duì)中的6個(gè)時(shí)間系數(shù)之間的頻率及其相位關(guān)系。對(duì)于一個(gè)復(fù)雜的旋渦流場(chǎng)來(lái)說(shuō),各階旋渦波動(dòng)共存,并且相互干涉。通過(guò)獲取各階波動(dòng)的關(guān)聯(lián)關(guān)系,如果能夠?qū)ο鄬?duì)較低能量的波動(dòng)實(shí)施干擾,控制流動(dòng)細(xì)微結(jié)構(gòu),從而實(shí)現(xiàn)旋渦運(yùn)動(dòng)的宏觀改善,這對(duì)于復(fù)雜流動(dòng)的主動(dòng)控制技術(shù)有一定的參考價(jià)值。

圖5 時(shí)間系數(shù)的利薩如圖形

2.3 基于POD結(jié)果的流場(chǎng)重構(gòu)

基于上述分析,對(duì)于本文所研究的低雷諾數(shù)條件下的方柱繞流,前4階POD模態(tài)表征了絕大部分的流場(chǎng)能量,為了做進(jìn)一步驗(yàn)證,下文通過(guò)前4階POD模態(tài)對(duì)流場(chǎng)進(jìn)行了重構(gòu),并比較了不同階數(shù)的POD模態(tài)對(duì)流場(chǎng)的重構(gòu)效果。

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

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

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

圖8 原始流場(chǎng)與不同階數(shù)POD重構(gòu)結(jié)果對(duì)比

3 結(jié) 論

本文的研究結(jié)果表明:

1) 應(yīng)用本文的數(shù)值分析手段,捕捉到了方柱繞流非定常流動(dòng)特征,且與公開(kāi)文獻(xiàn)結(jié)果吻合較好。

2) 對(duì)于本文的研究對(duì)象,前4階POD模態(tài)所占波動(dòng)總能量為99.4%,可以用于較準(zhǔn)確描述其非定常流場(chǎng)特征。

3) 各階POD模態(tài)的時(shí)間系數(shù)之間都有明確的頻率與相位關(guān)系,其中能量最大的第一對(duì)POD模態(tài)對(duì)方柱尾跡流動(dòng)的波動(dòng)頻率和幅值起決定作用,這對(duì)于流動(dòng)分析控制有一定的參考價(jià)值。

4) 在所研究的雷諾數(shù)下,前4階POD模態(tài)就可以很好地重構(gòu)流場(chǎng),這對(duì)以后低階模型的建立有一定的指導(dǎo)意義。

參考文獻(xiàn)：

[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] 張震宇. 風(fēng)力機(jī)翼型動(dòng)態(tài)失速的 POD 模型降階方法[J]. 南京航空航天大學(xué)學(xué)報(bào), 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] 倪振華, 江棹榮, 謝壯寧. 本征正交分解技術(shù)及其在預(yù)測(cè)屋蓋風(fēng)壓場(chǎng)中的應(yīng)用 [J]. 振動(dòng)工程學(xué)報(bào), 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òng)的二維數(shù)值分析[D]. 天津大學(xué), 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