基于浸入邊界-多松弛時間格子玻爾茲曼通量求解法的流固耦合算法研究?

吳曉笛 劉華坪 陳浮

(哈爾濱工業(yè)大學能源科學與工程學院,哈爾濱 150001)

基于浸入邊界-多松弛時間格子玻爾茲曼通量求解法的流固耦合算法研究?

吳曉笛 劉華坪?陳浮

(哈爾濱工業(yè)大學能源科學與工程學院,哈爾濱 150001)

(2017年5月19日收到;2017年8月22日收到修改稿)

針對流固耦合問題,發(fā)展了基于浸入邊界-多松弛時間格子玻爾茲曼通量求解法(immersed boundary method multi-relaxation-time lattice Boltzmann flux solver,IB-MRT-LBFS)的弱耦合算法.依據(jù)多尺度Chapman-Enskog展開,建立不可壓宏觀方程狀態(tài)變量和通量與格子玻爾茲曼方程中粒子密度分布函數(shù)之間的關系;采用強制浸入邊界法處理流固界面使固壁表面滿足無滑移邊界條件,根據(jù)修正的速度求解動量方程力源項;結(jié)構(gòu)運動方程采用四階龍格-庫塔法求解.格子模型與浸入邊界法的引入使流固耦合計算可以在笛卡爾網(wǎng)格下進行,無需生成貼體網(wǎng)格及運用動網(wǎng)格技術,簡化了計算過程.數(shù)值模擬了單圓柱橫向渦激振動、單圓柱及串列雙圓柱雙自由度渦激振動問題.結(jié)果表明,IB-MRT-LBFS能夠準確預測圓柱渦激振動的鎖定區(qū)間、振動響應、受力情況以及捕捉尾流場結(jié)構(gòu)形態(tài),驗證了該算法在求解流固耦合問題的有效性和可行性.

10.7498/aps.66.224702

流固耦合,浸入邊界法,格子玻爾茲曼通量求解法,渦激振動

1 引 言

流固耦合問題廣泛存在于航空航天、海洋水利工程、石油運輸及生物工程等領域中.依據(jù)耦合機理流固耦合問題可以分為兩類[1]:流體、固體區(qū)域部分或完全重疊在一起,如滲流問題;流體和固體的耦合作用僅發(fā)生在兩種介質(zhì)交界面處,如渦激振動.依據(jù)控制方程解法流固耦合問題又可分為強耦合和弱耦合[2,3].強耦合將流體、固體域和耦合作用構(gòu)造在同一控制方程中直接求解;弱耦合是在單一時間步內(nèi)分別求解流體域和固體域,通過兩相介質(zhì)界面的數(shù)據(jù)交換實現(xiàn)流固耦合的計算,基于其概念簡單,易于編程實現(xiàn),因此得到廣泛應用[4].

目前,對于流固耦合這種動邊界問題常采用邊界貼合方法處理,即根據(jù)每個時間步物體運動變形生成一套新的貼體網(wǎng)格,該方法能夠很好地捕捉附面層特性,有利于模擬高雷諾數(shù)流動問題,但是網(wǎng)格生成的計算量將大大增加,降低整體程序的計算效率[5].近年來,非邊界貼合方法越來越受到關注.該類方法往往通過對邊界周圍流場變量的恰當處理來實現(xiàn)邊界條件,因此不嚴格要求網(wǎng)格與物體邊界貼體,可以采用直角網(wǎng)格等規(guī)則的網(wǎng)格,使得網(wǎng)格的生成更容易.當物體運動時,可以直接在固定的網(wǎng)格上求解流動控制方程,從而無需在每個時間步上對網(wǎng)格進行變換、映射、重新生成等處理,大大提高了計算效率[6].因此,非邊界貼合方法在模擬復雜流場和運動邊界的流動時能夠真正做到簡便高效.

浸入邊界法(immersed boundary method)作為非邊界貼合方法之一,在流固耦合問題中得到了廣泛的應用.該方法需要兩套網(wǎng)格,流場使用固定直角坐標網(wǎng)格來求解,而浸入邊界則用拉格朗日點來標識,物體與流場的作用通過兩套網(wǎng)格之間的信息交互傳遞來完成.傳統(tǒng)的浸入邊界法主要分為直接力法[7]、懲罰力法[8]以及動量變化法[9].王文全等[10,11]應用浸入邊界法實現(xiàn)了對剛體定軸轉(zhuǎn)動的預測及水輪機活動導葉旋轉(zhuǎn)擺動繞流的數(shù)值計算;明平劍等[12,13]采用浸入邊界法及非結(jié)構(gòu)化網(wǎng)格求解了低雷諾數(shù)下振蕩圓柱繞流問題;李師堯等[14]采用浸沒邊界-格子玻爾茲曼格式模擬貫流式水輪機三維瞬變流.

