周 凱，王 震，陳維山，龍曉軍

（1.山東農(nóng)業(yè)大學(xué) 機(jī)械與電子工程學(xué)院，山東 泰安 271018；2.哈爾濱工業(yè)大學(xué) 機(jī)器人技術(shù)與系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150001）

0 引 言

格子Boltzmann方法是一種不同于傳統(tǒng)數(shù)值方法的流體建模和計算方法[1-3]。傳統(tǒng)的數(shù)值方法的基本思路是將宏觀控制方程進(jìn)行離散，然后利用數(shù)值求解的方法去求解該方程，例如有限差分、有限體積、有限元法等。格子Boltzmann方法的是基于分子運(yùn)動學(xué)和統(tǒng)計力學(xué)，在微觀粒子尺度上建立離散速度模型。在滿足質(zhì)量、動量和能量守恒的基本條件下，通過粒子遷移和碰撞兩個微觀行為來反應(yīng)宏觀特性，這使得算法簡化為對粒子遷移和碰撞的模擬。格子Boltzmann方法采用均勻的正交網(wǎng)格，使得網(wǎng)格算法的工作量大大簡化，特別是在處理復(fù)雜曲邊界時，結(jié)合更為準(zhǔn)確的邊界處理算法和多塊網(wǎng)格耦合算法，這種優(yōu)勢更加明顯。

多個圓柱的繞流問題在工程中是很常見的，例如海洋平臺支撐柱、橋墩和水底管路與流體之間的作用以及它們之間的相互作用。圓柱之間相互作用的存在使得此類問題比單圓柱繞流的情況要復(fù)雜得多。這類問題的早期研究多側(cè)重于對實(shí)驗(yàn)現(xiàn)象的觀察和圓柱力學(xué)參數(shù)的測量。近代以來數(shù)值模擬方法得到了巨大發(fā)展，此類問題可以用數(shù)值模擬的方法進(jìn)行研究，并取得了一定研究成果。

Zdravkovich[4-5]等對不同排列方式的雙圓柱繞流問題進(jìn)行了深入的實(shí)驗(yàn)研究。針對串列圓柱，研究發(fā)現(xiàn)了臨界間距的存在。當(dāng)兩圓柱間距小于該臨界值時，上游圓柱沒有明顯的脫渦現(xiàn)象。這一臨界值約為3.5倍直徑。上下游圓柱尾流場速度分布在該臨界間距附近突然發(fā)生改變，上游圓柱出現(xiàn)渦脫落現(xiàn)象，兩柱之間速度和尾流速度都變大。隨著流體力學(xué)理論的完善和計算機(jī)技術(shù)的發(fā)展，計算流體力學(xué)成為解決流動問題的重要工具。對串列圓柱繞流的數(shù)值模擬發(fā)現(xiàn)了與實(shí)驗(yàn)現(xiàn)象一致的結(jié)果，同時獲得了更多的流場細(xì)節(jié)，模擬結(jié)果也更接近實(shí)驗(yàn)結(jié)果[6]。

文章將基于格子Boltzmann方法，對流體的流動和流固耦合作用進(jìn)行數(shù)值模擬。對邊界的處理，采用經(jīng)典的粒子反彈格式，即流體粒子在碰到固體邊界后速度反向。針對格子Boltzmann方法的均勻規(guī)則網(wǎng)格在處理曲邊界問題時的缺陷，采用更為精確的邊界處理算法，讓流體粒子與固體的碰撞反彈過程在固體格點(diǎn)處發(fā)生。為提高計算效率，滿足固體周圍相對較高的網(wǎng)格密度要求，開發(fā)分塊網(wǎng)格耦合算法。以單個靜止圓柱繞流的已有研究結(jié)果為依據(jù)，對本文的數(shù)值算法進(jìn)行驗(yàn)證。最后，對不同間距下靜止串列雙圓柱繞流問題進(jìn)行數(shù)值模擬，研究不同情況下圓柱之間的相互作用以及尾流特征的變化。

1 數(shù)值方法

1.1 格子Boltzmann方法

在格子Boltzmann方法中，粒子的運(yùn)動由波爾茲曼方程描述

其中：u 為流速速度，c=δx/δt，ωα是加權(quán)系數(shù)，

最終，可以得到粒子的動力學(xué)演化方程：

（1）碰撞過程：

（2）遷移過程：

