吳曉笛 劉華坪 陳浮

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

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

吳曉笛 劉華坪?陳浮

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

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

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

10.7498/aps.66.224702

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

1 引 言

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

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

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

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

2 數(shù)值方法

2.1 2.1 Chapman-Enskog展開分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中

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

2.2 流體控制方程

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

其中F是由浸入邊界法求解的力源項(xiàng).

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

其中:

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

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

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

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

2.2.1 格子Boltzmann通量求解法

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

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

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

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

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

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

通量表示為

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

2.2.2 強(qiáng)制浸入邊界法

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

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

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

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

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

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

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

其中,

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

2.3 IB-MRT-LBFS算法流程

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

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

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

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

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

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

7)重復(fù)步驟2)—6)至計(jì)算程序收斂.

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

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

其中,X和Y分別為x,y方向上的無量綱位移,其一階導(dǎo)數(shù)和二階導(dǎo)數(shù)分別表示其相應(yīng)方向上的無量綱速度和加速度;ζ為結(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).升阻力由應(yīng)力積分法求解得表達(dá)式為

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

3.1 單圓柱橫向渦激振動

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

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

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

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

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

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

圖5 (網(wǎng)刊彩色)不同折合速度單圓柱渦激振動橫向升力系數(shù)與振幅的時(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)刊彩色)不同折合速度單圓柱橫向渦激振動瞬時(shí)渦量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條件下單圓柱雙自由度渦激振動,與文獻(xiàn)[24]中的計(jì)算參數(shù)相一致,取質(zhì)量比mr=10.0,阻尼系數(shù)ξ=0,計(jì)算域、邊界條件及網(wǎng)格與單圓柱橫向渦激振動相同.

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

圖8給出了不同雷諾數(shù)下單圓柱渦激振動瞬時(shí)渦量圖.當(dāng)Re=60,70和80時(shí),尾流模態(tài)為單行渦街;當(dāng)Re=90和100時(shí),尾流模態(tài)為雙行渦街;當(dāng)Re≥110時(shí),尾流又恢復(fù)到單行渦街形態(tài);當(dāng)雷諾數(shù)增加到140時(shí)尾流模態(tài)由單行渦街逐漸演化為雙行形態(tài).以上尾流場脫落渦結(jié)構(gòu)特性與文獻(xiàn)[24]中觀察到的尾流形態(tài)相一致,說明IB-MRT-LBFS能夠有效預(yù)測流場結(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ù)下單圓柱渦激振動瞬時(shí)渦量Fig.8.(color online)Instantaneous vorticity against Reynolds number for a VIV cylinder.

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

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

圖9 串列雙圓柱渦激振動計(jì)算域及邊界條件示意圖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)有文獻(xiàn)[25,26]中結(jié)果符合良好,說明IB-MRT-LBFS流固耦合算法能夠準(zhǔn)確地預(yù)測多鈍體結(jié)構(gòu)雙自由度渦激振動問題.

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

由圖10(b)可知,上游圓柱橫向振幅隨折合速度的增加先變大后減小直至趨于平緩;且最大振幅出現(xiàn)在Ur=5時(shí),約為0.6;對比單圓柱橫向振動結(jié)果,上游圓柱的振動響應(yīng)與單圓柱相似,但是對于串列雙圓柱系統(tǒng),在下游圓柱的作用下,上游圓柱的鎖定區(qū)域有所增加.由圖10(d)可知,下游圓柱橫向振幅隨折合速度變化趨勢與上游圓柱基本相同;當(dāng)Ur≥6時(shí)下游圓柱在上游尾跡的作用下振幅增加,均大于上游圓柱振幅,這是由于此折合速度下,下游圓柱升力均方根值大于上游圓柱(見圖11(c)和圖11(f));最大振幅出現(xiàn)在Ur=8時(shí),約為1.06;對于串列雙圓柱系統(tǒng),上游圓柱的作用推遲下游圓柱進(jìn)入鎖定區(qū)間,且擴(kuò)大了下游圓柱鎖定區(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ī)律:當(dāng)3≤Ur≤5時(shí),阻力系數(shù)均值及均方根隨折合速度的增加而增大;當(dāng)5＜Ur≤9時(shí),阻力系數(shù)隨折合速度的增加而減小;隨著折合速度的繼續(xù)增加,阻力系數(shù)趨于平緩且趨近于單圓柱阻力.由圖11(c)可知上游圓柱升力系數(shù)隨折合速度的變化規(guī)律:升力系數(shù)均方根隨折合速度的增加呈現(xiàn)先增大后減小直至Ur=8;當(dāng)Ur≥9時(shí),升力系數(shù)均方根趨于平緩且趨近于單圓柱升力.由圖11(d)—(f)可知,下游圓柱所受作用力的變化規(guī)律不同于上游圓柱,這是由于下游圓柱受到來自上游圓柱不同結(jié)構(gòu)尾跡形式的作用引起的.對比上下游圓柱所受阻力系數(shù)平均值,上游圓柱最大出現(xiàn)在Ur=5,而下游圓柱為Ur=7.

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

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

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

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

[1]Xing J T,Zhou S,Cui E J 1997Adv.Mech.27 19(in Chinese)[邢景棠,周盛,崔爾杰 1997力學(xué)進(jìn)展 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力學(xué)學(xué)報(bào)45 283]

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

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

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

[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北京理工大學(xué)學(xué)報(bào)37 151]

[11]Wang W Q,Su S Q,Yan Y 2015Chin.J.Comput.Mech.32 560(in Chinese)[王文全,蘇仕琪,閆妍2015計(jì)算力學(xué)學(xué)報(bào)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水動力研究與進(jìn)展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華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版)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

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

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