格子Boltzmann方程模拟高雷诺数三维方腔流

康海贵，李承功

（大连理工大学海岸和近海工程国家重点实验室，大连116023）

格子Boltzmann方程模拟高雷诺数三维方腔流

康海贵，李承功

（大连理工大学海岸和近海工程国家重点实验室，大连116023）

为分析方腔流内部流场的特性和验证格子Boltzmann方法模拟湍流的能力，应用标准Smagorinsky涡粘性模型与多松弛时间格子Boltzmann方程（Multiple Relaxation Time Lattice Boltzmann Equation，MRT-LBE）组合对高雷诺数（Re=10 000）三维方腔流进行数值研究，计算了时间平均量，如速度，均方根速度、雷诺应力以及中心断面（y=W/2）处的流线等高线。模拟结果与已有实验和数值模型结果比较可知，MRT-LBE能够精准地计算剪切驱动方腔内流场的变化。另外，将基于图形处理（graphic processor unit，GPU）的计算统一设备架构（Compute Unified Device Architecture，CUDA）并行技术引入到基于MRT-LBE的Smagorinsky模型以提高计算效率，计算效率提高达200倍。

多松弛时间格子Boltzmann方程；标准Smagorinsky涡粘性模型；高雷诺数三维方腔流；GPU；CUDA

格子Boltzmann方法［1］（Lattice Boltzmann Method，LBM）是克服了格子自动机方法（Lattice Gas Cellular Automata，LGCA）一系列缺点并继承其优点发展而来的一种新的数值方法。与传统的求解Navier-Stokes动量方程组的数值方法相比，LBM从介观的角度通过分布函数演化流体粒子的碰撞和迁移过程，从而描述流体的宏观运动，该方法具有边界处理简单，易于编程和适合并行的优点，因此引起了各领域学者的极大关注，并已经被成功应用于模拟众多复杂流体流动现象［2-3］。

近年来，LBM还被用来解决海洋和海岸工程领域的相关问题。例如Zhou等（2000）提出了浅水方程的LBM模型及包含湍流模型的形式，并利用该模型模拟了圆柱绕流，渠流和潮流等水动力问题［4］。张金凤和张庆河（2011）利用LBM对河口海岸细粒黏性泥沙进行数值研究，分析了三维分形絮团的沉降特性，对球形颗粒在静水中沉降引起的紊动流场进行了数值模拟［5］。在开发利用海洋能发电上，王树杰等（2011）应用LBM对水轮机的引导槽装置进行了水动力优化［6］。Janβen等（2012）应用基于完全非线性浅水波方程的LBM数值模拟了海洋长波的传播和爬坡过程［7］。最近，LBM与VOF方法相结合直接求解Navier-Stokes方程已经被用来计算自由表面流体流动问题，如破碎波［8］等。

方腔流通常是指顶部平板以恒定速度驱动规则区域内封闭的不可压流体（例如水）的流动，由于在方腔流的流动中可以观察到几乎所有可能发生在不可压流体中的流动现象（如涡，二次流，Taylor-Görtler-like（TGL）涡，流体分岔，不稳定，瞬变以及湍流等），因此对方腔流已经进行了广泛的实验和数值研究［9-21］。最初的方腔流研究始于Burggraf等（1966），他们应用解析和数值的方法计算了封闭区域内由具有恒定速度的顶部平板驱动的稳定流场［9］，随后Pan等（1967）在实验中应用摄影技术观测了矩形方腔内的流场结构，但是他们忽略了三维作用对流场的影响［10］。Koseff等（1984）最早应用可视化技术对雷诺数Re≤10 000的方腔流进行了深入的实验研究，详细讨论了腔内流场的三维特征，角涡，TGL涡，流场的不稳定状态，边墙摩擦对腔内流场的影响等［11-13］。Prasad和Koseff（1989）根据Koseff和Street的实验结果分析了表征湍流脉动强度的二阶统计量（均方根速度，雷诺应力），进一步考察了边墙，雷诺数和TGL涡对腔内流场的影响［14］。2000年Shankar和Deshpande对方腔流给出了详细的综述概括［15］。最近，Guermond等（2002）对三维方腔流（宽长高比为1：1：2，Re=1 000）的初始阶段流场进行了实验和数值研究［16］。Leriche等（2000）［17］和Bouffanais等（2007）［18］对Re=12 000的三维方腔流分别进行了直接数值模拟和大涡模拟研究。在应用LBM数值计算上，d′Humières等（2002）应用MRT-LBE对Re≤4 000三维对角驱动方腔流进行模拟证明了MRT模型的的稳定性要由于LBGK模型［19］。De等（2009）应用MRT-LBE对Re≤1 000的三维方腔流进行数值模拟，以研究在长宽比A=1条件下不同深宽比对腔内三维流场结构的影响和评估MRT模型计算三维流场的能力［20］。Premnath等（2009）应用广义格子Boltzmann方程模拟了Re=12 000的三维方腔流［21］。