粒子的運(yùn)動分解為碰撞和遷移兩個過程。在遷移過程中，粒子沿eα方向運(yùn)動，與相鄰網(wǎng)格點(diǎn)進(jìn)行數(shù)據(jù)交換；在碰撞過程中，粒子與來自不同方向的粒子進(jìn)行碰撞，在滿足質(zhì)量守恒、動量守恒和能量守恒的條件下，重新構(gòu)造粒子分布函數(shù)。通過（5）式可以看出，粒子的碰撞方程是顯式可解的，遷移過程的計算量也較小，算法的編制相對簡單。

在粒子的分布函數(shù)確定后，流場的宏觀量可以表示為

本文采用經(jīng)典二維九速度模型（D2Q9），將粒子速度離散為9個方向，即α的值有9個，如圖1所示。不同方向的速度矢量表達(dá)式為

圖1 二維9速度（D2Q9）模型Fig.1 The 2D,nine-velocity lattice(D2Q9)model

1.2 多塊網(wǎng)格耦合算法

不同于傳統(tǒng)的貼體網(wǎng)格處理方法，格子Boltzmann通常采用均勻結(jié)構(gòu)網(wǎng)格來劃分流場，網(wǎng)格算法相對更容易實(shí)現(xiàn)。然而，在處理流固耦合問題時，固體壁面周圍對網(wǎng)格的密度要求較高，而遠(yuǎn)離壁面處這種要求相對較低，均勻網(wǎng)格顯然會造成計算量的浪費(fèi)。因此，研究多塊網(wǎng)格的耦合算法是必要的。本文在保證質(zhì)量連續(xù)和應(yīng)力連續(xù)前提下，采用了多塊網(wǎng)格耦合的方法來處理流固耦合中的網(wǎng)格問題。

為簡單起見，文章以兩塊網(wǎng)格的耦合為例進(jìn)行分析。粗細(xì)網(wǎng)格的尺寸比m，而且在網(wǎng)格塊的計算中滿足 δx=δy=δt，粗細(xì)網(wǎng)格交界處的網(wǎng)格結(jié)構(gòu)如圖2所示。

前面已經(jīng)提到，運(yùn)動粘度和無量綱松弛時間應(yīng)滿足（8）式。根據(jù)（8）式，為保證在交界處運(yùn)動粘度的連續(xù)，不同網(wǎng)格塊的無量綱松弛時間應(yīng)該滿足（9）式。

圖2 不同網(wǎng)格塊交界處網(wǎng)格結(jié)構(gòu)Fig.2 Interface stucture between two blocks of different spacing

其中：下標(biāo)f表示細(xì)網(wǎng)格，c表示粗網(wǎng)格?？梢詫⑺俣确植己瘮?shù)分解為平衡部分和非平衡部分因?yàn)樵诰W(wǎng)格交界面上必須滿足速度和密度的連續(xù)，因此

圖3 網(wǎng)格耦合算法計算流程圖Fig.3 Flow chart of the computational procedure in the multi-block method

文章采用的多塊網(wǎng)格尺寸比m=2，相鄰網(wǎng)格會有重合。粗網(wǎng)格和細(xì)網(wǎng)格有部分節(jié)點(diǎn)是重合的，重合節(jié)點(diǎn)處可以利用（18）式和（19）式直接進(jìn)行數(shù)據(jù)的交換，其他未知點(diǎn)采用多點(diǎn)插值的方法得出，不同的插值方式對計算結(jié)果將產(chǎn)生不同的影響。計算流程如圖3所示，n表示當(dāng)前的時間步，箭頭后表示的是該步驟后得到的計算參數(shù)。

1.3 曲邊界處理方法和作用力的提取

在格子Boltzmann方法中，邊界條件的處理方法將直接關(guān)系到計算的精度和穩(wěn)定度。對于靜止的固體邊界，最常見的處理方法是切向速度為零的無滑移邊界條件。具體到格子Boltzmann方法中，入射粒子在碰到邊界后，速度方向變?yōu)橄喾捶较?，這種處理方法也被稱為反彈格式。反彈格式在處理靜止邊界問題上可以得到較好的效果。但是，由于格子Boltzmann方法多采用均勻結(jié)構(gòu)網(wǎng)格，邊界處會出現(xiàn)鋸齒，只能采用加密網(wǎng)格的方法來提高對邊界的擬合精度，因此無法準(zhǔn)確擬合曲邊界。郭照立等[8-9]結(jié)合平直邊界的非平衡態(tài)外推格式，得到了一種曲邊界處理格式，這種格式通過在固體格點(diǎn)位置處執(zhí)行碰撞過程來確定碰撞后的分布函數(shù)，可以獲得較好的效果。

