Lattice－Boltzmann兩種非均勻網格算法及其對突擴流的模擬

姚熊亮 朱永凱 張阿漫 楊樹濤

哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001

Lattice－Boltzmann兩種非均勻網格算法及其對突擴流的模擬

姚熊亮 朱永凱 張阿漫 楊樹濤

哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001

給出了Lattice－Boltzmann方法兩種非均勻網格算法（即區域分裂方法和坐標變換方法）及計算步驟。對典型的突擴流問題進行了模擬分析，將所得結果與均勻網格進行比較分析。并對兩種非均勻網格算法對流場模擬時進行了比較分析。區域分裂方法的優點在于區域劃分過程靈活，對結構形狀的要求低。坐標變換方法的優點在于只要對區域建立起合適的曲線坐標，計算過程簡單，更節省時間。

非均勻網格；Lattice-Boltzmann方法；突擴流；算法

1 引言

Lattice-Boltzmann方法（LBM）［1］是20世紀80年代中期提出研究流體流動的一種新方法，其基本理論來源于Lattice氣動機。與以往宏觀連續方程的離散化為基礎的傳統數值方法不同，Lattice-Boltzmann方法不是對宏觀的連續方程離散化，而是基于微觀的動力學模型，通過簡單的眾多粒子的微觀上的行為給出宏觀上的動力學方程。在應用Lattice-Boltzmann方法時由于受到網格劃分的限制，非均勻網格［2］的Lattice-Boltzmann方法越來越受到人們的重視。目前已經有幾種非均勻的Lattice-Boltzmann模型，本文將對下面兩種非均勻網格方法進行介紹。

第一種方法是將流場劃分為幾個子區域，在每個子區域使用均勻網格的Lattice-Boltzmann模型，在每個子區域上使用均勻網格的計算方法對Lattice-Boltzmann模型進行求解，在子區域相連的地方采用嵌套邊界的方法進行處理，使整個流場相互聯系起來，從而達到對流場進行模擬。

第二種方法是通過坐標變換的方式達到對區域的非均勻劃分［3］。計算步驟為：1）對流場區域進行計算網格和參考網格兩種劃分；2）設定參考網格上各種物理量的初始值，對應進行單位時間的傳遞，得到下一時刻參考網格上各物理量的值；3）利用插值的方法求解計算網格上各點的物理量；4）繼續在參考網格進行傳遞，得到下一時刻各點上的物理量，插值得到計算網格下的各物理量。

2 Lattice-Boltzmann模型和邊界處理

2.1 二維九點正方形網格模型

本文中采用二維九點正方形網格模型［4］，該模型是將物理流場離散為正方形網格（圖1）。每個網格節點同周圍的8個節點構成節點的9個可能的運動方向，離散的速度張量eα可以表示為：

圖1 二維九點模型速度張量圖

2.2 Lattice－Boltzmann方程Lattice－Boltzmann方法中的 Lattice－Boltmann方程［5］可以寫為如下一般形式：

式中，fα（x，t）和fα（eq）（x，t）分別為與第α個方向上的速度eα相對應的粒子速度分布函數和平衡分布函數，δt為時間步長，τ為無量綱松弛時間。

Lattice－Boltzmann方程中的平衡分布函數fα（eq）（x，t）可由Maxwell－Boltzmann分布函數對速度u的二階Taylor展開式［6］獲得，可以表示為：

其中，ωα為粒子分布函數的加權系數，

根據質量守恒與動量守恒定律有：

式中，ρ表示密度，u表示流速。

2.3 Lattice－Boltzmann方法的邊界處理

邊界條件的處理方法在Lattice－Boltzmann方法中起著重要的作用，對Lattice－Boltzmann模型的精度和穩定性都有很大影響。常用邊界處理［7］方法有：反彈格式、差值法、水動力條件方法等幾種方法，在本文中只對反彈格式的方法進行介紹，具體過程如下。

在該方法中，當粒子到達壁面格點后，粒子被沿原路返回到流體內部，且運動方向和入射方向相反。反彈格式非常容易實現，并且適用于處理幾何形狀非常復雜的系統，如多孔介質內的多相流動問題。

如圖2所示，以下壁面為例，對于固壁上的節點A，當經過流過程后，由A鄰近的B，C，D，E，F點可以分別確定格點A的5個方向的分布函數f8，f7，f3，f1，f4。

圖2 壁面反彈

根據壁面反彈的思想可知：

f2＝f4，f5＝f7，f6＝f8

則點A所有方向的分布函數可以確定。標準的壁面反彈方法同時還要求壁面節點不進行碰撞，以保證無滑移條件。然而這種處理方法得到的精度有限。如果讓固壁節點參與碰撞，可以看到反彈條件并不能保證壁面切向速度為零，因為f1＝f3。為了滿足無滑移條件，可以對部分分布函數fi進行如下重分配。

1）利用反彈邊界初步確定所有的分布函數，然后計算出ρ，u；