本文首先对三维方腔流初始流场进行求解以验证标准的Smagorinsky涡粘性模型［22］与MRT-LBE［19］组合的MRT-SMAG模型的有效性，其次为了评估LBM模拟湍流的能力为今后建立数值波浪水槽提供技术基础，利用该模型对高雷诺数三维方腔流进行数值实验，并与物理模型实验结果和其他数值结果进行比较验证，并考察了宽高比对腔内流场的影响。为了提高LBM求解三维流场的计算效率和针对LBM易于并行的特性，采用了最新的GPU-CUDA并行计算［23］以缩减模型中粒子的碰撞和迁移时间。

1 控制方程

1.1 MRT-LBE

对于D维空间具有Q个粒子运动方向的多松弛时间格子Boltzmann模型的演化方程为［19］

本文采用D3Q19格子模型［19］（图1），因此离散的粒子速度矢量为

1.2 基于Smagorinsky涡粘性模型的MRT-LBE模型

大涡模拟（Large Eddy Simulation，LES）［24］的主要原理是直接求解大尺度（或称可解尺度）的湍流运动，而小尺度（或称不可解尺度）的紊动作用则利用相应的涡粘性模型模拟出来以减少计算耗费。本文主要采用的Smagorinsky（1963）［22］所提出的涡粘性模型求解残余应力，该模型的主要形式为

将标准的Smagorinsky涡粘性模型与MRT-LBE组合的关键在于建立粒子分布函数与应变率张量之间的关系，即通过非平衡分布函数的二阶矩计算应变率张量［25］

1.3 边界条件

本文采用动量改进的半步反弹格式处理无滑移的墙边界［19］

2 数值模拟和结果讨论

2.1 数值验证

三维不可压方腔流流动的计算区域如图2所示，在长宽高分别为L，W和H的规则区域内充满了不可压缩流体（水），区域顶部平板以恒定的速度Up从沿x轴方向左向右移动，从而驱动区域内流体的流动。区域顶部为速度边界，其余为静止的墙边界，因此采用动量改进的半步反弹格式进行处理，具体形式为（14）式。应用LBM进行模拟计算通常需要将实际的物理量转化为无量纲的格子量［27］，下面定义r，t，u和υ分别表示长度，时间，速度和运动粘度系数，相应的大写字母L，T，U表示对应的特征量，角标为p的量是有量纲的物理量，而角标为lb的量表示无量纲的格子量。

本节首先对三维方腔流初始阶段（即从0时刻开始到第12周期时间内，并采用进行归一化≤12）流场进行模拟以验证D3Q19 MRT-SMAG模型的有效性，其中方腔的长宽高分别为=6.2 cm，=12.4 cm，=6.2 cm，顶部驱动速度为1.8 cm/s，雷诺数Re=1 000，根据雷诺数的定义可知水的运动粘度系数=1.116×10-6m2/s。对计算区域采用127×255×127的规则格子网格进行离散，模型中顶部驱动速度=0.1以确保（马赫数Ma=＜0.3），因此根据单位转换公式可以推导出其他格子量。Re=1 000的三维方腔流流动表现为层流状态，因此Smagorinsky常数取零，碰撞对角矩阵中的松弛率分别取为

由图3可知，方腔对称截面（y=W/2）上分别沿x=L/2和z=H/2的不同时刻（t=4，6，8，10，12）的瞬时水平和垂直速度的计算结果与实验和其他数值方法的结果符合较好，这表明了MRT-SMAG模型可以对流场的初始阶段进行求解。但是图3所示的数值结果（本文模型的和文献［16］中基于求解Navier-Stokes方程的）与实验结果都略有不同，这可能是由实验中顶部平板速度的微小变化所导致的。

2.2 Re=10 000三维方腔流数值结果及结果讨论

