剛性薄平板繞流特性的數(shù)值模擬

蔚杰, 武建軍

(蘭州大學(xué) 土木工程與力學(xué)學(xué)院,甘肅 蘭州 730000)

平板繞流是一類鈍體繞流現(xiàn)象,廣泛存在于自然界,如風(fēng)對植物葉片的繞流[1-2]。平板也是一種常見的工程構(gòu)件,大量應(yīng)用于船舶工程、能源動力工程、仿生設(shè)計(jì)等領(lǐng)域,如水面推進(jìn)器的旋轉(zhuǎn)葉片[3-4]、海洋平臺垂蕩板[5]、核反應(yīng)堆中的板狀燃料組件[6]、仿生飛行器的撲動翼板[7-8]。平板構(gòu)件周圍流體的流動特性直接影響著平板的穩(wěn)定性以及使役性能,甚至決定著裝置中其他結(jié)構(gòu)的設(shè)計(jì)。平板繞流系統(tǒng)還涉及復(fù)雜的轉(zhuǎn)捩過程與流動分離現(xiàn)象,尚未被完全揭示的流動機(jī)理,使得對該問題的研究具有廣泛的學(xué)術(shù)價(jià)值和工程實(shí)踐意義。與圓(橢圓)柱繞流的可變分離點(diǎn)不同,垂直平板繞流的分離點(diǎn)固定,位于平板邊緣。此外,由于剪切層不能重新附著于平板尾部,垂直平板近壁尾流的初始渦卷以及渦分離過程與其他鈍體繞流存在明顯差異,而且層流到湍流的轉(zhuǎn)捩過程發(fā)生在尾流區(qū)而非分離剪切層[9-10]。平板周圍的非定常流動包含2個主頻,除了與Krmn渦街對應(yīng)的頻率外,還具有低頻不穩(wěn)定的流動特征。平板扁平的幾何特征使得其具有較大的阻力系數(shù)[11],低頻周期內(nèi),平板在高平均阻力狀態(tài)和低平均阻力狀態(tài)之間反復(fù)切換[12]。剪切來流繞流各類鈍體時,來流與固壁邊界層相互作用,形成了多種尾渦脫落模式,影響了流場由層流發(fā)展到湍流的轉(zhuǎn)捩過程,而尾渦脫落頻率、渦流強(qiáng)度也依賴于來流剪切梯度[13-14]。此外,結(jié)構(gòu)承受多頻特征的流體載荷,甚至誘發(fā)復(fù)雜的自激振動現(xiàn)象[15-18]。剪切來流繞流鈍體的問題中包含了更為豐富的物理現(xiàn)象,其中所涉及的動力學(xué)問題值得深入的研究。另一方面,相比其他鈍體結(jié)構(gòu),平板繞流受到的關(guān)注較少,特別是平板弦厚比L/t極大的情況。

本文采用浸沒邊界-格子Boltzmann方法(IB-LBM)對垂直薄平板的繞流問題進(jìn)行了數(shù)值模擬。研究了雷諾數(shù)、來流剪切梯度對平板尾流渦街演化規(guī)律以及平板阻力特性的影響。

1 力學(xué)模型與數(shù)值方法

1.1 力學(xué)模型

如圖1所示,長度為L的剛性薄平板被固定在二維粘性不可壓流場中;假設(shè)薄平板厚度t遠(yuǎn)遠(yuǎn)小于長度。坐標(biāo)原點(diǎn)位于平板中點(diǎn),x軸與平板法向重合,z軸與平板切向重合。計(jì)算域42L×22L,平板距入口邊界11L。

圖1 計(jì)算域與邊界條件

1.2 控制方程

