徐 磊 宋安平 劉智翔 張 武

1(上海大學計算機工程與計算科學學院 上海 200444)2(上海海洋大學信息學院 上海 201306)

0 引 言

格子Boltzmann方法LBM是介于流體的微觀分子動力學模型和宏觀連續(xù)模型之間的介觀模型。由于Boltzmann方程自身的運動學特性，以及可以根據(jù)經(jīng)典的Chapman-Enskog展開從LBM得到Navier-Stokes方程，使得LBM比基于連續(xù)介質假設的Navier-Stokes方程包含了更多的物理內(nèi)涵。目前，LBM已經(jīng)應用在了很多流體問題(多項流、湍流、多孔介質以及微尺度流等)中[1-4]。LBM將流體視為比分子大，但在宏觀上無限小的一系列粒子。這些粒子按照一定物理規(guī)律在網(wǎng)格上進行演化計算，通過對反映粒子狀態(tài)的分布函數(shù)進行統(tǒng)計平均求得宏觀物理量。

文獻[5-6]首先分別提出了單松弛時間格子Boltzmann模型SRT-LBM(Single-Relaxation-Time LBM)和格子BGK模型LBGK(Lattice Bhatnagar-Gross-Krook)。它們采用單一松弛時間系數(shù)，計算效率得到極大提高。盡管SRT-LBM已經(jīng)被成功地應用于模擬各種流體問題，但是該模型仍然存在缺點：當無量綱松弛時間過于趨近0.5時，模型會產(chǎn)生數(shù)值不穩(wěn)定導致程序發(fā)散。針對這一問題，D′Humiéres[7]提出了廣義格子Boltzmann模型GLBE(Generalized Lattice Boltzmann Equation)，或稱為多松弛格子Boltzmann模型MRT-LBM(Multiple-Relaxation-Time LBM)。2000年-2002年間，文獻[8-10]對該模型對該模型二維和三維的物理原理、參數(shù)選取、數(shù)值穩(wěn)定性和計算效率等方面進行了詳細的理論分析，表明了MRT-LBM的穩(wěn)定性和精確性都要優(yōu)于SRT-LBM。

相對于SRT-LBM，MRT-LBM穩(wěn)定性更高，但是并未克服松馳時間趨近于0.5時計算不穩(wěn)定導致程序發(fā)散的問題。隨著雷諾數(shù)的增大，MRT-LBM和SRT-LBM雖然都可通過增加網(wǎng)格數(shù)保持計算穩(wěn)定，但是當雷諾數(shù)更大時，利用SRT-LBM和MRT-LBM直接數(shù)值模擬流場全部動態(tài)信息所需要的內(nèi)存和求解時間非常巨大[11]。為了在模擬湍流場時節(jié)省計算耗費并在有限的計算機硬件條件下盡可能完整地給出瞬時流場信息，可在格子Boltzmann模型中引入大渦模擬技術LES(Large Eddy Simulation)[12]。Hou等[13-14]通過非平衡粒子分布函數(shù)的二階矩計算了應變率張量，從而將Smagorinsky渦粘性模型引入SRT-LBM中。而Kxafczyk等[15]則利用傳遞矩陣將粒子分布函數(shù)的二階矩從速度空間傳遞到矩空間計算禍粘性系數(shù)，將Smagorinsky渦粘性模型與D3Q15 MRT-LBM組合。文獻[16-17]通過Chapman-Enskog分析推導應變率張量，分別將Smagroinsky模型引入到D3Q19 MRT-LBM和包含外力項的D3Q19 MRT-LBM中，得到MRT-LBM-LES(MRT-LBM with LES)模型。