2）對具有切向分量的流入流體內部的分布函數fi做如下分配。

其中，a為常數，對于d2q9模型，a＝1／2，以下壁面為例，可得

代入方程中，可知u＝0，v＝0，滿足無滑移條件。上述過程就是壁面反彈邊界處理邊界的過程。

3 非均勻Lattice－Boltzmann模型

近年來，Lattice－Boltzmann方法已越來越引起人們的重視。然而，在構造Lattice－Boltzmann模型時，卻會常常受到規則網格的限制。因此非規則網格的Lattice－Boltzmann方法顯得越來越重要。非規則網格得到了更廣泛的應用，對不需過細劃分的地方采用大網格，可以節省計算空間和時間。

3.1 區域分裂的非均勻Lattice－Boltzmann方法

首先我們根據實際情況對流場［8］進行區域劃分，在不同的區域用不同的網格尺度求解。在不同區域內部我們可以用均勻Lattice求解，相鄰區域的粗網格區域設為Ωc，相應粗網格上的網格間距Δc，松弛時間設為τc；細網格區域設為Ωx，相應細網格上的網格間距Δx，松弛時間設為τx；粗細網格相交的部分稱為虛擬邊界Γcx。這里

從式（10）可以看出，在粗網格上傳遞一次時，細網格相應傳遞n次，這就需要對粗細網格相接的地方進行處理，這里采用插值的方法計算。考慮到Lattice－Boltzmann方法是二階精度格式，這里使用二階的插值方法，具體過程如下：

1）對于對應點在粗網格上節點的插值（如圖3上b1點的插值），則有

2）對于對應點不在粗網格上節點的插值（如圖3上b2點的插值），則有

其中，k為a1點與a2點之間的距離和a2點與a3點之間距離的比值。

圖3 粗細網格相接處示意圖

3.2 坐標變換的非均勻 Lattice－Boltzmann方法

首先，流場區域建立直角坐標系（x，y），并對其進行網格劃分。區域建立合理的曲線坐標系（ε，η），利用曲線坐標系對流場區域進行非均勻網格劃分。接下來利用均勻網格給定流場初始分布函數fα，進行單位時間的網格傳遞，得到下一時刻均勻網格上的分布函數。利用插值的方法，將均勻網格上的分布函數轉換為曲線坐標系非均勻網格上的分布函數。由于Lattice－Boltzmann方法是二階精度格式，這里使用插值方法也要達到二階精度以上，下面介紹一種插值方法。

1）如圖4所示，當非均勻網格的節點a1在均勻網格坐標系上任意點時，a1在均勻網格坐標系中對應的坐標為 （x1，y1），b1－b4的坐標分別為（xb1，yb1）、（xb2，yb2）、（xb3，yb3）、（xb4，yb4），則有a1上的分布函數為：

其中，

2）當非均勻網格的節點a2在均勻網格的網格線上時，a2在均勻網格中對應的坐標為（x2，y2），c1、c4的坐標分別為（xc1，yc1）、（xc2，yc2），則有a2上的分布函數為：

其中，

3）當非均勻網格的節點a3在均勻網格的網格節點上時，a3的分布函數的值就等于該點上的分布函數的值。

圖4 節點插值示意圖

對于其它的插值方法，只要保證達到二階精度就可以適當采用，對其討論在此略去。

4 數值模擬及結果分析

本文將利用上述非均勻網格模型對突擴流進行數值模擬。流場區域如圖5所示，將其劃分為流場前段的區域1和流場后段的區域2兩部分。初始時，流場左側入口處節點速度為0.2，其余節點速度為0，密度設定為1。

圖5 流場示意圖

4.1 區域分裂的非均勻 Lattice－Boltzmann方法模型

對于圖5所示的模型進行網格劃分，區域1為粗網格，劃分為20×17；區域2為細網格，為81×81。進行數值計算時，雷諾數設定為300，粗細網格的Lattice間距為n＝Δc／Δx＝2。得出區域分裂的非均勻網格方法的流場向量圖如圖6所示，對應向量圖相應畫出漩渦渦心處速度等值線，如圖8所示。

用該方法所得的速度等值線與均勻網格方法（該模擬過程較為簡單，在此略去）所得的等值線相吻合，并用該方法對漩渦渦心位置（由于該結構是對稱結構，故只取一個渦心的位置比較）與均勻網格進行比較如表1所示。通過比較可以看出，這種方法所得的結果與均勻網格方式所得的結果相吻合，應用該方法對流場進行模擬是可行的。

4.2 坐標變換的非均勻 Lattice－Boltzmann方法模型

該方法需要用到計算網格和參考網格兩種坐標，參考網格進行均勻網格劃分，流場區域1為40×30，區域2為72×72；計算網格［9］區域1進行如下劃分為40×30的網格，網格節點的位置為：

區域2進行如下劃分為72×72的網格，區域2網格節點的位置為：

進行數值計算時，雷諾數設定為300，得出坐標變換的非均勻網格方法的向量圖如圖7所示，對應向量圖相應畫出漩渦渦心處速度等值線如圖9所示。