對于長度L遠(yuǎn)遠(yuǎn)大于厚度t的剛性平板,由于繞流結(jié)構(gòu)最小尺寸的限制,若采用全因素建模并結(jié)合貼體網(wǎng)格的計(jì)算方式,則需要在平板周圍劃分了大量的計(jì)算網(wǎng)格。浸沒邊界-格子Boltzmann方法(IB-LBM)通過引入體力項(xiàng)fv描述固體邊界與流體間的相互作用,將邊界產(chǎn)生的拉格朗日力插值并分配給周圍的 Euler節(jié)點(diǎn),并采用LBM對流場域進(jìn)行求解。如圖2所示,受文獻(xiàn)[19]的啟發(fā),以一維固壁邊界等效薄平板并將其離散為一組拉格朗日點(diǎn)集。薄平板繞流問題即可轉(zhuǎn)化為對含外力項(xiàng)LBE方程的求解。此外,可以用幾何形狀規(guī)則的Cartesian網(wǎng)格對計(jì)算域進(jìn)行離散。該方法有效地簡化了此類平板繞流問題的求解過程。

圖2 平板繞流問題浸沒邊界原理局部示意

含外力項(xiàng)的LBE方程[20]:

fα(r+eαδt,t+δt)-fα(r,t)=

(1)

(2)

本文采用D2G9模型[21]作為離散速度模型,該模型以壓力作為獨(dú)立變量,速度配置為:

(3)

平衡態(tài)分布函數(shù)為:

(4)

(5)

解得粒子分布函數(shù)后,速度和壓力可得:

(6)

(7)

式中密度ρ0為常數(shù)。

常用的浸沒邊界法[22-24]存在無法精確滿足無滑移邊界條件的固有缺陷。Wu 等[25]提出了基于隱式速度修正的浸沒邊界格子Boltzmann方法 (implicit velocity correction-based immersed boundary-lattice boltzmann method)。將式 (6)中的速度u分解為:

u=u*+δu

(8)

式中:u*為未經(jīng)修正的速度;δu為修正速度。該方法將體力fv作為未知量,并通過隱式求解速度u的方式使得拉格朗日邊界點(diǎn)滿足無滑移邊界條件:

(9)

ρ0δu=(fvδt)/2

(10)

對于圖2中,拉格朗日邊界Γ上的離散點(diǎn)XB(sk)(k=0,1,…,N),引入修正量δuB并將其設(shè)為未知量,Euler節(jié)點(diǎn)的修正速度δu可由δuB插值得到:

(11)

(12)

式中:D為Delta函數(shù);核函數(shù)δ(r)取為:

(13)

將式(12)代入式 (11)得:

(14)

拉格朗日邊界點(diǎn)速度UB由Delta函數(shù)插值得到,且需滿足無滑移邊界條件:

(15)

將式 (8)、(14)代入式 (15)可得:

(16)

作用于固壁邊界的流體力fb為:

fb=-2ρ0δuB/δt

(17)

由此可得結(jié)構(gòu)的阻力、升力為:

(18)

進(jìn)一步,阻力系數(shù)與升力系數(shù)分別被定義為[26]:

(19)

1.3 計(jì)算域離散

標(biāo)準(zhǔn)格子Boltzmann方法只允許在均勻網(wǎng)格系統(tǒng)中完成遷移-碰撞步,存在網(wǎng)格生成量大、計(jì)算精度不可調(diào)的缺點(diǎn)。為了提高計(jì)算效率,本文以非均勻網(wǎng)格對計(jì)算域進(jìn)行離散,并根據(jù)泰勒展開和最小二乘格子Boltzmann方法(Taylor series expansion-and least square-based lattice Boltzmann method, TLLBM)對粒子分布函數(shù)進(jìn)行求解,TLLBM的原理與實(shí)現(xiàn)方法可參照文獻(xiàn)[27]。

如圖3所示,直線x/L=±1,z/L=±1將計(jì)算域分為9個子區(qū)域,以映射網(wǎng)格技術(shù)生成網(wǎng)格。圖中分以數(shù)字(1)～(9)對子區(qū)域1～9進(jìn)行標(biāo)識。子區(qū)域5采用均勻網(wǎng)格,Δx0=Δz0=L/80;浸沒邊界位于(5)中。其余各子區(qū)域的最小網(wǎng)格尺寸Δxmin=Δzmin=Δx0;子區(qū)域1、子區(qū)域2、子區(qū)域4、子區(qū)域7、子區(qū)域8最大網(wǎng)格尺寸Δxmax=Δzmax=2.2Δx0,子區(qū)域3、子區(qū)域6、子區(qū)域9最大網(wǎng)格尺寸Δxmax=Δzmax=3.0Δx0。計(jì)算域網(wǎng)格數(shù):1 860×1 160。