如今，計算流體力學研究對象的規(guī)模和復雜程度向更大更深的方向發(fā)展，單個計算節(jié)點不能滿足計算需求，需要借助于高性能計算機滿足計算的需要。大規(guī)模并行計算在一定程度上可以解決計算需要和超級計算性能之間的矛盾。LBM中格點演化主要分為碰撞和遷移兩個過程，碰撞是各個格點獨自同時進行的，只與格點自身相關。遷移時格點上的粒子也只是與離它們最近格點上的粒子進行信息交換，所以LBM非常適合并行計算。諸多研究人員對LBM的并行性能進行了研究，設計了高可擴展高效率的并行算法。Pan等[19]在不同并行計算平臺比較區(qū)域劃分方法的性能，利用D3Q15模型模擬了單向和多項流在多孔介質中的流動情況。Wu等[20]對高雷諾數(shù)二維方腔流，采用域分解方法對D2Q9 MRT-LBM和SRT-LBM進行了比較。Velivelli等[21]利用CPU緩存的優(yōu)勢，將計算區(qū)域分塊循環(huán)計算，加快了計算速度。Schepke等[22]采用塊劃分策略對SRT-LBM進行了并行性能的分析，該方法使得每個進程都分配到整個計算區(qū)域中的一個子計算區(qū)域，各個子計算區(qū)域與相鄰的子計算區(qū)域進行數(shù)據(jù)的傳遞。文獻[23]針對多孔介質內(nèi)流體流動提出了負載均衡策略減少計算時間。本文主要研究了MRT-LBM引入Smagorinsky渦粘性模型后在大規(guī)模并行計算機上的并行性能。

1 MRT-LBM-LES以及邊界條件

1.1 MRT-LBM-LES

MRT-LBM將速度空間的碰撞通過線性變換轉化為矩空間的碰撞，使得作用后得到的量具有物理意義。在矩空間完成碰撞后，再進行逆變換到速度空間，進行遷移的過程。MRT-LBM的碰撞過程是在矩空間中進行的，它的演化方程為：

f(x+ciδt,t+δt)-f(x,t)=-M-1S[m(x,t)-meq(x,t)]

(1)

式中：f(x,t)是粒子分布函數(shù)，M是變換矩陣，m(x,t)=Mf(x,t)是矩空間的粒子分布函數(shù)，meq(x,t)是矩空間的平衡態(tài)分布函數(shù)。S是松弛對角陣，并且S對角線各項滿足0≤Si<2。

式(1)可以分為碰撞和遷移兩個過程:

f(x,t)=f(x,t)-M-1S[m(x,t)-meq(x,t)]

(2)

f(x+ciδt,t+δt)=f(x,t)

(3)

常用的二維LBM模型為D2Q9(如圖1所示)。它的離散速度為:

(4)

圖1 D2Q9

平衡態(tài)分布函數(shù):

(5)

(6)

D2Q9對應的變換矩陣:

(7)

矩空間的平衡態(tài)分布函數(shù)meq為:

meq=ρ(1,-2+3u2,α+βu2,ux,-ux,uy,-uy,

(8)

式中:α和β為自由參數(shù)，ρ為密度，u為速度。松弛矩陣diag(0,se,sε,0,sq,0,sq,sv,sv)。剪切粘性和體粘性系數(shù)分別為：

(9)

(10)

宏觀量為:

ρ=∑fi

(11)

(12)

為了能夠模擬湍流，Hou等[13]將Smagroinsky模型引入到D2Q9中。在引入的Smagorinsky模型中，無量綱的運動粘性系數(shù)v由分子運動粘性v0和渦粘性vt組成:

υ=v0+vt

(13)

式中：vt由應變率張量Sαβ，參數(shù)Cs以及格子步長δx決定,其表達式如下：

vt=(Csδx)2|Sαβ|

(14)

1.2 邊界條件

在處理曲面復雜邊界時，如果使用平直邊界的邊界處理方法，精度較低。為此，Yu等[24]提出了一種曲面邊界的處理方法。如圖2所示，實心圓表示固體點，空心圓表示流場點，曲線表示實際的物理邊界，它將流場分為固體點和流場點，方格點為邊界節(jié)點。

圖2 曲面邊界

這里將流場點記為xf，固體點記為xb，與xf相鄰的流場點記為xff，邊界點記為xw。則流場點xf的分布函數(shù)為：

(15)

(16)

對于平直邊界，本文采用Guo等[25]提出的非平衡態(tài)外推方法。

2 格子類型的判斷

當流場中物體存在復雜邊界時，需要判斷物體附近格子點的類型，通過不同類型的格子點計算出物體邊界處的格子點粒子分布函數(shù)。本文通過射線法判斷格子的類型。從待判斷點的某一個方向引射線，計算和物體邊界的交點個數(shù)，通過交點個數(shù)的奇偶性得到格子點的類型。本文使用圓柱和Rae2822翼型為例判斷格子點的類型，如圖3和圖4所示。該方法可以準確的判斷出圓柱和翼型流場中的格子類型。

圖3 圓的格子類型(包括固體點和邊界點)

圖4 Rae2822翼型的格子類型(包括固體點和邊界點)

3 并行策略

對于二維問題，流場可以沿著一個方向或兩個方向進行劃分。當并行環(huán)境中核數(shù)增多時，沿著一個方向進行劃分并不符合實際問題的需要，所以本文選擇沿著兩個方向對流場進行劃分。在進行并行處理時，每個MPI進程負責處理一個流場的子區(qū)域，每個子區(qū)域用一個二元組(i,j)標記。進程總數(shù)記為n。二元組可以通過下式獲得：

i=mod(n,nx)

(17)

(18)

式中：nx為沿著x方向的劃分數(shù)。每個子區(qū)域的范圍為:

(19)

(20)

式中：nLx是x方向格子點數(shù)，xmin是子區(qū)域的起始格點位置，xmax是子區(qū)域的結束格點位置。子區(qū)域在y方向的范圍可以類似計算得到。

圖5 進程(i,j)傳遞部分數(shù)據(jù)給進程(i+1,j)

圖6 流場中各個子區(qū)域信息傳遞

4 數(shù)值實驗

4.1 數(shù)值驗證

為了驗證本文并行算法的正確性，我們使用圓柱和翼型在流場中的情況分別進行了測試。流場的速度為U=0.115 5。圖7中，(a)和(b)是圓柱繞流雷諾數(shù)Re=20和40時的流線圖，從圖中可以看出，在圓柱后產(chǎn)生了一對對稱的渦，隨著雷諾數(shù)的增大，渦也在逐漸變大。(c)和(d)是雷諾數(shù)Re=100和150時的流線圖，圓柱后方兩側的渦隨著雷諾數(shù)的增大開始出現(xiàn)震蕩，當超過臨界值時，圓柱后的渦逐漸分離，形成周期性交替脫落、旋轉方向相反的非對稱渦[26]。

圖7 圓柱繞流流線圖

圖8中的(a)和(b)是翼型繞流雷諾數(shù)Re=100 000時的流線圖，(a)的攻角為0°，(b)的攻角為3°。可以看出，實驗結果與文獻[27]一致。

圖8 翼型繞流渦量圖

4.2 效率驗證

本節(jié)基于128核集群進行并行性能的測試，該集群包含1個登錄節(jié)點、8個計算節(jié)點，每個節(jié)點有16個核，配置了兩塊2×Intel Xeon E5-2660 CPU。每個節(jié)點的內(nèi)存為128 GB。節(jié)點之間通過56 GB Infiniband線連接。

為了統(tǒng)計程序運行時間，進行了500步迭代。采用了1 000×1 000和1 500×1 500兩種格子規(guī)模進行了加速比和效率的測試，比較了兩種不同規(guī)模下的加速比和效率。

分別在核數(shù)為8、16、32、64以及128核上進行了時間的統(tǒng)計。圖9為兩種不同格子規(guī)模下加速比的比較，圓形為格子規(guī)模1 500×1 500的加速比，正方形為格子規(guī)模1 000×1 000的加速比，黑色直線為理想加速比。可以看出，兩種格子規(guī)模的加速比都呈直線上升，格子規(guī)模較大的加速比要略高。由于進程之間需要進行邊界上格子信息的通信，加速比要低于理想加速比。因為進程間的通信量只發(fā)生在相鄰進程間，并沒有大規(guī)模的主從式通信，所以加速比基本為直線。

圖9 加速比

圖10為兩種不同規(guī)模下效率的比較，圓形為格子規(guī)模1 500×1 500的效率，正方形為格子規(guī)模1 000×1 000的效率?？梢钥闯?，兩種格子規(guī)模的效率都在緩慢下降，隨著核數(shù)的增加，格子規(guī)模較大的效率下降趨勢更為緩慢，當核數(shù)達到128核時，計算效率可以達到86.5%。從圖9和圖10可以看出，我們所實現(xiàn)的并行算法有良好的可擴展性。

圖10 效率

5 結 語

本文提出了一種高可擴展的格子Boltzmann方法。該方法可以有效判斷格子類型，模擬高雷諾數(shù)下流場中的湍流，采用雙向的區(qū)域分解方法，并設計了詳細的劃分策略以及數(shù)據(jù)交換策略，大大降低了通信量，減少了通信時間。數(shù)值結果表明，該算法具有良好的可擴展性，效率可以在128個核心上達到86.5%。