剛性薄平板繞流特性的數值模擬

蔚杰, 武建軍

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

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

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

1 力學模型與數值方法

1.1 力學模型

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

圖1 計算域與邊界條件

1.2 控制方程

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

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

含外力項的LBE方程[20]:

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

(1)

(2)

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

(3)

平衡態分布函數為:

(4)

(5)

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

(6)

(7)

式中密度ρ0為常數。

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

u=u*+δu

(8)

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

(9)

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

(10)

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

(11)

(12)

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

(13)

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

(14)

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

(15)

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

(16)

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

fb=-2ρ0δuB/δt

(17)

由此可得結構的阻力、升力為:

(18)

進一步,阻力系數與升力系數分別被定義為[26]:

(19)

1.3 計算域離散

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

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

圖3 計算域分區策略與局部網格示意

1.4 數值驗證

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

圖4 圓柱繞流局部流線圖

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

表1 圓柱繞流Re=20, 40, 100阻力系數、斯特羅哈數計算結果

2 雷諾數對平板繞流的影響

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

模型的邊界條件為:

入口(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 雷諾數對平板尾流的影響

本文通過模擬得出,對于均勻來流條件下的平板繞流系統,由定常轉變為非定常狀態的臨界雷諾數11015,漩渦由P模式向具有更大的空間尺度的P+S模式轉變。二次渦街尚未充分發展的流向范圍內(如Re=400, 1420),可以明顯觀察到以2S漩渦模式為主的尾流結構;該區域內頻繁的出現失配單渦,這些失配單渦具有不規則的運動軌跡與空間分布特征。此外,尾流區出現一次躍遷、二次躍遷的流向位置具有時變特征;隨著Re的增大,出現躍遷的流向位置向上游移動。

圖5 瞬時渦量云圖

圖6 監測點u(t*)功率譜密度函數

2.2 雷諾數對平板阻力的影響

圖7 平板阻力系數功率譜密度

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

3 來流剪切參數對平板繞流的影響

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

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

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

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

3.1 來流剪切參數對平板尾流的影響

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

圖9 平板繞流系統St隨β變化規律

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

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

3.2 來流剪切參數對平板阻力的影響

圖12 平板阻力系數功率譜密度

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

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

4 結論

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

2)在主渦街與雙排渦內,漩渦運動軌跡是確定的,流動是周期性的。二次躍遷伴隨著漩渦失配-重組現象以及漩渦模式由P模式向大尺度的P+S或2S模式的轉變。二次渦街充分發展的區域,大尺度渦結構以及失配單渦的運動軌跡均是不規則的,引發流場的低頻、非周期性振蕩。

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

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

上述研究結果可為流體機械葉片設計、平板結構減阻、平板結構尾跡控制等問題的理論分析與工程應用提供理論支撐。