圖6 區域分裂的非均勻網格方法的流場速度向量圖

圖7 坐標變換的非均勻網格方法的流場向量圖

該方法所得的速度等值線與均勻網格方法所得的等值線相吻合，并用該方法漩渦渦心位置（同上，只取一個渦心的位置進行比較）與均勻網格進行比較如表2所示。通過比較可以看出，這種方法所得的結果與均勻網格方式所得的結果相吻合，應用該方法對流場進行模擬也是可行的。

圖8 區域分裂方法漩渦渦心處速度等值線圖

圖9 坐標變換方法漩渦渦心處速度等值線圖

表1 區域分裂方法與均勻網格方法的比較

4.3 區域分裂與坐標變換兩種方法的比較

通過上面的數值模擬可以看出，上面兩種非均勻網格方法成功地解決了對區域某些地方需要使用較細的網格離散，而遠離該區域的地方采用較粗的網格離散。通過應用兩種方法的數值計算過程和所得結果的比較可以得出如下結論。

區域分裂方法的優點在于可以對區域任意的部分進行粗細劃分，與坐標變換方法比較，該方法區域劃分過程靈活，對流場結構形狀的要求更低；坐標變換方法的優點在于只要選取合適的曲線坐標，該方法的計算過程簡單，更節省時間，而該方法的缺點在于要為流場區域選取合適的曲線坐標，如果不能將流場劃分為合適的曲線坐標，那么這種方法將不能應用。

5 結束語

本文詳細介紹了Lattice－Boltzmann方法的兩種非均勻網格方法，并應用兩種方法對突擴流進行模擬，與均勻網格進行了比較。通過兩種方法對流場的模擬過程的比較，得出了兩種方法對流場模擬時各自的優缺點。文中只對結構簡單的突擴流進行了模擬，接下來將嘗試把兩種方法結合起來對流場進行非均勻網格模擬。

［1］CHEN S Y，DOOLEN G D.Lattice Boltzmann methods for fluid flows［J］.Annu Rev Fluid Mech，1998，30：329.

［2］WOLFRAM S.Cellular Automaton Fluids 1：Basic Theory［J］.J Stat Phys，1986，45（3/4）：471－526.

［3］HE Xiao－yi，LUO Li－shi，DEMBO M.Some Prograss in Lattice Boltzmann method.part I：nonuniform mesh grids［J］.J Comput Phys，1996，129：357－363.

［4］BHATNAGAR P L，GROSS E P，KROOK M.A model for collision processes in gases I：small amplitude processes in charged and neutral one－component systems［J］.Phys Rev，1954，94（3）：511－525.

［5］QIAN Y H.d′BUMIèRES D，LALLEMEND P.Lattice BGK models for Navier－Stokes equation［J］.Europhys Lett，1992，17（6）：479－484.

［6］NOBLE D R，GEORGIADIS J G，BUCKIUS R O.Comparison of accuracy and performance for Lattice Boltzmann and finite difference simulation of steady viscous flow［J］.Int J Numer Meth Fluids，1996，23（1）：1－18.

［7］聶德明，林建忠.Lattice－Boltzmann方法中的邊界條件［J］.計算物理，2004，21（1）：21－26.

［8］郭照立，施保昌，王能超.基于區域分裂的非均勻Lattice Boltzmann方法［J］.計算物理，2001，18（2）：181－184.

［9］程永光.基于插值的Lattice Boltzmann方法非均勻網格算法［J］.武漢水利電力大學學報，2000，33（5）：26-31.

Simulation of 2D Sudden-expansion Flow Based on Two Algorithms of Non-uniform Mesh Grids Using Lattice-Boltzmann Method

Yao Xiong－liang Zhu Yong－kai Zhang A－man Yang Shu－tao
College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China

Two non－uniform mesh grid algorithms （viz the methods based on domain decomposition technique and coordinate transformation technique）for the Lattice-Boltzmann method based on interpolation technique，with detailed calculation process，are presented.The simulation of the two－dimensional sudden－expansion flow was carried out in the current study，and the simulated results were then compared with the results calculated with the uniform mesh grid algorithm.At the same time，the contrast between the two methods were also conducted during the process of simulation.The advantages of the method based on domain decomposition technique are that the flow domain can be divided flexibly and the method can be employed to structures with various shapes，while to the method based on coordinate transformation technique，a suitable curvilinear coordinate system should be established，and the process of calculation is simple，highly efficient in computation.

non－uniform mesh grid；Lattice-Boltzmann method；sudden－expansion flow；algorithm

O35

A

1673－3185（2009）02－15－05

2008-11-11

姚熊亮（1963－），男，教授，博士生導師。研究方向：船舶與海洋工程結構動力學。E-mail：saibei8411＠163.com

朱永凱（1982－），男，碩士研究生。研究方向：船舶與海洋工程結構動力學。E-mail：zhu432122＠163.com