曲邊界處網(wǎng)格結(jié)構(gòu)如圖4所示。以s點(diǎn)代表固體格點(diǎn)，f代表流體格點(diǎn)，w代表固體邊界點(diǎn)。本文所采用的方法假定碰撞在s處執(zhí)行，分析以一個方向的處理為例，其他方向的處理方法類似。在固體點(diǎn)s處，將碰撞前的分布函數(shù)分解為平衡態(tài)feq和非平衡態(tài)fneq兩部分

圖4 曲邊界示意圖Fig.4 Schematic description of curved wall boundary

其中，平衡態(tài)部分采用粒子的平衡分布函數(shù)近似

對于（21）式中的非平衡部分，可表示為

至此，碰撞前固體點(diǎn)的分布函數(shù)就確定了。然后，根據(jù)LBGK碰撞模型進(jìn)行標(biāo)準(zhǔn)碰撞，可以得到碰撞后的分布函數(shù)

在處理流固耦合問題中，固體的受力參數(shù)往往是問題的目標(biāo)參數(shù)，因此流體與固體相互作用力的提取是必要的。本文采用動量交換法來計算流體與圓柱之間的相互作用力，即通過計算流體粒子與壁面發(fā)生碰撞前后的動量變化量，根據(jù)動量定理來計算粒子所受的流體作用力。在格子Boltzmann方法中，在與壁面碰撞前粒子的動量可以表示為

與壁面碰撞后粒子速度反向，其動量可以表示為

其中的負(fù)號表示與α方向的相反方向。

由此，可以計算出流體粒子所受的作用力。根據(jù)經(jīng)典牛頓力學(xué)中相互作用力的性質(zhì)，二者大小相

同而方向相反。因此，通過作用力沿壁面的積分就可以得到固體所受的作用力，具體形式為

進(jìn)而可以得出圓柱的阻力系數(shù)Cd和升力系數(shù)Cl：

其中：FD和FL分別是圓柱所受的流體作用力在順流方向和橫流方向的分量，U為來流速度。

2 數(shù)值方法驗(yàn)證

為驗(yàn)證算法的可靠性，本文在耦合網(wǎng)格下對Re=200靜止單圓柱繞流問題進(jìn)行了數(shù)值模擬，其中雷諾數(shù)Re=DU/v，D為圓柱直徑，U為來流速度，v是流體的運(yùn)動粘度。計算區(qū)域的劃分和網(wǎng)格劃分情況如圖5所示。在網(wǎng)格密度要求較高的圓柱周圍采用尺寸較小的網(wǎng)格塊，其他區(qū)域采用尺寸相對較大的網(wǎng)格塊。 細(xì)網(wǎng)格區(qū)域取即一倍直徑D分為50個格子，粗網(wǎng)格區(qū)域網(wǎng)格尺寸是細(xì)網(wǎng)格區(qū)域的2倍，即網(wǎng)格尺寸比m=2。網(wǎng)格塊尺寸間距和時間步大小的確定兼顧了計算的精度和效率。靜止單個圓柱繞流問題已經(jīng)有較多的理論研究、實(shí)驗(yàn)測量和數(shù)值模擬研究成果，通過與已有結(jié)果對比可以驗(yàn)證本文算法的可靠性。

本文數(shù)值模擬結(jié)果如圖6所示。圓柱表面發(fā)生了周期性的脫渦現(xiàn)象，正負(fù)渦從圓柱上下表面交替脫落。由圓柱升阻力系數(shù)曲線可以看出，其升阻力系數(shù)曲線呈現(xiàn)周期性震蕩，Cd的頻率是Cl頻率的兩倍，這是因?yàn)槊總€周期脫落了一個正渦和一個負(fù)渦。

圖5 計算區(qū)域的劃分Fig.5 The division of computational domain

圖6 單網(wǎng)格下圓柱繞流計算結(jié)果Fig.6 The results of single cylinder with single grid block