若采用浸入邊界法實現(xiàn)流固耦合問題的數(shù)值模擬,則需要一種有效的流場求解器.現(xiàn)階段流場計算主要是通過求解Navier-Strokes(N-S)宏觀方程或者格子玻爾茲曼方法(LBM).針對不可壓控制方程,由于連續(xù)方程中壓力變量的缺失,求解N-S方程的困難是對壓力場的預測,需引入壓力泊松方或者壓力修正方法來求解,會影響整體的計算效率.而LBM中壓力場可以由粒子分布函數(shù)直接給出,且具有算法簡單、易于并行的優(yōu)勢,但還存在一些缺點,由于格子速度模型的均勻性且網(wǎng)格尺度與時間推進步長的關聯(lián)性,造成該方法不能有效地計算工程應用問題,如計算湍流問題,邊界層問題等;LBM中松弛時間項與黏性系數(shù)有關,因而該方法只能求解黏性流動問題,而N-S方程既可以求解黏性流動又可以求解無黏流動.基于以上分析,Shu等[15]開發(fā)了一種基于單松弛時間格式的格子玻爾茲曼通量求解器,其思想是基于有限體積法空間離散宏觀方程,狀態(tài)變量及通量由格子玻爾茲曼模型中粒子密度分布函數(shù)重新構(gòu)造.該方法有效結(jié)合了兩種求解器的優(yōu)點.基于該思想,本文發(fā)展了一種基于多松弛時間格式的玻爾茲曼通量求解器(multi-relaxation-time lattice Boltzmann flux solver).將該求解流場的方法與Wang等[16,17]開發(fā)的強制浸入邊界法相結(jié)合,開發(fā)了一種基于浸入邊界-多松弛時間格子玻爾茲曼通量求解法(IB-MRT-LBFS)的弱耦合算法.該方法求解動邊界問題無需使用動網(wǎng)格技術,節(jié)約了網(wǎng)格重新生成的計算量;強制浸入邊界法的應用嚴格保證了無滑移邊界條件;多松弛格子玻爾茲曼求解法改善了傳統(tǒng)流體動力學計算方法必須花費大量時間求解泊松方程而得到壓力場的限制,并且通量表達式直接由平衡分布函數(shù)給出,簡化了流場的計算過程.

2 數(shù)值方法

2.1 2.1 Chapman-Enskog展開分析

在低馬赫數(shù)條件下,借助Chapman-Enskog展開可以完成宏觀Navier-Stokes方程與介觀Boltzmann方程的相互轉(zhuǎn)化,基本思想是將時空變量用多層次時空尺度表示,使其在各級尺度上物理量的量級保持一致.

宏觀不可壓Navier-Stokes控制方程表達式為:

其中,ρ,u,p和μ分別代表流體密度、速度、壓力及動力黏性系數(shù).

采用 Bhatnagar-Gross-Krook(BGK)碰撞模型的多松弛時間格子玻爾茲曼方程表示如下:

其中,α表示不同離散速度粒子;eα是粒子速度;fα是沿α方向的離散速度密度分布函數(shù);r代表物理位置;M是正交轉(zhuǎn)換矩陣,其作用是將速度空間密度分布函數(shù)fα及其平衡分布函數(shù)轉(zhuǎn)換為多松弛時間尺度下的動量空間分布函數(shù)mα=Mfα及其平衡分布函數(shù)S=diag(s0,s1,s2,s3,s4,s5,s6,s7,s8)為對角碰撞松弛矩陣,由于在碰撞過程中為了保證質(zhì)量和動量守恒,s0=s3=s5=0,s7=s8=1/τ,松弛時間τ是與動力黏性系數(shù)μ相關的物理量,取其中,cs為LBGK模型的聲速.

采用D2Q9格子玻爾茲曼模型,離散速度定義為

其中,c通常取值為1.平衡分布函數(shù)定義為

其中,w0=4/9,w1=w2=w3=w4=1/9,

宏觀方程中密度與動量可以用密度分布函數(shù)表示為

通過多尺度展開,分布函數(shù)、時間及空間導數(shù)可以表示為:

其中,ε為與克努森數(shù)成正比的參量.

將方程(3)進行二階泰勒展開,得到以下方程:

將方程(7)—(10)代入方程(11)得到

其中I為單位矩陣.根據(jù)方程(12)和(13)推導出:

通過對方程(13)和(14)零階矩和一階矩的求和得到以下方程:

其中

通過以上公式的推導,建立了宏觀方程通量與格子Boltzmann分布函數(shù)的關系,如(17)和(18)式所示.

2.2 流體控制方程

基于浸入邊界法的宏觀不可壓Navier-Stokes控制方程表達式為

其中F是由浸入邊界法求解的力源項.

根據(jù)2.1節(jié)中建立的宏觀方程狀態(tài)變量和通量與格子玻爾茲曼方程中密度分布函數(shù)之間的關系((6),(17),(18)式),采用IB-MRT-LBFS的不可壓控制方程表達式為:

其中:

求解流體控制方程(23)和(24)分為兩個步驟:

1)預測步,不考慮物面邊界,采用格子Boltzmann通量求解法(LBFS)求解流場變量,即密度ρn+1和瞬時速度u?,