方腔内的流动形式主要取决于方腔几何尺寸（如宽高比）和雷诺数（在方腔尺寸和封闭流体粘度系数不变的情况下反映顶部驱动速度的大小）。通常认为雷诺数为2 000时，腔内流场为稳定的层流，当雷诺数在2 000～3 000时，方腔底部的角涡表现出不稳定状态，随着雷诺数的增大，腔内流体的不稳定状态更加剧烈，当雷诺数在6 000～8 000时，腔内流场首次表现出湍流特征（如扩散现象），当雷诺数大于等于10 000时，方腔底部的角涡完全发展为湍流［13-14］。因此，为了评估MRT-SMAG模型模拟湍流的能力为今后建立基于LBM的数值波浪水槽提供技术基础，本文利用D3Q19 MRT-SMAG模型对Re=10 000的三维方腔流进行数值计算，对比了方腔宽高比K=W/H=1和3两种情况下的时间平均速度以分析边墙对内部流场的影响，给出了K=1方腔流的时间平均均方根速度以及雷诺应力曲线，并与实验和其他数值结果进行比较。

腔的长和高分别为L=15 cm，H=15 cm，根据宽高比可知方腔宽分别取为W=15 cm和45 cm，腔内封闭的流体为不可压缩的水，其运动粘度系数取，由雷诺数定义可以推导出方腔顶部驱动速度≈0.067 m/s。对计算区域分别采用127×127×127（K=W/H=1）和100×300×100（K=W/H=3）的规则格子网格进行离散，并定义格子系统中顶部驱动速度0.1，再将实际物理量转换为无量纲的格子量。Smagorinsky常数=0.12，碰撞对角矩阵中的松弛率分别取为。为了获得精确的时间平均结果，首先选取作为初始计算阶段以确保腔内流场完全发展，再对150的瞬时结果进行时间平均以获得稳定的时均结果。选取不同的初始计算时间（最大为）和时间平均周期（最多平均1 000的瞬时结果）以分析其对时间平均量的影响，分析结果表明选取更长的初始计算时间和平均时间对结果的影响可以忽略不计。

图4给出了利用D3Q19 MRT-SMAG模型模拟Re=10 000三维方腔流得到的时间平均水平速度＜u＞和垂直速度＜w＞，MRT-SMAG模型结果与实验［13-14］和直接数值模拟结果［17］符合较好。另外，无论是实验和直接数值模拟结果还是MRT-SMAG模型计算结果都表明了中心对称断面处宽高比K=3的时间平均速度要略大于K=1的结果，这是由边墙对腔内流体的摩擦作用造成，即宽高比越大，边墙摩擦对中心对称断面的影响逐越低，对湍流脉动的掺混能力越弱。

本文还计算了腔内流场的二阶统计量以分析湍流脉动强度分布，这些量均以和行归一化。图5-a表示在中心对称断面（y=W/2）上与顶部驱动方向一致的时间平均均方根速度，即在x= L/2处沿垂向方向（z轴）的分布，而图5-b表示在中心对称断面（y= W/2）上与顶部驱动方向垂直的时间平均均方根速度，即在z=H/2处沿水平方向（x轴）的分布。另外，图6-a和6-b分别给出了在中心对称断面（y=W/2）上时间平均雷诺应力在x=L/2处沿垂直方向（z轴）的分布和在z=H/2处沿水平方向（x轴）的分布。这些结果表明了腔内湍流是沿墙边界产生，下游边界（x=L）附近的湍流脉动要比上游边界（x=0）附近的强烈，而底部强边界（z=0）附近的湍流脉动强度是最大的。这与Leriche等（2000）对三维方腔流的讨论是一致的。因此腔内湍流强度的分布可由腔内动量传递过程解释，即速度边界首先将动量传递到顶部边界附近的粘性层当中使得顶部边界的下游处压力增大，增大的压力抑制了粘性层的运动并驱动流体沿着下游墙边界（x=L）向下运动，类似于沿平直壁面的射流，该射流沿着下游墙边界的中心线分解为两个近似椭圆的自由射流，随后以一个非常小的角度冲击到底部墙边界导致湍流脉动的产生［17］。通过与已有数据进行对比可知，应用MRT-SMAG模型计算的结果与实验和直接数值模拟以及大涡模拟的结果在峰值上略有偏差但是总体趋势符合较好，表明了该模型有能力对湍流场进行求解，其中模拟结果与已有数据的偏差可能是由顶部速度边界的不同造成的，由于在方腔流的实验中运动平板与腔内流体接触会带走少量流体，造成顶部边界的速度分布无法确定，但是也可假设为常数，但是文献［17］和［18］在对方腔流的数值计算中采用高阶的关于坐标的多项式模拟顶部的速度分布，而本文采用的是恒定的速度边界。

2.3 计算效率

