吳曉笛 劉華坪 陳浮

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

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

吳曉笛 劉華坪?陳浮

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

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

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

10.7498/aps.66.224702

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

1 引 言

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

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

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

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

2 數值方法

2.1 2.1 Chapman-Enskog展開分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中

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

2.2 流體控制方程

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

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

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

其中:

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

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

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

其中,網格界面通量重構與邊界處理是數值計算中兩個關鍵部分.

2.2.1 格子Boltzmann通量求解法

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

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

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

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

圖1 界面通量重構示意圖Fig.1.Local flux reconstruction at cell interface.

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

通量表示為

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

2.2.2 強制浸入邊界法

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

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

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

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

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

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

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

其中,

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

2.3 IB-MRT-LBFS算法流程

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

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

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

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

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

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

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

3 數值計算與分析

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

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

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

3.1 單圓柱橫向渦激振動

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

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

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

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

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

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

圖5 (網刊彩色)不同折合速度單圓柱渦激振動橫向升力系數與振幅的時程變化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 (網刊彩色)不同折合速度單圓柱橫向渦激振動瞬時渦量Fig.6.(color online)Instantaneous vorticity against reduced velocity for a VIV cylinder in transverse direction.

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

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

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

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

圖7 (網刊彩色)振幅及作用力系數隨雷諾數的變化情況 (a)橫向最大振幅;(b)軸向振幅均方根;(c)升力系數均方根;(d)阻力系數均方根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 (網刊彩色)不同雷諾數下單圓柱渦激振動瞬時渦量Fig.8.(color online)Instantaneous vorticity against Reynolds number for a VIV cylinder.

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

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

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

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

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

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

圖10 (網刊彩色)最大振幅隨折合速度的變化情況 (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 (網刊彩色)作用力系數隨折合速度的變化情況 (a)上游圓柱平均阻力系數;(b)上游圓柱阻力系數均方根;(c)上游圓柱升力系數均方根;(d)下游圓柱平均阻力系數;(e)下游圓柱阻力系數均方根;(f)下游圓柱升力系數均方根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)可知上游圓柱阻力系數的變化規(guī)律:當3≤Ur≤5時,阻力系數均值及均方根隨折合速度的增加而增大;當5＜Ur≤9時,阻力系數隨折合速度的增加而減小;隨著折合速度的繼續(xù)增加,阻力系數趨于平緩且趨近于單圓柱阻力.由圖11(c)可知上游圓柱升力系數隨折合速度的變化規(guī)律:升力系數均方根隨折合速度的增加呈現先增大后減小直至Ur=8;當Ur≥9時,升力系數均方根趨于平緩且趨近于單圓柱升力.由圖11(d)—(f)可知,下游圓柱所受作用力的變化規(guī)律不同于上游圓柱,這是由于下游圓柱受到來自上游圓柱不同結構尾跡形式的作用引起的.對比上下游圓柱所受阻力系數平均值,上游圓柱最大出現在Ur=5,而下游圓柱為Ur=7.

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

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

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

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

4 結 論

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

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

[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空間結構 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