2)修正步,采用強制浸入邊界法(IBM)求解修正后的速度變量un+1,進一步求解邊界作用力,

其中,網(wǎng)格界面通量重構(gòu)與邊界處理是數(shù)值計算中兩個關鍵部分.

2.2.1 格子Boltzmann通量求解法

為了求解方程(27)和(28)中ρn+1和u?,采用LBFS計算相鄰單元界面的通量.基于有限體積法的控制方程離散形式:

其中W=[ρ,ρu]T;dVi為控制體的體積,dSk為控制體第k個控制面,n為控制面上的法向向量.

基于格子玻爾茲曼D2Q9模型的界面通量重構(gòu)示意圖見圖1.已知相鄰兩個網(wǎng)格單元和界面的物理位置,分別定義為ri,ri+1以及r,則有

其中ψ代表流體變量ρ或u.

圖1 界面通量重構(gòu)示意圖Fig.1.Local flux reconstruction at cell interface.

基于泰勒展開計算圖1中離散點的速度,重構(gòu)界面處密度和動量,表達式為:

通量表示為

其中,平衡項可以根據(jù)界面處宏觀變量求解,而非平衡項可以用平衡項近似求解,表達式為:

2.2.2 強制浸入邊界法

浸入邊界法在數(shù)學建模過程中采用歐拉變量描述流場網(wǎng)格點,采用拉格朗日變量表示物面邊界,如圖2所示.應用強制浸入邊界法進行流場速度的修正,使得物面滿足無滑移邊界條件.u?為預測步求解的流場瞬時速度,?u為修正速度,則修正后的歐拉點速度表達式為

圖2 (網(wǎng)刊彩色)計算域二維示意圖Fig.2.(color online)A solid boundary immersed in a two-dimensional computational domain.

歐拉點上修正速度可以通過Dirac delta函數(shù)插值求解,表達式如下:

其中rj為歐拉點坐標;分別為拉格朗日點坐標以及修正速度;D是連續(xù)Kernel函數(shù),表達式為:

為了保證無滑移邊界條件,拉格朗日點的速度等價于同一位置流場中的速度,通過插值來求解拉格朗日點的流動速度,即

其中,l=1,···,N;i=1,···,M;h為歐拉點網(wǎng)格尺寸,N和M分別為拉格朗日點與歐拉點的個數(shù),則拉格朗日點上的速度為

將上述方程寫成矩陣的形式:

其中,

求解矩陣方程AX=B獲得拉格朗日點的修正速度通過(38)式求得歐拉點的修正速度?u,進而求得修正后的流場瞬時速度un+1及邊界作用力F.

2.3 IB-MRT-LBFS算法流程

IB-MRT-LBFS流固耦合算法步驟歸納如下:

1)給定初始狀態(tài)值,選取遷移時間步長δt,計算碰撞松弛矩陣S;

3)根據(jù)(32)和(33)式求解界面處流場宏觀變量,進而計算

4)根據(jù)(34)式計算界面通量,求解預測步中網(wǎng)格中心宏觀變量;

5)依據(jù)不同情況下剛體運動方程,計算修正步中拉格朗日點的運動速度及坐標解矩陣方程(43)獲得拉格朗日點修正速度;

6)根據(jù)(38)式求解流場修正速度,將其代入(37)式計算出修正后的流場宏觀速度變量;根據(jù)(29)式求解邊界作用力,完成修正步的計算;

7)重復步驟2)—6)至計算程序收斂.

3 數(shù)值計算與分析

利用IB-MRT-LBFS流固耦合算法開展單圓柱及串列雙圓柱渦激振動問題的數(shù)值模擬研究.設圓柱為彈性支撐,雙自由度無量綱運動方程表達式:

其中,X和Y分別為x,y方向上的無量綱位移,其一階導數(shù)和二階導數(shù)分別表示其相應方向上的無量綱速度和加速度;ζ為結(jié)構(gòu)阻尼系數(shù);質(zhì)量比mr=4m/(πρD2);折合固有頻率fr=fnD/U∞,fn為結(jié)構(gòu)固有頻率;折合速度Ur=U∞/(fnD);升阻力系數(shù){CL,CD}={FL,FD}/(0.5ρu2D).升阻力由應力積分法求解得表達式為

然而,當采用浸入邊界法處理運動剛體的物面條件時,由于剛體內(nèi)部受到內(nèi)部流體的作用力,因此應該考慮消除內(nèi)部作用力的影響[18],(49)式修正為

3.1 單圓柱橫向渦激振動

