姚熊亮 朱永凯 张阿漫 杨树涛

哈尔滨工程大学船舶工程学院，黑龙江哈尔滨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