求解格子Boltzmann模型主要包括初始化变量，粒子的碰撞，迁移和边界处理，以及每隔1 000个迭代步判断结果是否收敛三个步骤。为了改进LBM求解三维问题的计算效率，采用GPU-CUDA并行计算技术加速D3Q19 MRT-SMAG模型，具体的实现过程参见Li等的工作［28］。本文使用的支持CUDA计算机的具体配置为：AMD Phenom II 1100T 3.3GHz处理器，16GB内存，以及具有1.5GB显存的GTX580显卡。软件环境为Windows XP操作系统，CUDA Toolkit版本为3.2，GPU驱动263.06。基于GPU的MRT-SMAG模型与仅用CPU进行计算的效率进行比较可知相对于仅用CPU处理，应用GPU对模型进行加速可以得到较高的计算比（见表1），这对应用LBM处理实际问题具有非常重要的意义。

3 结论

本文利用D3Q19 MRT-SMAG模型对三维方腔内的流场进行了数值模拟研究，并采用最新的基于GPU的CUDA并行技术以提高MRT-SMAG模型在模拟三维问题时的计算效率，将模型计算结果与实验和其他数值结果进行比较验证。可以得到如下结论：

（1）为了验证模型的有效性，采用本文所建立的MRT-SMAG模型对初始阶段的三维方腔流进行模拟，计算结果很好的符合实验数据表明MRT-SMAG可以对三维初始流场进行精确的预测。

（2）为了评估LBM模拟湍流的能力，采用MRT-SMAG模型求解了Re=10 000的三维方腔流，计算结果与已有数据符合较好，通过对比宽高比K=1和K=3的时间平均速度表明宽高比越大，边墙的摩擦作用对腔内湍流脉动的掺混能力越弱。另外通过分析表征湍流脉动强度的二阶统计量可知腔内湍流是沿墙边界产生，下游边界（x=L）附近的湍流脉动要比上游边界（x=0）附近的强烈，而底部强边界（z=0）附近的湍流脉动强度是最大的。

（3）利用GPU-CUDA并行技术对MRT-SMAG模型进行计算加速可以得到较高的加速比，约为200倍。

综上所述，MRT-SMAG模型能够较好的对湍流进行求解，而GPU-CUDA技术又能很好的加速模型的计算效率，因此将两者结合并引入可处理自由表面流动的VOF技术以建立数值波浪水槽具有非常好的应用前景和科研价值。

［1］Succi S.The lattice Boltzmann equation for fluid dynamics and beyond［M］.New York：Oxford University Press，2001.

［2］Kang S K，Hassan Y A.The effect of lattice models within the lattice Boltzmann method in the simulation of wall-bounded turbulent flows［J］.Journal of Computational Physics，2013，232：100-117.

［3］Thpmmes G，Becker J，Junk M，et al.A lattice Boltzmann method for immiscible multiphase flow simulations using the level set method［J］.Journal of Computational Physics，2009，228：1 139-1 156.

［4］Zhou J G.A lattice Boltzmann model for the shallow water equations［J］.Computer Methods in Applied Mechanics and Engineering，2002，191：3 527-3 539.

［5］Zhang JF，Zhang Q H.Lattice Boltzmann simulation of the flocculation process of cohesive sediment due to differential settling［J］.Continental Shelf Research，2011，31：S94-S105.

［6］Wang S，Xu C，Yuan P，et al.Hydrodynamic optimization of channeling device for hydro turbine based on lattice Boltzmann method［J］.Computers and Mathematics with Applications，2011，61：3 722-3 729.

［7］JANBEN C F，GRILLIST，KRAFCZYK M.Efficient simulations of long wave propagation and runup using a LBM approach on GPGPU hardware［C］//Chung J S.Proceedings of the Twenty-second International Offshore and Polar Engineering Conference. Rhodes：ISOPE，2012：17-22.

［8］JANSSEN C F，GRILLI S T，KRAFCZYK M.Modeling of wave breaking and wave-structure interactions by coupling of fully nonlinear potential flow and lattice-Boltzmann models［C］//Chung JS.Proceedings of the Twentieth International Offshore and Polar Engineering Conference.Beijing：ISOPE，2010：686-693.

［9］Burggraf O R.Analytical and numerical studies of the structure of steady separated flows［J］.Journal of Fluid Mechanics，1966，24：113-151.

［10］Pan F，Acrivos A.Steady flows in rectangular cavities［J］.Journal of Fluid Mechanics，1967，28（4）：643-655.

［11］Koseff JR，Street R L.Visualization studies of a shear driven three-dimensional recirculating flow［J］.ASME Journal of Fluids Engineering，1984，106：21-27.

［12］Koseff JR，Street R L.On end wall effects in a lid driven cavity flow［J］.ASME Journal of Fluids Engineering，1984，106：385-389.