單圓柱鈍體結(jié)構(gòu)橫向渦激振動問題作為經(jīng)典驗證算例已有較多的文獻報道,其中Ahn和Kallinderis[19]采用ALE動網(wǎng)格技術數(shù)值模擬了流固耦合問題;Borazjani等[20]采用浸入邊界法中的直接力法求解曲線邊界的流固耦合問題;Jiang等[21]采用外插法處理運動邊界問題;Wang等[22]采用ALE動網(wǎng)格與雙重網(wǎng)格相結(jié)合的方法求解動邊界的問題.應用IB-MRT-LBFS模擬雷諾數(shù)為150的單圓柱單自由度渦激振動,取質(zhì)量比mr=8/π,為了使結(jié)構(gòu)振動幅度最大,阻尼系數(shù)設置為ξ=0,計算折合速度為3≤Ur≤8.計算域及邊界條件如圖3所示,計算域入口距圓柱中心距離為40D且為Dirichlet型邊界條件,密度為ρ=1.0,速度為U∞=0.1;上下邊界距圓柱中心位置均為40D且為無滑移邊界條件;出口距圓柱中心位置均為60D且為Neumann型邊界條件;圓柱表面為無滑移邊界條件,用均勻分布150個拉格朗日點代替,浸沒在160×240的均勻歐拉網(wǎng)格區(qū)域中;計算域總直角網(wǎng)格數(shù)為547×681.

圖3 單圓柱橫向渦激振動計算域及邊界條件示意圖Fig.3.The computational domain and boundary conditions for a VIV cylinder in transverse direction.

圖4給出了在不同折合速度下IB-MRT-LBFS流固耦合算法計算的無量綱橫向最大振幅與文獻數(shù)值結(jié)果的比較.由圖可知,本文結(jié)果與現(xiàn)有文獻[19—22]中采用不同數(shù)值方法的計算結(jié)果趨勢相一致.當Ur=4時,單圓柱單自由度渦激振動進入“鎖定區(qū)間”,尾流脫落渦頻率接近圓柱固有頻率,發(fā)生共振現(xiàn)象,振幅增加;當Ur=8時,單圓柱渦激振動超出“鎖定區(qū)間”,圓柱振幅降低.

圖4 (網(wǎng)刊彩色)不同折合速度單圓柱橫向渦激振動最大振幅Fig.4.(color online)Maximum amplitude against reduced velocity for a VIV cylinder in transverse direction.

不同的數(shù)值計算方法對單圓柱渦激振動橫向最大振幅結(jié)果影響較大,尤其是在非鎖定區(qū)間,例如在Ur=8時文獻[19,20]中計算結(jié)果相對誤差達到40%以上.由圖4可以看出,數(shù)值算法模擬結(jié)果介于不同文獻數(shù)值結(jié)果之間,且可以準確地捕捉鎖定區(qū)間范圍及預測共振區(qū)域的振幅值.表1給出了單圓柱橫向渦激振動最大振幅的誤差分析,可以看出本文數(shù)值結(jié)果與文獻結(jié)果相比較,多數(shù)誤差在10%以內(nèi),且鎖定區(qū)間內(nèi)的誤差較小,說明IB-MRT-LBFS流固耦合算法在計算渦激振動問題上的有效性.

圖5給出了不同折合速度下單圓柱橫向渦激振動升力系數(shù)與無量綱最大振幅隨時間歷程的變化情況.當Ur=3,4,5時,升力與振幅相位相同;而當Ur=6,7,8時,二者相位相反.圖6給出了不同折合速度單圓柱橫向渦激振動瞬時渦量圖,可以清晰地觀察到尾流場脫落渦結(jié)構(gòu)特性.當Ur=3時,尾流模態(tài)為單行渦街;當Ur=4時,尾流模態(tài)呈現(xiàn)不規(guī)律性;隨著折合速度進一步增加到Ur=5時,尾流模態(tài)出現(xiàn)雙行渦街;當Ur≥6時,尾流又恢復到單行渦街形態(tài).綜上所述,單圓柱橫向渦激振動響應規(guī)律以及尾流特性與文獻[23]中相一致,說明IB-MRT-LBFS流固耦合新型算法能準確的求解單圓柱單自由度渦激振動問題.

圖5 (網(wǎng)刊彩色)不同折合速度單圓柱渦激振動橫向升力系數(shù)與振幅的時程變化Fig.5.(color online)Lift coefficient and amplitude against reduced velocity for a VIV cylinder in transverse direction.

表1 單圓柱橫向渦激振動最大振幅誤差分析Table 1.The error analysis on maximum amplitude for a VIV cylinder in transverse direction.

圖6 (網(wǎng)刊彩色)不同折合速度單圓柱橫向渦激振動瞬時渦量Fig.6.(color online)Instantaneous vorticity against reduced velocity for a VIV cylinder in transverse direction.

3.2 單圓柱雙自由度渦激振動

采用IB-MRT-LBFS數(shù)值算法模擬60≤Re≤140條件下單圓柱雙自由度渦激振動,與文獻[24]中的計算參數(shù)相一致,取質(zhì)量比mr=10.0,阻尼系數(shù)ξ=0,計算域、邊界條件及網(wǎng)格與單圓柱橫向渦激振動相同.