表1是本文結(jié)果與文獻(xiàn)研究成果的比較，可以看出本文結(jié)果和文獻(xiàn)結(jié)果符合得較好。其中的Strouhal數(shù)是表征圓柱脫渦情況的重要參數(shù)，定義為

其中：f為圓柱的脫渦頻率，U為自由來流速度。

表1 結(jié)果對比Tab.1 Comparison of results

3 串列雙圓柱相互影響研究

對于靜止單個圓柱的繞流問題，流場特征和圓柱受力情況主要依賴于雷諾數(shù)Re。對于串列雙圓柱繞流問題，除了雷諾數(shù)之外，兩圓柱中心間距比L/D的影響也是關(guān)鍵的（L為前后圓柱中心間距，D為圓柱直徑）。

3.1 計算區(qū)域和網(wǎng)格劃分

本文在耦合網(wǎng)格下對Re=200串列雙圓柱繞流問題進(jìn)行了數(shù)值模擬，計算區(qū)域的劃分如圖7所示。計算區(qū)域分為兩部分，兩個圓柱周圍采用尺寸較小的細(xì)網(wǎng)格，離圓柱較遠(yuǎn)的區(qū)域采用尺寸較大的粗網(wǎng)格。粗網(wǎng)格區(qū)域取粗細(xì)網(wǎng)格尺寸比 m=2。

圖7 計算區(qū)域的劃分Fig.7 The division of computational domain

3.2 計算結(jié)果分析

本文對圓柱中心間距比L/D=1.5、2.0、3.0和4.0的情況進(jìn)行了數(shù)值模擬。不同間距比下的上下游圓柱升力系數(shù)曲線、阻力系數(shù)曲線、流場渦量等值線和流場流線如圖8～11所示。

圖8 L=1.5D情況下雙圓柱計算結(jié)果Fig.8 The results of two cylinders with L=1.5D

圖11 L=4.0D情況下雙圓柱計算結(jié)果Fig.11 The results of two cylinders with L=4.0D

本文將不同間距下上下游圓柱升力系數(shù)、阻力系數(shù)和Strouhal數(shù)進(jìn)行了匯總，并和文獻(xiàn)結(jié)果進(jìn)行了對比，如表2所示。這三個參數(shù)是表征圓柱受力情況的重要參數(shù)，可以作為判斷數(shù)值模擬結(jié)果可靠性的依據(jù)?？梢钥闯觯疚挠嬎憬Y(jié)果與文獻(xiàn)計算值基本相符。

表2 結(jié)果對比Tab.2 Comparison of results

綜合數(shù)值模擬結(jié)果可以看出臨界間距的存在（3.0＜L/D＜4.0）。當(dāng)間距小于該臨界值時，上游圓柱的分離剪切層附著在下游圓柱上，只有下游圓柱有脫渦現(xiàn)象。此時，下游圓柱的阻力系數(shù)較小且小于零。當(dāng)兩圓柱間距超過該臨界值時，上下游圓柱都產(chǎn)生明顯的渦脫落，且升力系數(shù)最大幅值都變大，平均阻力系數(shù)也突然變大，特別是對下游圓柱而言，由一個負(fù)值變到一個較大的正值。

具體來看，當(dāng)L/D=1.5、2.0（小于臨界值）時，下游圓柱的阻力系數(shù)為負(fù)值，即下游圓柱并沒有受到與來流方向相同的阻力，而是受到了與來流方向相反的推力。L/D=1.5、2.0時，兩柱的升力系數(shù)幅值都比較小，尤其上游圓柱尾流受到限制，其升力系數(shù)幾乎為零。當(dāng)間距比增大到4.0（大于臨界值）時，兩圓柱阻力系數(shù)和升力系數(shù)都突然增大。上游柱的升力系數(shù)振幅突然增大至0.8左右，下游圓柱的升力系數(shù)振幅增至更大的值。下游圓柱的阻力系數(shù)由負(fù)值變?yōu)檩^大的正值，但其依然比上游柱的阻力系數(shù)小，每個圓柱的阻力系數(shù)也都小于單柱繞流時的阻力系數(shù)。由升阻力系數(shù)變化曲線還可以看出，升力振動正負(fù)幅值相等，即圓柱所受升力合力為零。另外，上、下游圓柱的渦脫落頻率保持相等，但都較單圓柱繞流小。從Strouhal數(shù)的變化趨勢可以看出，當(dāng)間距比跨越臨界間距值時Strouhal數(shù)會突然增大，而隨著兩圓柱間距的繼續(xù)增大，Strouhal數(shù)接近于單一圓柱繞流情況。