圖3 計(jì)算域分區(qū)策略與局部網(wǎng)格示意

1.4 數(shù)值驗(yàn)證

本文采用C++結(jié)合OpenMP并行技術(shù)編寫計(jì)算程序,為了驗(yàn)證數(shù)值方法與程序的可靠性,對二維圓柱繞流進(jìn)行模擬。圓柱繞流是鈍體繞流中的經(jīng)典問題,可在大量文獻(xiàn)中獲得實(shí)驗(yàn)及數(shù)值結(jié)果。圖4給出了基于上述算法得到的定常、非定常狀態(tài)下的局部流線圖,圓圈為圓柱邊界。定常狀態(tài)下2個反向旋轉(zhuǎn)的Fopple渦對稱的附著于圓柱尾部;隨著Re增大,Fopple渦交替脫落并向下游運(yùn)動,流動呈現(xiàn)非定常特征。

圖4 圓柱繞流局部流線圖

由表1所示,本文得到的結(jié)果與文獻(xiàn)結(jié)果吻合較好,當(dāng)前數(shù)值方法和程序是有效、可靠的,且離散方案能夠保證計(jì)算精度。通過修改Lagrange點(diǎn)集的分布函數(shù)即可實(shí)現(xiàn)對薄平板的繞流問題的求解。

表1 圓柱繞流Re=20, 40, 100阻力系數(shù)、斯特羅哈數(shù)計(jì)算結(jié)果

2 雷諾數(shù)對平板繞流的影響

圖1所示的平板繞流問題,雷諾數(shù)被定義為Re=U∞L/ν。本文的模擬中,100≤Re≤1 200。隨著Re的進(jìn)一步增加,考察繞流結(jié)構(gòu)的平均積分參數(shù),如:阻力系數(shù)、斯特羅哈數(shù),參數(shù)對Re的變化并不敏感[9]。本文在近壁尾流區(qū)(x/L<20)對尾流的2次失穩(wěn)躍遷機(jī)理開展了分析與研究。

模型的邊界條件為:

入口(x=-11L):均勻來流,u=U∞;w=0。

出口(x=31L):自由出流,?(u,w,p)/?x=0。

上、下邊界(z=11L,z=-11L):對稱邊界,?(u,p)/?z=0;w=0。

2.1 雷諾數(shù)對平板尾流的影響

本文通過模擬得出,對于均勻來流條件下的平板繞流系統(tǒng),由定常轉(zhuǎn)變?yōu)榉嵌ǔ顟B(tài)的臨界雷諾數(shù)11015,漩渦由P模式向具有更大的空間尺度的P+S模式轉(zhuǎn)變。二次渦街尚未充分發(fā)展的流向范圍內(nèi)(如Re=400, 1420),可以明顯觀察到以2S漩渦模式為主的尾流結(jié)構(gòu);該區(qū)域內(nèi)頻繁的出現(xiàn)失配單渦,這些失配單渦具有不規(guī)則的運(yùn)動軌跡與空間分布特征。此外,尾流區(qū)出現(xiàn)一次躍遷、二次躍遷的流向位置具有時變特征;隨著Re的增大,出現(xiàn)躍遷的流向位置向上游移動。

圖5 瞬時渦量云圖

圖6 監(jiān)測點(diǎn)u(t*)功率譜密度函數(shù)

2.2 雷諾數(shù)對平板阻力的影響

圖7 平板阻力系數(shù)功率譜密度

圖8 瞬時局部渦量場與監(jiān)測點(diǎn)p0(2L,0)功率譜密度

3 來流剪切參數(shù)對平板繞流的影響

剪切來流條件下,平板繞流問題的邊界條件:

入口(x=-11L):速度剪切來流,線性剪切段u(z)=Uc+Gz,G為剪切梯度。

邊界x=31L,z=11L,z=-11L處的邊界條件同上。

定義剪切參數(shù)β=GL/Uc,本節(jié)中,0<β<0.166 7;此外,Re=UcL/ν。

3.1 來流剪切參數(shù)對平板尾流的影響

圖9給出了平板尾流渦脫頻率St隨剪切參數(shù)β的變化規(guī)律。Re=100時,均勻來流條件下的平板繞流系統(tǒng)處于定常狀態(tài),而在剪切來流的作用下,尾流出現(xiàn)流動分離現(xiàn)象,系統(tǒng)進(jìn)入非定常狀態(tài)。值得注意的是,剪切來流對尾渦脫落頻率具有一定的抑制作用,特別是在0.06<β<0.12范圍內(nèi)。

圖9 平板繞流系統(tǒng)St隨β變化規(guī)律

圖10 平板繞流瞬時渦量云圖

圖11 p0(2L,0)點(diǎn)流速功率譜

3.2 來流剪切參數(shù)對平板阻力的影響

圖12 平板阻力系數(shù)功率譜密度

本文采用C-C法對Cd(t*)時間序列進(jìn)行相空間重構(gòu)并由小數(shù)據(jù)量法計(jì)算最大Lyapunov指數(shù)[33-34],用以對平板Cd(t*)的混沌特性進(jìn)行識別,并據(jù)此對Cd(t*)的混沌程度進(jìn)行定量描述。如圖13,本文預(yù)測Re=200,阻力系數(shù)Cd(t*)轉(zhuǎn)變?yōu)榛煦缣卣鲿r間序列的臨界剪切參數(shù)βc≈0.16。Re=400, 600, 800對應(yīng)的臨界剪切參數(shù)分別為βc=0.14, 0.10, 0.03。Re越大,對應(yīng)的的臨界剪切參數(shù)越小。在工程實(shí)踐中,減小雷諾數(shù)或來流剪切梯度均可起到抑制平板流向載荷向混沌特征轉(zhuǎn)變的作用。

圖13 平板阻力系數(shù)的最大Lyapunov指數(shù)λ1隨來流剪切參數(shù)β的變化

4 結(jié)論

1)雷諾數(shù)Re與來流剪切梯度β對垂直平板尾流渦街的影響體現(xiàn)為:對尾流區(qū)兩次失穩(wěn)躍遷的流向位置以及對主渦街、雙排渦、二次渦街演化過程的影響。

2)在主渦街與雙排渦內(nèi),漩渦運(yùn)動軌跡是確定的,流動是周期性的。二次躍遷伴隨著漩渦失配-重組現(xiàn)象以及漩渦模式由P模式向大尺度的P+S或2S模式的轉(zhuǎn)變。二次渦街充分發(fā)展的區(qū)域,大尺度渦結(jié)構(gòu)以及失配單渦的運(yùn)動軌跡均是不規(guī)則的,引發(fā)流場的低頻、非周期性振蕩。

3)Re或β的增大均可導(dǎo)致尾流區(qū)一次躍遷與二次躍遷的流向位置向上游移動。

5)剪切來流與平板尾渦結(jié)構(gòu)相互作用,誘發(fā)了多頻特征的流向載荷。隨著來流剪切參數(shù)的增大,雙排渦結(jié)構(gòu)逐漸失穩(wěn)甚至消失,平板近尾流受大尺度低頻結(jié)構(gòu)主導(dǎo),阻力系數(shù)Cd(t*)具有低頻、混沌的特征。在工程實(shí)踐中,減小雷諾數(shù)或來流剪切梯度均可起到抑制平板流向載荷向混沌特征轉(zhuǎn)變的作用。

上述研究結(jié)果可為流體機(jī)械葉片設(shè)計(jì)、平板結(jié)構(gòu)減阻、平板結(jié)構(gòu)尾跡控制等問題的理論分析與工程應(yīng)用提供理論支撐。