圖7給出了在不同雷諾數(shù)條件下振幅及作用力的數(shù)值計算結(jié)果,且與文獻[24]結(jié)果相對比,可以看出各項結(jié)果符合良好;對比圖7(a)與圖7(b),可以看出橫向振幅值遠遠大于軸向振幅,當Re=90時,單圓柱進入鎖定區(qū)間,產(chǎn)生共振,振幅達到最大值;隨著雷諾數(shù)的增加,振幅值逐漸減小;當Re=130時,圓柱駛出鎖定區(qū)間.從圖7(c)與圖7(d)可以看出升阻力均方根值的變化情況,升力系數(shù)均方根值與橫向振幅的變化規(guī)律不同,首先隨雷諾數(shù)的增大而增加;當進入“鎖定區(qū)間”時(Re=90)達到最大值,隨之逐漸減小;當駛出“鎖定區(qū)間”時(Re=130)時,升力均方根值突然增大;阻力系數(shù)均方根值與軸向振幅的變化規(guī)律相一致.表2給出了橫向最大振幅的誤差結(jié)果,可以看出相對誤差結(jié)果均在10%以下,進一步說明IB-MRT-LBFS數(shù)值算法能準確地預測低雷諾數(shù)條件下圓柱渦激振動問題.

圖8給出了不同雷諾數(shù)下單圓柱渦激振動瞬時渦量圖.當Re=60,70和80時,尾流模態(tài)為單行渦街;當Re=90和100時,尾流模態(tài)為雙行渦街;當Re≥110時,尾流又恢復到單行渦街形態(tài);當雷諾數(shù)增加到140時尾流模態(tài)由單行渦街逐漸演化為雙行形態(tài).以上尾流場脫落渦結(jié)構(gòu)特性與文獻[24]中觀察到的尾流形態(tài)相一致,說明IB-MRT-LBFS能夠有效預測流場結(jié)構(gòu)特性.

圖7 (網(wǎng)刊彩色)振幅及作用力系數(shù)隨雷諾數(shù)的變化情況 (a)橫向最大振幅;(b)軸向振幅均方根;(c)升力系數(shù)均方根;(d)阻力系數(shù)均方根Fig.7.(color online)Amplitude and force coefficients against Reynolds number:(a)Max transverse amplitude;(b)the rms value of in-line amplitude;(c)the rms value o flift coefficient;(d)the rms value of drag coefficient.

表2 單圓柱雙自由度渦激振動橫向最大振幅誤差分析Table 2.The error analysis on maximum transverse amplitude for a VIV cylinder.

圖8 (網(wǎng)刊彩色)不同雷諾數(shù)下單圓柱渦激振動瞬時渦量Fig.8.(color online)Instantaneous vorticity against Reynolds number for a VIV cylinder.

3.3 串列雙圓柱雙自由度渦激振動

串列雙圓柱雙自由度渦激振動計算域及邊界條件如圖9所示,兩圓柱的質(zhì)量和直徑相同,圓心間距L/D=5;雷諾數(shù)為150;質(zhì)量比mr=8/π;阻尼系數(shù)ξ=0.邊界條件與單圓柱渦激振動相同,兩圓柱表面均為無滑移邊界條件,固壁邊界同樣用150個拉格朗日點代替,浸沒在640×320的均勻歐拉網(wǎng)格區(qū)域中;計算域總直角網(wǎng)格數(shù)為1060×740.分別計算折合速度3≤Ur≤11情況下串列雙圓柱雙自由度渦激振動問題.

圖9 串列雙圓柱渦激振動計算域及邊界條件示意圖Fig.9.The computational domain and boundary conditions for two VIV cylinders in tandem.

圖10給出了不同折合速度下串列雙圓柱系統(tǒng)上下游圓柱軸向及橫向最大振幅結(jié)果;圖11給出了不同折合速度下串列雙圓柱系統(tǒng)上下游圓柱受力情況.由圖可知,本文數(shù)值方法模擬結(jié)果與現(xiàn)有文獻[25,26]中結(jié)果符合良好,說明IB-MRT-LBFS流固耦合算法能夠準確地預測多鈍體結(jié)構(gòu)雙自由度渦激振動問題.

由圖10(a)可知,在折合速度3≤Ur≤11情況下,上游圓柱軸向振動最大振幅小于0.05,且在Ur=4時軸向位移最大;當Ur≥6時隨著折合速度的增加,軸向位移呈現(xiàn)先減小后增加的趨勢.由圖10(c)可知,下游圓柱軸向位移無顯著變化規(guī)律,這是由于上游圓柱形成的尾跡脫落渦結(jié)構(gòu)作用在下游圓柱表面,導致下游圓柱所受阻力有所變化.當Ur≥9時,下游圓柱軸向振幅突然增加,最大值達到0.34;說明在較高折合速度條件下,上游圓柱尾跡作用使得下游圓柱軸向振幅加大.