4 結(jié) 論

本文基于格子Boltzmann方法，結(jié)合多塊網(wǎng)格耦合算法和曲邊界處理算法，對雷諾數(shù)Re=200條件下的靜止串列雙圓柱繞流進(jìn)行了數(shù)值模擬。數(shù)值模擬得到了和已有研究成果一致的結(jié)果，驗(yàn)證了臨界間距的存在。當(dāng)兩圓柱間距小于該臨界值時，上游圓柱沒有明顯脫渦現(xiàn)象，下游圓柱出現(xiàn)周期性脫渦；當(dāng)間距大于該臨界值時，上游圓柱開始出現(xiàn)周期性脫渦，而下游圓柱所受的升、阻力系數(shù)明顯提高，并且下游圓柱阻力方向發(fā)生了突變。不同間距下兩圓柱受力參數(shù)和脫渦情況與已有研究成果符合得較好。

[1]Chen S,Doolen G D.Lattice Boltzmann method for fluid flows[J].Annual Review of Fluid Mechanics,1983,30:329-364.

[2]韓文驥,顏 開,褚學(xué)森,張 珂.LBM方法數(shù)值模擬飛濺流動的影響因素研究[J].船舶力學(xué),2015,19(4):356-362.Han Wenji,Yan Kai,Chu Xueseng,Zhang Ke.Influence factors to numerical simulation of splash flow based on LBM[J].Journal of Ship Mechanics,2015,19(4):356-362.

[3]黃橋高,潘 光.疏水表面流體流動特性的格子Boltzmann方法模擬[J].船舶力學(xué),2016,20(10):1211-1218.Huang Qiaogao,Pan Guang.Lattice Boltzmann simulation of liquid flow characteristics of hydrophobic surfaces[J].Journal of Ship Mechanics,2016,20(10):1211-1218.

[4]Zdravkovich M M.Review of flow interference between two circular cylinders in various arrangements[J].ASME Journal of Fluids Engineering,1977,99:618-633.

[5]Zdrvakoeieh M M.The effets of interference between circular cylinders in a cross flow[J].Journal of Fluids and Structures,1987,l:239-261.

[6]梁亮文.低雷諾數(shù)下圓柱橫向受迫振蕩和渦激運(yùn)動的數(shù)值分析[D].上海:上海交通大學(xué),2009.Liang Liangwen.Numerical analysis of forced oscillation and vortex-induced motion of circular cylinder in cross flow with low Reynolds number[D].Shanghai:Shanghai Jiao Tong University,2009.

[7]Yu D,Mei R,Shyy W.A multi-block lattice Boltzmann method for viscous fluid flows[J].International Journal for Numerical Methods in Fluids,2002,39:99-120.

[8]Guo Z L,Zheng C G,Shi B C.An extrapolation method for boundary conditions in lattice Boltzmann method[J].Physics of Fluids,2002,14:2007-2010.

[9]龔 帥,郭照立.流向振蕩圓柱繞流的格子Boltzmann方法模擬[J].力學(xué)學(xué)報,2011,43(1):11-17.Gong Shuai,Guo Zhaoli.Lattice Boltzmann simulation of the flow around a circular cylinder oscillating streamwisely[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(1):11-17.

[10]Chan C T,Anastasiou K.Solution of incompressible flows with or without a free surface using the finite volume method on unstructured triangular meshes[J].International Journal for Numerical Methods in Fluids,1995,29:35-57.

[11]Lecointe Y,Piquet J.On the use of several compact methods for the study of unsteady incompressible viscous flow round a circular cylinder[J].Computers&Fluids,1984,12(4):255-280.

[12]Rosenfeld M,Kwak D,Vinokur M.A solution method for the unsteady and in compressible Navier-Stokes equations in generalized coordinate systems[R].AIAA Paper 88-0718,1988.

[13]劉 松,符 松.串列雙圓柱繞流問題的數(shù)值模擬[J].計算力學(xué)學(xué)報,2000,17(3):260-266.Liu Song,Fu Song.Numerical simulation of flow past two cylinders in tandem arrangement[J].Chinese Journal of Computational Mechanics,2000,17(3):260-266.