［13］Koseff JR，Street R L.The lid-driven cavity flow：a synthesis of qualitative and quantitative observations［J］.ASME Journal of Fluids Engineering，1984，106：390-398.

［14］Prasad A K，Koseff JR.Reynolds number and end-wall effects on a lid driven cavity flow［J］.Physics of Fluids A，1989，1：208-218.

［15］Shankar PN，Deshpande M D.Fluid mechanics in the driven cavity［J］.Annual Review of Fluid Mechanics，2000，32：93-136.

［16］Guermond JL，Migeon C，Pineau G，etal.Start-up flows in a three-dimensional rectangular driven cavity of aspect ratio 1：1：2 at Re=1 000［J］.Journal of Fluid Mechanics，2002，450：169-199.

［17］Leriche E，Gavrilakis S.Direct numerical simulation of the flow in a lid-driven cubical cavity［J］.Physics of Fluids，2000，12（6）：1 363-1 376.

［18］Bouffanais R，Deville M O.Large eddy simulation of the flow in a lid-driven cubical cavity［J］.Physics of Fluids，2007，19：055108.

［19］D′Humières D，Ginzburg I，Krafczyk M，et al.Multiple-relaxation-time lattice Boltzmann models in three dimensions［J］. Philosophical Transactions of the Royal Society A，2002，360：437-451.

［20］De S，Nagendra K，Lakshmisha K N.Simulation of laminar flow in a three dimensional lid driven cavity by lattice Boltzmann method［J］.International Journal of NumericalMethods for Heat&Fluid Flow，2009，19（6）：790-815.

［21］Premnath K N，Pattison M J，Banerjee S.Generalized lattice Boltzmann equation with forcing term for computation of wallbounded turbulent flows［J］.Physical Review E，2009，79：026703.

［22］Smagorinsky J.General circulation experiments with the primitive equations：I.The basic equations［J］.Monthly Weather Review，1963，91：99-164.

［23］Kirk D B，Hwu WW.Programmingmassively parallel processors：a hands-on approach［M］.San Francisco：Morgan Kaufmann，2010.

［24］Pope S.Turbulent flows［M］.New York：Cambridge University Press，2000.

［25］Hou S L，Sterling J，Chen S Y，et al.A lattice Boltzmann subgrid model for high Reynolds number flows［J］.In：Pattern Formation and Lattice Gas Automata，in：A.T.Lawniczak，R.Kapral（Eds.），Fields Inst Comm，vol.6，AMS，Providence，1996，6：151-166.

［26］Krafczyk M，Tölke J，Luo L S.Large-eddy simulations with a multiple-relaxation-time LBE model［J］.International Journal of Modern Physics B，2003，17：33-39.

［27］Lätt J.Hydrodynamic limitof lattice Boltzmann equation［D］.Switzerland：University of Geneva，2007.

［28］LiCG，Maa JPY，Kang H G.Solving generalized lattice Boltzmannmodel for 3-D cavity flows using CUDA-GPU［J］.Science China Physics，Mechanics&Astronomy，2012，10：1 894-1 904.

Simulation of three-dimensional cavity flow with high Reynolds num ber using lattice Boltzmann equation

KANG Hai-gui,LICheng-gong
(State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116023,China)

In order to analyze the features of the flow field in the cavity and access the capability of simulating turbulent flow of the lattice Boltzmann method(LBM),the multiple relaxation time lattice Boltzmann equation(MRT-LBE)with the standard Smagorinsky eddy viscosity model was used for the simulation of three dimensional cavity flow with high Reynolds number(Re=10 000).The time averaged variables,e.g.,the velocity, root mean square velocity,Reynolds stress,and stream line contours on the symmetry plane(y=W/2)were computed.By the comparisons of the present results with the experimental and numerical results,it shows that the changing of the turbulent flow field in the cavity driven by the shear plane could be performed accurately by the MRT-LBE with the Smagorinsky model.In additional,the compute unified device architecture(CUDA)parallel technique on graphic processor unit(GPU)was introduced into MRT-LBE with the Smagorinsky model for improving the computational efficiency,up to 200 times.

multiple relaxation time lattice Boltzmann equation;standard Smagorinsky eddy viscosity model;high Reynolds number three dimensional cavity flow;GPU;CUDA

2012-04-07；

2013-05-12

国家自然科学基金（50679008）

康海贵（1945-），男，辽宁省大连人，教授，主要从事海岸工程方面的研究。

Biography：KANG Hai-gui（1945-），male，professor.

TV 131；O 242.1

A

1005-8443（2013）05-0453-08