由圖10(b)可知,上游圓柱橫向振幅隨折合速度的增加先變大后減小直至趨于平緩;且最大振幅出現(xiàn)在Ur=5時,約為0.6;對比單圓柱橫向振動結(jié)果,上游圓柱的振動響應與單圓柱相似,但是對于串列雙圓柱系統(tǒng),在下游圓柱的作用下,上游圓柱的鎖定區(qū)域有所增加.由圖10(d)可知,下游圓柱橫向振幅隨折合速度變化趨勢與上游圓柱基本相同;當Ur≥6時下游圓柱在上游尾跡的作用下振幅增加,均大于上游圓柱振幅,這是由于此折合速度下,下游圓柱升力均方根值大于上游圓柱(見圖11(c)和圖11(f));最大振幅出現(xiàn)在Ur=8時,約為1.06;對于串列雙圓柱系統(tǒng),上游圓柱的作用推遲下游圓柱進入鎖定區(qū)間,且擴大了下游圓柱鎖定區(qū)間的范圍.

圖10 (網(wǎng)刊彩色)最大振幅隨折合速度的變化情況 (a)上游圓柱(軸向);(b)上游圓柱(橫向);(c)下游圓柱(軸向);(d)下游圓柱(橫向)Fig.10.(color online)Max amplitude against reduced velocity:(a)Upstream cylinder(in-line direction);(b)upstream cylinder(transverse direction);(c)downstream cylinder(in-line direction);(d)downstream cylinder(transverse direction).

圖11 (網(wǎng)刊彩色)作用力系數(shù)隨折合速度的變化情況 (a)上游圓柱平均阻力系數(shù);(b)上游圓柱阻力系數(shù)均方根;(c)上游圓柱升力系數(shù)均方根;(d)下游圓柱平均阻力系數(shù);(e)下游圓柱阻力系數(shù)均方根;(f)下游圓柱升力系數(shù)均方根Fig.11.(color online)Force coefficients against reduced velocity:(a)The mean drag coefficient,upstream cylinder;(b)the rms value of drag coefficient,upstream cylinder;(c)the rms value of lift coefficient,upstream cylinder;(d)the mean drag coefficient,downstream cylinder;(e)the rms value of drag coefficient,downstream cylinder;(f)the rms value of lift coefficient,downstream cylinder.

由圖11(a)和圖11(b)可知上游圓柱阻力系數(shù)的變化規(guī)律:當3≤Ur≤5時,阻力系數(shù)均值及均方根隨折合速度的增加而增大;當5＜Ur≤9時,阻力系數(shù)隨折合速度的增加而減小;隨著折合速度的繼續(xù)增加,阻力系數(shù)趨于平緩且趨近于單圓柱阻力.由圖11(c)可知上游圓柱升力系數(shù)隨折合速度的變化規(guī)律:升力系數(shù)均方根隨折合速度的增加呈現(xiàn)先增大后減小直至Ur=8;當Ur≥9時,升力系數(shù)均方根趨于平緩且趨近于單圓柱升力.由圖11(d)—(f)可知,下游圓柱所受作用力的變化規(guī)律不同于上游圓柱,這是由于下游圓柱受到來自上游圓柱不同結(jié)構(gòu)尾跡形式的作用引起的.對比上下游圓柱所受阻力系數(shù)平均值,上游圓柱最大出現(xiàn)在Ur=5,而下游圓柱為Ur=7.

通過圖12不同折合速度串列雙圓柱渦激振動瞬時渦量可以觀察串列雙圓柱尾跡模態(tài).在圓心間距L/D=5的情況下,上游圓柱在不同折合速度條件下均可以形成周期性的脫落渦結(jié)構(gòu),但是由于折合速度與脫落渦頻率有關,因此形成的脫落渦的頻率與渦長有所改變,作用在下游圓柱表面的渦結(jié)構(gòu)不同進而影響下游圓柱表面作用力的分布及大小;觀察下游圓柱尾跡結(jié)構(gòu),當Ur=3時,下游圓柱尾跡為雙行渦街;當Ur=4時下游圓柱尾跡模態(tài)為單行渦街;當Ur=5,6時尾流模態(tài)呈現(xiàn)不規(guī)律性;隨著折合速度進一步的增加,尾流模態(tài)出現(xiàn)雙行渦街,當Ur=11時尾跡結(jié)構(gòu)呈不規(guī)律性.

圖13為串列雙圓柱系統(tǒng)上下游圓柱位移軌跡對比結(jié)果;紅色曲線表示上游圓柱運動軌跡,黑色曲線表示下游圓柱運動軌跡.由圖可知,當Ur=3時上下游圓柱運動軌跡均呈現(xiàn)“8”字形,在此折合速度條件下,上游圓柱未進入“鎖頻區(qū)”,振動較弱,而下游圓柱在上游圓柱形成的周期性尾跡的作用下振動有所加強,雙自由度的幅值均大于上游圓柱;當Ur=4時,上下游圓柱運動軌跡均呈現(xiàn)不規(guī)律性,在此折合速度條件下,上游圓柱進入“鎖頻區(qū)”,發(fā)生共振現(xiàn)象,振幅增加,此時上游圓柱脫落渦接近圓柱固有頻率,因此下游圓柱在尾跡的作用下振動加強;隨著折合速度的增加至Ur=8,上游圓柱在“鎖頻區(qū)”且均呈現(xiàn)“8”字形軌跡;當Ur=7,8時,下游圓柱軌跡為“8”字形;隨著折合速度進一步增加,上游圓柱駛出“鎖頻區(qū)”,振動減弱,而下游圓柱在上游圓柱尾跡的作用下保持較強的振動狀態(tài),運動軌跡呈現(xiàn)非規(guī)律性.

