格子Boltzmann方法在串列雙圓柱繞流數(shù)值模擬中的應(yīng)用研究

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

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

0 引 言

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

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

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

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

1 數(shù)值方法

1.1 格子Boltzmann方法

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

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

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

（1）碰撞過(guò)程：

（2）遷移過(guò)程：

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

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

本文采用經(jīng)典二維九速度模型（D2Q9），將粒子速度離散為9個(gè)方向，即α的值有9個(gè)，如圖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)格來(lái)劃分流場(chǎng)，網(wǎng)格算法相對(duì)更容易實(shí)現(xiàn)。然而，在處理流固耦合問(wèn)題時(shí)，固體壁面周?chē)鷮?duì)網(wǎng)格的密度要求較高，而遠(yuǎn)離壁面處這種要求相對(duì)較低，均勻網(wǎng)格顯然會(huì)造成計(jì)算量的浪費(fèi)。因此，研究多塊網(wǎng)格的耦合算法是必要的。本文在保證質(zhì)量連續(xù)和應(yīng)力連續(xù)前提下，采用了多塊網(wǎng)格耦合的方法來(lái)處理流固耦合中的網(wǎng)格問(wèn)題。

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

前面已經(jīng)提到，運(yùn)動(dòng)粘度和無(wú)量綱松弛時(shí)間應(yīng)滿(mǎn)足（8）式。根據(jù)（8）式，為保證在交界處運(yùn)動(dòng)粘度的連續(xù)，不同網(wǎng)格塊的無(wú)量綱松弛時(shí)間應(yīng)該滿(mǎn)足（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)格交界面上必須滿(mǎn)足速度和密度的連續(xù)，因此

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 單網(wǎng)格下圓柱繞流計(jì)算結(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為自由來(lái)流速度。

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

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

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

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

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

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

3.2 計(jì)算結(jié)果分析

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

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

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

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

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

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

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

4 結(jié) 論

本文基于格子Boltzmann方法，結(jié)合多塊網(wǎng)格耦合算法和曲邊界處理算法，對(duì)雷諾數(shù)Re=200條件下的靜止串列雙圓柱繞流進(jìn)行了數(shù)值模擬。數(shù)值模擬得到了和已有研究成果一致的結(jié)果，驗(yàn)證了臨界間距的存在。當(dāng)兩圓柱間距小于該臨界值時(shí)，上游圓柱沒(méi)有明顯脫渦現(xiàn)象，下游圓柱出現(xiàn)周期性脫渦；當(dāng)間距大于該臨界值時(shí)，上游圓柱開(kāi)始出現(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]韓文驥,顏 開(kāi),褚學(xué)森,張 珂.LBM方法數(shù)值模擬飛濺流動(dòng)的影響因素研究[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]黃橋高,潘 光.疏水表面流體流動(dòng)特性的格子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)動(dòng)的數(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é)報(bào),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]劉 松,符 松.串列雙圓柱繞流問(wèn)題的數(shù)值模擬[J].計(jì)算力學(xué)學(xué)報(bào),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.