圖12 (網(wǎng)刊彩色)不同折合速度串列雙圓柱渦激振動瞬時渦量Fig.12.(color online)Instantaneous vorticity against reduced velocity for two VIV cylinders in tandem.

圖13 (網(wǎng)刊彩色)串列雙圓柱上下游圓柱位移軌跡Fig.13.(color online)Displacement trajectories of the upstream and downstream cylinder centers with various reduced velocities.

4 結(jié) 論

發(fā)展了一種基于浸入邊界-多松弛時間格式格子玻爾茲曼通量求解法的弱耦合流固算法(IBMRT-LBFS).對于流體控制方程,基于多尺度Chapman-Enskog展開,建立宏觀方程狀態(tài)變量及通量與格子玻爾茲曼方程中密度分布函數(shù)之間的關系,開發(fā)了多松弛時間格式的格子玻爾茲曼通量求解器;在流體方程中引入力源項,在保證無滑移邊界條件下采用強制浸入邊界法求解邊界作用力大小;對于結(jié)構(gòu)振動方程,采用四階龍格-庫塔法求解.該算法采用有限體積法進行空間離散,使得計算過程可以在非均勻網(wǎng)格下進行,對比傳統(tǒng)的格子Boltzmann方法減少了網(wǎng)格量;控制方程中狀態(tài)變量及壓力項由粒子密度分布函數(shù)直接給出,網(wǎng)格界面通量重構(gòu)使得宏觀方程中對流項與黏性項可以同時由平衡分布函數(shù)求得,通量求解過程得以簡化;強制邊界法處理運動邊界問題時,使得計算網(wǎng)格是固定的,無需再每一時間步更新網(wǎng)格的尺寸及位置,減少網(wǎng)格重新生成的計算時間.

數(shù)值計算和分析單圓柱橫向渦激振動、單圓柱及串列雙圓柱雙自由度渦激振動問題.對于單圓柱渦激振動,橫向振動響應大于軸向振動,當振動發(fā)生在鎖定區(qū)間,尾流脫落渦頻率接近圓柱固有頻率,產(chǎn)生共振,振幅較大.對于串列雙圓柱渦激振動,在較大折合速度條件下,下游圓柱振動響應大于上游圓柱,上游圓柱的作用力接近于單圓柱振動情況;上游圓柱推遲下游圓柱進入鎖定區(qū)間且增加區(qū)間范圍.通過對比現(xiàn)有文獻數(shù)值結(jié)果,表明IB-MRT-LBFS能夠準確預測圓柱渦激振動的鎖定區(qū)間、振動響應、受力情況以及捕捉尾流場結(jié)構(gòu)形態(tài),驗證了文中所提出的數(shù)值算法在計算流固耦合問題上的準確性與可行性.

[1]Xing J T,Zhou S,Cui E J 1997Adv.Mech.27 19(in Chinese)[邢景棠,周盛,崔爾杰 1997力學進展 27 19]

[2]Qian R J,Dong S L,Yuan X F 2008Spat.Struct.14 3(in Chinese)[錢若軍,董石麟,袁行飛 2008空間結(jié)構(gòu) 14 3]

[3]Guo P,Liu J,Wu W H 2013Chin.J.Theor.Appl.Mech.45 283(in Chinese)[郭攀,劉君,武文華 2013力學學報45 283]

[4]Zhou D,He T,Tu J H 2012Chin.J.Theor.Appl.Mech.44 494(in Chinese)[周岱,何濤,涂佳黃 2012力學學報44 494]

[5]Zhong G H,Liang A,Sun X F 2007J.Eng.Thermophys.28 399(in Chinese)[鐘國華,梁岸,孫曉峰 2007工程熱物理學報28 399]

[6]Liu Q Y 2012M.S.Dissertation(Nanjing:Nanjing University of Aeronautics and Astronautics)(in Chinese)[劉齊迎2012碩士學位論文(南京:南京航空航天大學]

[7]Luo H X,Dai H,Ferreira D S,Paulo J S A,Yin B 2012Comput.Fluids56 61

[8]Feng Z,Michaelides E 2004J.Comput.Phys.195 602

[9]Chen Y,Cai Q D,Xia Z H,Wang M,Chen S Y 2013Phys.Rev.E88 013303

[10]Wang W Q,Zhang G W,Yan Y 2017Trans.Beijing Inst.Tech.37 151(in Chinese)[王文全,張國威,閆妍2017北京理工大學學報37 151]

[11]Wang W Q,Su S Q,Yan Y 2015Chin.J.Comput.Mech.32 560(in Chinese)[王文全,蘇仕琪,閆妍2015計算力學學報32 560]

[12]Ming P J,Zhang W P 2009Chin.J.Aeronaut.22 480

[13]Ming P J,Zhang W P,Lu X Q,Zhu M G 2010Chin.J.Hydrodyn.25 321(in Chinese)[明平劍,張文平,盧熙群,朱明剛2010水動力研究與進展25 321]

[14]Li S Y,Cheng Y G,Zhang C Z 2016J.Huazhong Univ.Sci.Tech.(Natural Science Edition)44 122(in Chinese)[李師堯,程永光,張春澤2016華中科技大學學報(自然科學版)44 122]

[15]Shu C,Wang Y,Teo C J,Wu J 2014Adv.Appl.Math.Mech.6 436

[16]Wang Y,Shu C,Teo C J,Wu J 2015J.Fluids Struct.54 440

[17]Wang Y,Shu C,Yang L M,Sun Y 2017Int.J.Numer.Meth.Fluids83 331

[18]Suzuki K,Inamuro T 2011Comput.Fluids49 173

[19]Ahn H T,Kallinderis Y 2006J.Comput.Phys.219 671

[20]Borazjani I,Ge L,Sotiropoulos F 2008J.Comput.Phys.227 7587

[21]Jiang R J,Lin J Z,Chen Z L 2013Phys.Rev.E88 023009

[22]Wang C L,Tang H,Duan F,Yu S C M 2016J.Fluids Struct.60 160

[23]Han Z L,Zhou D,Tu J H 2014J.Eng.Mech140 04014059

[24]Prasanth A K,Mittal S 2008J.Fluids Mech.594 463

[25]Bao Y,Huang C,Zhou D,Tu J H,Han Z L 2012J.Fluids Struct.35 50

[26]Yu K R,Etienne S Scolan,Y M,Hay A,Fontaine E,Pelletier D 2016J.Fluids Struct.60 37

PACS:47.63.mf,47.11.–j,47.11.Qr,43.40.–rDOI:10.7498/aps.66.224702

*Project supported by the National Natural Science Foundation of China(Grant No.51306042).

?Corresponding author.E-mail:hgdlhp@163.com

A method combined immersed boundary with multi-relaxation-time lattice Boltzmann flux solver for fluid-structure interaction?

Wu Xiao-DiLiu Hua-Ping?Chen Fu

(School of Energy Science and Engineering,Harbin Institute of Technology,Harbin 150001,China)

19 May 2017;revised manuscript

22 August 2017)

This paper performs a newly developed method,which combines the immersed boundary method(IBM)with multi-relaxation-time lattice Boltzmann flux solver(MRT-LBFS),for solving fluid-structure interaction problems.Finite volume discretization is used to solve the macroscopic governing equations with the flow variables de fined at cell centers.Based on the multi-scale Chapman-Enskog expansion analysis,LBFS builds a relationship between the variables and fluxes in incompressible Navier-Stokes equations and density distribution functions in lattice Boltzmann equation.In order to ensure no-slip boundary condition,boundary condition-enforced immersed boundary method is used to treat the fluid-structure interface.The restoring force can be resolved by making a velocity correction in the flow field.The four-stage RungeKutta scheme is used to solve the motion equation of structure.Using the lattice model and immersed boundary method, fluid-structure coupling calculation can be implemented in a Cartesian grid,without generating the body- fitted mesh and using moving mesh technique.Therefore,the computational process is considerably simpli fied.In order to verify the validity and feasibility of IB-MRT-LBFS to solve fluid-structure interaction problems,both oneand two-degree of freedom vortex-induced vibrations(VIV)of a circular cylinder and two-degree of freedom VIV of two cylinders in a tandem arrangement are simulated by this proposal method.For a VIV cylinder system,the transverse vibration response is much stronger than the axial response.When the vibration occurs in the range of lock-in regime,the shedding vortex frequency of the wake is close to natural frequency of the cylinder so that resonance appears,consequently causing larger amplitude.For two VIV cylinders in a tandem arrangement,the dynamic behavior of each cylinder is signi ficantly di ff erent from that of a single cylinder.The gap spacing between the two cylinder centers is a signi ficant parameter which e ff ects vibration characteristics and the spacing is fixed in the simulations of two tandem cylinders.With the e ff ects of upstream cylinder wake,the axial and transverse amplitudes of downstream cylinder obviously increase with adding the reduced velocity.The downstream cylinder is delayed,coming into lock-in regime,and the range of lock-in regime is expanded under the e ff ects of the wake of the upstream cylinder.As the reduced velocity is relatively large,the vibration response of the upstream cylinder is close to a single cylinder and the vibration response of the downstream cylinder is more intense than the upstream cylinder.Compared with the existing literature results,our result illustrates that IB-MRT-LBFS owns the ability to correctly predict the lock-in regime,dynamic response and the forces of vortex-induced vibrations of cylinders.And this method can accurately capture the wake structures.

fluid-structure interaction,immersed boundary method,lattice Boltzmann flux solver,vortex-induced vibration

?國家自然科學基金(批準號:51306042)資助的課題.

?通信作者.E-mail:hgdlhp@163.com